1#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
2#if !defined(NWAD_PRINT)
3C> \ingroup nwxc
4C> @{
5C>
6C> \file nwxc_x_sogga.F
7C> The SOGGA, SOGGA11 and SOGGA-X exchange functionals
8C>
9C> @}
10#endif
11#endif
12C>
13C> \ingroup nwxc_priv
14C> @{
15C>
16C> \brief Evaluate the SOGGA, SOGGA11 and SOGGA11-X exchange functionals
17C>
18C> The SOGGA, SOGGA11 and SOGGA11-X functionals are GGA functionals that
19C> are formulated and optimized such that the coefficient of the
20C> gradient correction term matches that of the appropriate expansion ofC> the electronic energy [1,2,3].
21C>
22C> ### References ###
23C>
24C> [1] Y. Zhao, D.G. Truhlar,
25C>     "Construction of a generalized gradient approximation by
26C>     restoring the density-gradient expansion and enforcing a tight
27C>     Lieb-Oxford bound", J. Chem. Phys. <b>128</b> (2008) 184109,
28C>     DOI:
29C>     <a href="https://doi.org/10.1063/1.2912068">
30C>     10.1063/1.2912068</a>.
31C>
32C> [2] R. Peverati, Y. Zhao, D.G. Truhlar,
33C>     "Generalized gradient approximation that recovers the
34C>     second-order density-gradient expansion with optimized
35C>     across-the-board performance", J. Phys. Chem. Lett. <b>2</b>
36C>     (2011) 1991-1997, DOI:
37C>     <a href="https://doi.org/10.1021/jz200616w">
38C>     10.1021/jz200616w</a>.
39C>
40C> [3] R. Peverati, D.G. Truhlar,
41C>     "Communication: A global hybrid generalized gradient
42C>     approximation to the exchange-correlation functional that
43C>     satisfies the second-order density-gradient constraint and has
44C>     broad applicability in chemistry", J. Chem. Phys. <b>135</b>
45C>     (2011) 191102, DOI:
46C>     <a href="https://doi.org/10.1063/1.3663871">
47C>     10.1063/1.3663871</a>.
48C>
49#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
50#if defined(NWAD_PRINT)
51      Subroutine nwxc_x_sogga_p(param, tol_rho, ipol, nq, wght, rho,
52     &                          rgamma, ffunc)
53#else
54      Subroutine nwxc_x_sogga(param, tol_rho, ipol, nq, wght, rho,
55     &                        rgamma, ffunc)
56#endif
57#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV)
58      Subroutine nwxc_x_sogga_d2(param, tol_rho, ipol, nq, wght, rho,
59     &                           rgamma, ffunc)
60#else
61      Subroutine nwxc_x_sogga_d3(param, tol_rho, ipol, nq, wght, rho,
62     &                           rgamma, ffunc)
63#endif
64c
65c$Id$
66c
67c**********************************************************************c
68c                                                                      c
69c  SOGGA11X evaluates the exchange part of the SOGGA, SOGGA11          c
70c  and SOGGA11-X functionals on the grid.                              c
71c                                                                      c
72c     a) Zhao and Truhlar, J.Chem.Phys., 128, 184109 (2008)            c
73c     b) Peverati, Zhao and Truhlar, J.Phys.Chem.Lett, 2, 1991 (2011)  c
74c     c) Peverati and Truhlar, J.Chem.Phys, 135, 191102 (2011)         c
75c                                                                      c
76c      ijzy = 1 - SOGGA functional (a) - it requires PBE correlation   c
77c      ijzy = 2 - SOGGA11 functional (b)                               c
78c      ijzy = 3 - SOGGA11-X functional (c)                             c
79c                                                                      c
80c Coded by Roberto Peverati (12/11)                                    c
81c                                                                      c
82c**********************************************************************c
83c
84#include "nwad.fh"
85c
86      implicit none
87#include "nwxc_param.fh"
88c
89c     Input and other parameters
90c
91#if defined(NWAD_PRINT)
92#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
93      type(nwad_dble)::param(*)!< [Input] Parameters of the functional,
94      type(nwad_dble)::CxA0,CxA1,CxA2,CxA3,CxA4,CxA5
95      type(nwad_dble)::CxB0,CxB1,CxB2,CxB3,CxB4,CxB5
96      type(nwad_dble)::FA0
97      type(nwad_dble)::FB0
98#else
99      double precision param(*)!< [Input] Parameters of the functional,
100      double precision CxA0,CxA1,CxA2,CxA3,CxA4,CxA5
101      double precision CxB0,CxB1,CxB2,CxB3,CxB4,CxB5
102      double precision FA0
103      double precision FB0
104#endif
105#else
106      double precision param(*)!< [Input] Parameters of the functional,
107                               !< see Table 1 of [1] and Table 1 of [2].
108                               !< - param(1): \f$ a_0 \f$
109                               !< - param(2): \f$ a_1 \f$
110                               !< - param(3): \f$ a_2 \f$
111                               !< - param(4): \f$ a_3 \f$
112                               !< - param(5): \f$ a_4 \f$
113                               !< - param(6): \f$ a_5 \f$
114                               !< - param(7): \f$ b_0 \f$
115                               !< - param(8): \f$ b_1 \f$
116                               !< - param(9): \f$ b_2 \f$
117                               !< - param(10): \f$ b_3 \f$
118                               !< - param(11): \f$ b_4 \f$
119                               !< - param(12): \f$ b_5 \f$
120      double precision CxA0,CxA1,CxA2,CxA3,CxA4,CxA5
121      double precision CxB0,CxB1,CxB2,CxB3,CxB4,CxB5
122      double precision FA0
123      double precision FB0
124#endif
125      double precision tol_rho !< [Input] The lower limit on the density
126      integer ipol             !< [Input] The number of spin channels
127      integer nq               !< [Input] The number of points
128      double precision wght    !< [Input] The weight of the functional
129c
130c     Charge Density
131c
132      type(nwad_dble)::rho(nq,*)    !< [Input] The density
133c
134c     Charge Density Gradient
135c
136      type(nwad_dble)::rgamma(nq,*) !< [Input] The norm of the density gradients
137c
138c     Sampling Matrices for the XC Potential
139c
140      type(nwad_dble)::ffunc(nq)    !< [Output] The value of the functional
141c     double precision Amat(nq,*)   !< [Output] The derivative wrt rho
142c     double precision Cmat(nq,*)   !< [Output] The derivative wrt rgamma
143c
144      double precision pi
145c
146c     Intermediate derivative results, etc.
147c
148      integer n
149c
150      type(nwad_dble)::rho43, rho13, rhoo
151c
152      double precision AS, ASO, AX, DELOCDR
153      double precision DFA1DG, DFA1DR, DFA1DY
154      double precision DFA2DG, DFA2DR, DFA2DY
155      double precision DFA3DG, DFA3DR, DFA3DY
156      double precision DFA4DG, DFA4DR, DFA4DY
157      double precision DFA5DG, DFA5DR, DFA5DY
158      double precision DFB1DG, DFB1DR, DFB1DY
159      double precision DFB2DG, DFB2DR, DFB2DY
160      double precision DFB3DG, DFB3DR, DFB3DY
161      double precision DFB4DG, DFB4DR, DFB4DY
162      double precision DFB5DG, DFB5DR, DFB5DY
163      double precision DFEXPDPON, DFFRACDPON, DFGGAXDG, DFGGAXDR
164      double precision DYDG, DYDR, DTOL
165      double precision MU
166      type(nwad_dble)::FA1, FA2, FA3, FA4, FA5
167      type(nwad_dble)::FB1, FB2, FB3, FB4, FB5
168c     type(nwad_dble)::Gam12, Gam, X, FEXP, FFRAC, FGGAX, S, Y, ELOC
169      type(nwad_dble)::Gam, X, FEXP, FFRAC, FGGAX, Y, ELOC
170      type(nwad_dble)::PON
171c
172      double precision f1,f2,f3,f4,f5,f8
173      double precision F1o3,F4o3,F48
174      parameter( F1=1.0D+00,  F2=2.0D+00,  F3=3.0D+00,
175     $           F4=4.0D+00,  F5=5.0D+00,  F8=8.0D+00,
176     $           F48=48.0D+00)
177c
178      pi=acos(-1d0)
179c
180c     if (ijzy.eq.1) then
181c SOGGA
182c      CxA0 = 0.5d0
183c      CxA1 = 0.276d0
184c      CxA2 = 0d0
185c      CxA3 = 0d0
186c      CxA4 = 0d0
187c      CxA5 = 0d0
188c      CxB0 = 0.5d0
189c      CxB1 = 0.276d0
190c      CxB2 = 0d0
191c      CxB3 = 0d0
192c      CxB4 = 0d0
193c      CxB5 = 0d0
194c     elseif (ijzy.eq.2) then
195c SOGGA11
196c      CxA0 =  0.50000d0
197c      CxA1 = -2.95535d0
198c      CxA2 =  15.7974d0
199c      CxA3 = -91.1804d0
200c      CxA4 =  96.2030d0
201c      CxA5 =  0.18683d0
202c      CxB0 =  0.50000d0
203c      CxB1 =  3.50743d0
204c      CxB2 = -12.9523d0
205c      CxB3 =  49.7870d0
206c      CxB4 = -33.2545d0
207c      CxB5 = -11.1396d0
208c     elseif (ijzy.eq.3) then
209c SOGGA11-X
210c      CxA0 =  2.99250d-01
211c      CxA1 =  3.21638d+00
212c      CxA2 = -3.55605d+00
213c      CxA3 =  7.65852d+00
214c      CxA4 = -1.12830d+01
215c      CxA5 =  5.25813d+00
216c      CxB0 =  2.99250d-01
217c      CxB1 = -2.88595d+00
218c      CxB2 =  3.23617d+00
219c      CxB3 = -2.45393d+00
220c      CxB4 = -3.75495d+00
221c      CxB5 =  3.96613d+00
222c     endif
223      CxA0 = param(1)
224      CxA1 = param(2)
225      CxA2 = param(3)
226      CxA3 = param(4)
227      CxA4 = param(5)
228      CxA5 = param(6)
229      CxB0 = param(7)
230      CxB1 = param(8)
231      CxB2 = param(9)
232      CxB3 = param(10)
233      CxB4 = param(11)
234      CxB5 = param(12)
235c
236      DTol = tol_rho
237      F1o3 = F1/F3
238      F4o3 = F4/F3
239      Pi   = ACos(-F1)
240      AsO  = (F48*PI*PI)**F1o3
241      As   = F1/AsO
242      Ax   = -(F3/F2) * (F4o3*Pi)**(-F1o3)
243      mu = 0.2236536053d0
244c
245      if (ipol.eq.1 )then
246c
247c        ======> SPIN-RESTRICTED <======
248c                     or
249c                SPIN-UNPOLARIZED
250c
251c
252         do 10 n = 1, nq
253            if (rho(n,R_T).lt.DTol) goto 10
254            rhoo = rho(n,R_T)/F2
255            rho43 = rhoo**F4o3
256            rho13 = rho43/rhoo
257            Gam = rgamma(n,G_TT)/F4
258c           Gam =(delrho(n,1,1)*delrho(n,1,1) +
259c    &              delrho(n,2,1)*delrho(n,2,1) +
260c    &              delrho(n,3,1)*delrho(n,3,1))/F4
261c           Gam12 = sqrt(Gam)
262c           if(gam12.lt.dtol) goto 10
263c
264            Eloc = Ax*Rho43
265c           x = Gam12/Rho43
266c           s = As*x
267c           y = s*s
268            y = As*As*Gam/(Rho43*Rho43)
269            PON = mu*y
270            Ffrac = F1-F1/(F1+PON)
271            Fexp  = F1-exp(-PON)
272            fa0 = CxA0
273            fa1 = CxA1 *Ffrac
274            fa2 = CxA2 *Ffrac**F2
275            fa3 = CxA3 *Ffrac**F3
276            fa4 = CxA4 *Ffrac**F4
277            fa5 = CxA5 *Ffrac**F5
278            fb0 = CxB0
279            fb1 = CxB1 *Fexp
280            fb2 = CxB2 *Fexp**F2
281            fb3 = CxB3 *Fexp**F3
282            fb4 = CxB4 *Fexp**F4
283            fb5 = CxB5 *Fexp**F5
284c
285            Fggax = fa0+fa1+fa2+fa3+fa4+fa5 +
286     $              fb0+fb1+fb2+fb3+fb4+fb5
287C
288C     1st derivatives.
289C
290
291c           dElocdR=Ax*F4o3*Rho13
292c           dydR = -(F8/F3)*y/Rhoo
293c           dydG   = y/Gam
294c           dFfracdPON = F1/((F1+PON)**F2)
295c           dFexpdPON  = exp(-PON)
296c           dfa1dy = CxA1 *mu*dFfracdPON
297c           dfa2dy = CxA2 *mu*F2*Ffrac*dFfracdPON
298c           dfa3dy = CxA3 *mu*(F3 *Ffrac**F2)*dFfracdPON
299c           dfa4dy = CxA4 *mu*(F4 *Ffrac**F3)*dFfracdPON
300c           dfa5dy = CxA5 *mu*(F5 *Ffrac**F4)*dFfracdPON
301c           dfa1dR = dfa1dy *dydR
302c           dfa2dR = dfa2dy *dydR
303c           dfa3dR = dfa3dy *dydR
304c           dfa4dR = dfa4dy *dydR
305c           dfa5dR = dfa5dy *dydR
306c           dfa1dG = dfa1dy *dydG
307c           dfa2dG = dfa2dy *dydG
308c           dfa3dG = dfa3dy *dydG
309c           dfa4dG = dfa4dy *dydG
310c           dfa5dG = dfa5dy *dydG
311c           dfb1dy = CxB1 *mu*dFexpdPON
312c           dfb2dy = CxB2 *mu*F2*Fexp*dFexpdPON
313c           dfb3dy = CxB3 *mu*(F3 *Fexp**F2)*dFexpdPON
314c           dfb4dy = CxB4 *mu*(F4 *Fexp**F3)*dFexpdPON
315c           dfb5dy = CxB5 *mu*(F5 *Fexp**F4)*dFexpdPON
316c           dfb1dR = dfb1dy *dydR
317c           dfb2dR = dfb2dy *dydR
318c           dfb3dR = dfb3dy *dydR
319c           dfb4dR = dfb4dy *dydR
320c           dfb5dR = dfb5dy *dydR
321c           dfb1dG = dfb1dy *dydG
322c           dfb2dG = dfb2dy *dydG
323c           dfb3dG = dfb3dy *dydG
324c           dfb4dG = dfb4dy *dydG
325c           dfb5dG = dfb5dy *dydG
326c
327c           dFggaxdR = dfa1dR+dfa2dR+dfa3dR+dfa4dR+dfa5dR +
328c    $                 dfb1dR+dfb2dR+dfb3dR+dfb4dR+dfb5dR
329c
330c           dFggaxdG = dfa1dG+dfa2dG+dfa3dG+dfa4dG+dfa5dG +
331c    $                 dfb1dG+dfb2dG+dfb3dG+dfb4dG+dfb5dG
332c
333c          Ex = Ex +F2*(Eloc*Fggax)*qwght(n)
334           ffunc(n)=ffunc(n)+F2*(Eloc*Fggax)*wght
335c          Amat(n,D1_RA)=Amat(n,D1_RA)
336c    $                  +(dElocdR*Fggax+Eloc*dFggaxdR)*wght
337c          Cmat(n,D1_GAA)=Cmat(n,D1_GAA)+Eloc*dFggaxdG*wght
338c
33910      continue
340c
341c UUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUnrestricted
342      else  ! ipol=2
343c
344c        ======> SPIN-UNRESTRICTED <======
345
346c
347c  use spin density functional theory ie n-->2n
348c                               Ex=(1/2)Ex[2*alpha] + (1/2)Ex[2*beta]
349c
350c     Alpha            ALPHA               ALPHA
351c
352         do 20 n = 1, nq
353           if (rho(n,R_A)+rho(n,R_B).lt.DTol) goto 20
354           if (rho(n,R_A).lt.DTol) goto 25
355           rhoo  = rho(n,R_A)
356           rho43 = rhoo**F4o3
357           rho13 = rho43/rhoo
358           Gam = rgamma(n,G_AA)
359c          Gam =(delrho(n,1,1)*delrho(n,1,1) +
360c    &           delrho(n,2,1)*delrho(n,2,1) +
361c    &           delrho(n,3,1)*delrho(n,3,1))
362c          Gam12 = sqrt(Gam)
363c          if(gam12.lt.dtol) goto 25
364c
365           Eloc = Ax*Rho43
366c          x = Gam12/Rho43
367c          s = As*x
368c          y = s*s
369           y = As*As*Gam/(Rho43*Rho43)
370           PON = mu*y
371           Ffrac = F1-F1/(F1+PON)
372           Fexp  = F1-exp(-PON)
373           fa0 = CxA0
374           fa1 = CxA1 *Ffrac
375           fa2 = CxA2 *Ffrac**F2
376           fa3 = CxA3 *Ffrac**F3
377           fa4 = CxA4 *Ffrac**F4
378           fa5 = CxA5 *Ffrac**F5
379           fb0 = CxB0
380           fb1 = CxB1 *Fexp
381           fb2 = CxB2 *Fexp**F2
382           fb3 = CxB3 *Fexp**F3
383           fb4 = CxB4 *Fexp**F4
384           fb5 = CxB5 *Fexp**F5
385c
386           Fggax = fa0+fa1+fa2+fa3+fa4+fa5 +
387     $             fb0+fb1+fb2+fb3+fb4+fb5
388C
389C     1st derivatives.
390C
391c          dElocdR=Ax*F4o3*Rho13
392c          dydR = -(F8/F3)*y/Rhoo
393c          dydG   = y/Gam
394c          dFfracdPON = F1/((F1+PON)**F2)
395c          dFexpdPON  = exp(-PON)
396c          dfa1dy = CxA1 *mu*dFfracdPON
397c          dfa2dy = CxA2 *mu*F2*Ffrac*dFfracdPON
398c          dfa3dy = CxA3 *mu*(F3 *Ffrac**F2)*dFfracdPON
399c          dfa4dy = CxA4 *mu*(F4 *Ffrac**F3)*dFfracdPON
400c          dfa5dy = CxA5 *mu*(F5 *Ffrac**F4)*dFfracdPON
401c          dfa1dR = dfa1dy *dydR
402c          dfa2dR = dfa2dy *dydR
403c          dfa3dR = dfa3dy *dydR
404c          dfa4dR = dfa4dy *dydR
405c          dfa5dR = dfa5dy *dydR
406c          dfa1dG = dfa1dy *dydG
407c          dfa2dG = dfa2dy *dydG
408c          dfa3dG = dfa3dy *dydG
409c          dfa4dG = dfa4dy *dydG
410c          dfa5dG = dfa5dy *dydG
411c          dfb1dy = CxB1 *mu*dFexpdPON
412c          dfb2dy = CxB2 *mu*F2*Fexp*dFexpdPON
413c          dfb3dy = CxB3 *mu*(F3 *Fexp**F2)*dFexpdPON
414c          dfb4dy = CxB4 *mu*(F4 *Fexp**F3)*dFexpdPON
415c          dfb5dy = CxB5 *mu*(F5 *Fexp**F4)*dFexpdPON
416c          dfb1dR = dfb1dy *dydR
417c          dfb2dR = dfb2dy *dydR
418c          dfb3dR = dfb3dy *dydR
419c          dfb4dR = dfb4dy *dydR
420c          dfb5dR = dfb5dy *dydR
421c          dfb1dG = dfb1dy *dydG
422c          dfb2dG = dfb2dy *dydG
423c          dfb3dG = dfb3dy *dydG
424c          dfb4dG = dfb4dy *dydG
425c          dfb5dG = dfb5dy *dydG
426c
427c          dFggaxdR = dfa1dR+dfa2dR+dfa3dR+dfa4dR+dfa5dR +
428c    $                dfb1dR+dfb2dR+dfb3dR+dfb4dR+dfb5dR
429c
430c          dFggaxdG = dfa1dG+dfa2dG+dfa3dG+dfa4dG+dfa5dG +
431c    $                dfb1dG+dfb2dG+dfb3dG+dfb4dG+dfb5dG
432c
433c          Ex = Ex + (Eloc*Fggax)*qwght(n)
434           ffunc(n)=ffunc(n)+ Eloc*Fggax*wght
435c          Amat(n,D1_RA) = Amat(n,D1_RA) + dElocdR*Fggax*wght
436c    $                                   + Eloc*dFggaxdR*wght
437c          Cmat(n,D1_GAA)=  Cmat(n,D1_GAA)   + Eloc*dFggaxdG*wght
438c
43925         continue
440c
441c     Beta               BETA           BETA
442c
443           if (rho(n,R_B).lt.DTol) goto 20
444           rhoo  = rho(n,R_B)
445           rho43 = rhoo**F4o3
446           rho13 = rho43/rhoo
447c
448           Gam = rgamma(n,G_BB)
449c          Gam =(delrho(n,1,2)*delrho(n,1,2) +
450c    &           delrho(n,2,2)*delrho(n,2,2) +
451c    &           delrho(n,3,2)*delrho(n,3,2))
452c          Gam12 = sqrt(Gam)
453c          if(gam12.lt.dtol) goto 20
454c
455           Eloc = Ax*Rho43
456c          x = Gam12/Rho43
457c          s = As*x
458c          y = s*s
459           y = As*As*Gam/(Rho43*Rho43)
460           PON = mu*y
461           Ffrac = F1-F1/(F1+PON)
462           Fexp  = F1-exp(-PON)
463           fa0 = CxA0
464           fa1 = CxA1 *Ffrac
465           fa2 = CxA2 *Ffrac**F2
466           fa3 = CxA3 *Ffrac**F3
467           fa4 = CxA4 *Ffrac**F4
468           fa5 = CxA5 *Ffrac**F5
469           fb0 = CxB0
470           fb1 = CxB1 *Fexp
471           fb2 = CxB2 *Fexp**F2
472           fb3 = CxB3 *Fexp**F3
473           fb4 = CxB4 *Fexp**F4
474           fb5 = CxB5 *Fexp**F5
475c
476           Fggax = fa0+fa1+fa2+fa3+fa4+fa5 +
477     $             fb0+fb1+fb2+fb3+fb4+fb5
478C
479C     1st derivatives.
480C
481
482c          dElocdR=Ax*F4o3*Rho13
483c          dydR = -(F8/F3)*y/Rhoo
484c          dydG   = y/Gam
485c          dFfracdPON = F1/((F1+PON)**F2)
486c          dFexpdPON  = exp(-PON)
487c          dfa1dy = CxA1 *mu*dFfracdPON
488c          dfa2dy = CxA2 *mu*F2*Ffrac*dFfracdPON
489c          dfa3dy = CxA3 *mu*(F3 *Ffrac**F2)*dFfracdPON
490c          dfa4dy = CxA4 *mu*(F4 *Ffrac**F3)*dFfracdPON
491c          dfa5dy = CxA5 *mu*(F5 *Ffrac**F4)*dFfracdPON
492c          dfa1dR = dfa1dy *dydR
493c          dfa2dR = dfa2dy *dydR
494c          dfa3dR = dfa3dy *dydR
495c          dfa4dR = dfa4dy *dydR
496c          dfa5dR = dfa5dy *dydR
497c          dfa1dG = dfa1dy *dydG
498c          dfa2dG = dfa2dy *dydG
499c          dfa3dG = dfa3dy *dydG
500c          dfa4dG = dfa4dy *dydG
501c          dfa5dG = dfa5dy *dydG
502c          dfb1dy = CxB1 *mu*dFexpdPON
503c          dfb2dy = CxB2 *mu*F2*Fexp*dFexpdPON
504c          dfb3dy = CxB3 *mu*(F3 *Fexp**F2)*dFexpdPON
505c          dfb4dy = CxB4 *mu*(F4 *Fexp**F3)*dFexpdPON
506c          dfb5dy = CxB5 *mu*(F5 *Fexp**F4)*dFexpdPON
507c          dfb1dR = dfb1dy *dydR
508c          dfb2dR = dfb2dy *dydR
509c          dfb3dR = dfb3dy *dydR
510c          dfb4dR = dfb4dy *dydR
511c          dfb5dR = dfb5dy *dydR
512c          dfb1dG = dfb1dy *dydG
513c          dfb2dG = dfb2dy *dydG
514c          dfb3dG = dfb3dy *dydG
515c          dfb4dG = dfb4dy *dydG
516c          dfb5dG = dfb5dy *dydG
517c
518c          dFggaxdR = dfa1dR+dfa2dR+dfa3dR+dfa4dR+dfa5dR +
519c    $                dfb1dR+dfb2dR+dfb3dR+dfb4dR+dfb5dR
520c
521c          dFggaxdG = dfa1dG+dfa2dG+dfa3dG+dfa4dG+dfa5dG +
522c    $                dfb1dG+dfb2dG+dfb3dG+dfb4dG+dfb5dG
523c
524c          Ex = Ex + (Eloc*Fggax)*qwght(n)
525           ffunc(n)=ffunc(n)+ Eloc*Fggax*wght
526
527c          Amat(n,D1_RB) = Amat(n,D1_RB) + dElocdR*Fggax*wght
528c    $                                   + Eloc*dFggaxdR*wght
529c
530c          Cmat(n,D1_GBB)=  Cmat(n,D1_GBB)  + Eloc*dFggaxdG*wght
531c
53220      continue
533      endif
534      return
535      end
536#ifndef NWAD_PRINT
537#define NWAD_PRINT
538c
539c     compile again for 2nd derivatives
540c
541#include "nwxc_x_sogga.F"
542#endif
543#ifndef SECOND_DERIV
544#define SECOND_DERIV
545c
546c     compile again for 2nd derivatives
547c
548#include "nwxc_x_sogga.F"
549#endif
550#ifndef THIRD_DERIV
551#define THIRD_DERIV
552c
553c     compile again for 3rd derivatives
554c
555#include "nwxc_x_sogga.F"
556#endif
557#undef NWAD_PRINT
558C>
559C> @}
560