1#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
2#if !defined(NWAD_PRINT)
3C> \ingroup nwxc
4C> @{
5C>
6C> \file nwxc_c_vs98.F
7C> Implementation of the VS98 correlation functional
8C>
9C> @}
10#endif
11#endif
12C>
13C> \ingroup nwxc_priv
14C> @{
15C>
16C> \brief The Voorhis and Scuseria correlation functional
17C>
18C> The VS98 functional [1,2] is a meta-GGA. This routine implements
19C> the correlation component.
20C>
21C> Due to the form of the meta-GGAs we need to screen on the kinetic
22C> energy density to ensure that LDA will be obtained when the kinetic
23C> energy density goes to zero [3].
24C>
25C> ### References ###
26C>
27C> [1] T van Voorhis, GE Scuseria,
28C>     "A novel form for the exchange-correlation energy functional",
29C>     J.Chem.Phys. <b>109</b>, 400-410 (1998), DOI:
30C>     <a href="https://doi.org/10.1063/1.476577">
31C>     10.1063/1.476577</a>.
32C>
33C> [2] T van Voorhis, GE Scuseria,
34C>     Erratum: "A novel form for the exchange-correlation energy
35C>     functional",
36C>     J.Chem.Phys. <b>129</b>, 219901-219901 (2008), DOI:
37C>     <a href="https://doi.org/10.1063/1.3005348">
38C>     10.1063/1.3005348</a>.
39C>
40C> [3] J. Gr&auml;fenstein, D. Izotov, D. Cremer,
41C>     "Avoiding singularity problems associated with meta-GGA exchange
42C>     and correlation functionals containing the kinetic energy
43C>     density", J. Chem. Phys. <b>127</b>, 214103 (2007), DOI:
44C>     <a href="https://doi.org/10.1063/1.2800011">
45C>     10.1063/1.2800011</a>.
46C>
47c    VSXC correlation functional
48c           META GGA
49C         utilizes ingredients:
50c                              rho   -  density
51c                              delrho - gradient of density
52c                              tau (tauN)- K.S kinetic energy density
53c                              ijzy - 1  VS98
54c                              ijzy - 2  M06-L
55c                              ijzy - 3  M06-HF
56c                              ijzy - 4  M06
57c                              ijzy - 5  M06-2X
58c
59#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
60#if defined(NWAD_PRINT)
61      Subroutine nwxc_c_vs98_p(param, tol_rho, ipol, nq, wght, rho,
62     &                         rgamma, tau, func)
63#else
64      Subroutine nwxc_c_vs98(param, tol_rho, ipol, nq, wght, rho,
65     &                       rgamma, tau, func)
66#endif
67#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV)
68      Subroutine nwxc_c_vs98_d2(param, tol_rho, ipol, nq, wght, rho,
69     &                          rgamma, tau, func)
70#else
71      Subroutine nwxc_c_vs98_d3(param, tol_rho, ipol, nq, wght, rho,
72     &                          rgamma, tau, func)
73#endif
74c
75c$Id$
76c
77c  Reference
78c   [a] T. V. Voorhis and G. E. Scuseria, J. Chem. Phys. 109, 400 (1998).
79c   [b] Y. Zhao and D. G. Truhlar, J. Chem. Phys. 125, 194101 (2006).
80
81c
82#include "nwad.fh"
83c
84      implicit none
85c
86#include "intf_nwxc_c_lsda.fh"
87#include "intf_nwxc_vs98ss.fh"
88c
89#include "nwxc_param.fh"
90c
91c     Input and other parameters
92c
93#if defined(NWAD_PRINT)
94#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
95      type(nwad_dble)::param(*) !< [Input] Parameters of functional (table 1)
96      type(nwad_dble)::r7, r8, r9, r10, r11, r12
97      type(nwad_dble)::r13, r14, r15, r16, r17, r18
98#else
99      double precision param(*) !< [Input] Parameters of functional (table 1)
100      double precision r7, r8, r9, r10, r11, r12
101      double precision r13, r14, r15, r16, r17, r18
102#endif
103#else
104      double precision param(*) !< [Input] Parameters of functional (table 1)
105                                !< - param(1): \f$ a_{\sigma\sigma'} \f$
106                                !< - param(2): \f$ b_{\sigma\sigma'} \f$
107                                !< - param(3): \f$ c_{\sigma\sigma'} \f$
108                                !< - param(4): \f$ d_{\sigma\sigma'} \f$
109                                !< - param(5): \f$ e_{\sigma\sigma'} \f$
110                                !< - param(6): \f$ f_{\sigma\sigma'} \f$
111                                !< - param(7): \f$ a_{\sigma\sigma} \f$
112                                !< - param(8): \f$ b_{\sigma\sigma} \f$
113                                !< - param(9): \f$ c_{\sigma\sigma} \f$
114                                !< - param(10): \f$ d_{\sigma\sigma} \f$
115                                !< - param(11): \f$ e_{\sigma\sigma} \f$
116                                !< - param(12): \f$ f_{\sigma\sigma} \f$
117      double precision r7, r8, r9, r10, r11, r12
118      double precision r13, r14, r15, r16, r17, r18
119#endif
120      double precision tol_rho !< [Input] The lower limit on the density
121      integer nq               !< [Input] The number of points
122      integer ipol             !< [Input] The number of spin channels
123      double precision wght    !< [Input] The weight of the functional
124c
125c     Charge Density
126c
127      type(nwad_dble)::rho(nq,*) !< [Input] The density
128c
129c     Charge Density Gradient Norm
130c
131      type(nwad_dble)::rgamma(nq,*) !< [Input] The density gradient norm
132c
133c     Kinetic Energy Density
134c
135      type(nwad_dble)::tau(nq,*) !< [Input] The kinetic energy density
136c
137c     Functional values
138c
139      type(nwad_dble)::func(*) !< [Output] The functional value
140c
141c     Sampling Matrices for the XC Potential
142c
143c     double precision Amat(nq,*) !< [Output] Derivative wrt density
144c     double precision Cmat(nq,*) !< [Output] Derivative wrt rgamma
145c     double precision Mmat(nq,*) !< [Output] Derivative wrt tau
146c
147c     Threshold parameters
148c
149      double precision DTol,F1, F2, F3, F4, gab, cf
150      Data F1/1.0d0/,F2/2.0d0/,
151     & F3/3.0d0/,F4/4.0d0/,gab/0.00304966d0/,
152     & cf/9.115599720d0/
153c
154c     Local
155c
156      integer n
157      double precision tauN
158
159c    call to the vs98css subroutine
160      type(nwad_dble)::PA,GAA,TA,FA,EUA,ZA,ChiA
161      double precision FPA,FGA,FTA,EUEGA,EUPA,ChiAP,ChiAG,
162     &                 ZAP,ZAT
163      type(nwad_dble)::PB,GBB,TB,FB,EUB,ZB,ChiB
164      double precision FPB,FGB,FTB,EUEGB,EUPB,ChiBP,ChiBG,
165     &                 ZBP,ZBT
166c
167      type(nwad_dble)::RS,Zeta,gcab,P,PotLC,EUEG,ZAB
168      double precision Pi, F43, F13, Pi34, F6,
169     &RSP,dZdA,dZdB,dLdS,dLdZ,d2LdSS,d2LdSZ,d2LdZZ
170      type(nwad_dble)::XAB, kab, xk, zk
171      double precision dgdx,dgdz,dgdPA,dgdGA,dgdTA,dgdPB,dgdGB,dgdTB
172      double precision EUEGPA,EUEGPB
173
174
175c
176c     ======> BOTH SPIN-RESTRICETED AND UNRESTRICTED <======
177c
178c     DTol=1.0d-7
179      dtol=tol_rho
180C     Parameters for VS98 / M06-L / M06-HF / M06 / M06-2X
181      r7=   param(1)
182      r8=   param(2)
183      r9=   param(3)
184      r10=  param(4)
185      r11=  param(5)
186      r12=  param(6)
187      r13=  param(7)
188      r14=  param(8)
189      r15=  param(9)
190      r16=  param(10)
191      r17=  param(11)
192      r18=  param(12)
193c
194      Pi = F4*ATan(F1)
195      F6=6.0d0
196      F43 = F4 / F3
197      Pi34 = F3 / (F4*Pi)
198      F13 = F1 / F3
199c
200      do 20 n = 1, nq
201        if (ipol.eq.1) then
202          if (rho(n,R_A).lt.DTol) goto 20
203        else
204          if (rho(n,R_A)+rho(n,R_B).lt.DTol) goto 20
205        endif
206        if (ipol.eq.1) then
207c
208c         get the density, gradient, and tau for the alpha spin
209c         from the total
210c
211          PA   = rho(n,R_A)/F2
212          PB   = 0.0d0
213          EUA  = 0.0d0
214          ChiA = 0.0d0
215          EUB  = 0.0d0
216          ChiB = 0.0d0
217          GAA  = rgamma(n,G_AA)/4.0d0
218c         In the bc95css subroutine, we use 2*TA as the tau, so we do
219c         not divide the tau by 2 here
220
221          TA = tau(n,T_A)
222          if(ta.lt.dtol) goto 20
223
224#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
225#if defined(NWAD_PRINT)
226          Call nwxc_vs98ss_p(tol_rho,PA,GAA,TA,FA,FPA,FGA,FTA,EUA,ZA,
227     &                     ChiA,EUPA,ChiAP,ChiAG,ZAP,ZAT,
228     &                     r13,r14,r15,r16,r17,r18)
229#else
230          Call nwxc_vs98ss(tol_rho,PA,GAA,TA,FA,FPA,FGA,FTA,EUA,ZA,
231     &                     ChiA,EUPA,ChiAP,ChiAG,ZAP,ZAT,
232     &                     r13,r14,r15,r16,r17,r18)
233#endif
234#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV)
235          Call nwxc_vs98ss_d2(tol_rho,PA,GAA,TA,FA,FPA,FGA,FTA,EUA,ZA,
236     &                     ChiA,EUPA,ChiAP,ChiAG,ZAP,ZAT,
237     &                     r13,r14,r15,r16,r17,r18)
238#else
239          Call nwxc_vs98ss_d3(tol_rho,PA,GAA,TA,FA,FPA,FGA,FTA,EUA,ZA,
240     &                     ChiA,EUPA,ChiAP,ChiAG,ZAP,ZAT,
241     &                     r13,r14,r15,r16,r17,r18)
242#endif
243          PB = PA
244          GBB = GAA
245          TB = TA
246          FB = FA
247          FPB = FPA
248          FGB = FGA
249          FTB = FTA
250          EUB = EUA
251          ZB = ZA
252          ChiB = ChiA
253          EUPB = EUPA
254          ChiBP = ChiAP
255          ChiBG = ChiAG
256          ZBP = ZAP
257          ZBT = ZAT
258
259          func(n)        = func(n)        + 2.0d0*FA*wght
260c         Amat(n,D1_RA)  = Amat(n,D1_RA)  + FPA*wght
261c         Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + FGA*wght
262c         Mmat(n,D1_TA)  = Mmat(n,D1_TA)  + FTA*wght
263
264
265c UUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUnrestricted
266      else  ! ipol=2
267c
268c        ======> SPIN-UNRESTRICTED <======
269c
270c
271c       alpha
272c
273
274         PA   = 0.0d0
275         GAA  = 0.0d0
276         TA   = 0.0d0
277         EUA  = 0.0d0
278         ChiA = 0.0d0
279         ZA   = 0.0d0
280         if (rho(n,R_A).le.0.5d0*DTol) go to 25
281         PA  = rho(n,R_A)
282c
283c        In the bc95css subroutine, we use 2*TA as the tau
284c
285         if (tau(n,T_A).le.0.5d0*DTol) go to 25
286         GAA = rgamma(n,G_AA)
287         TA = tau(n,T_A)*2.0d0
288
289#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
290#if defined(NWAD_PRINT)
291         Call nwxc_vs98ss_p(tol_rho,PA,GAA,TA,FA,FPA,FGA,FTA,EUA,ZA,
292     &                    ChiA,EUPA,ChiAP,ChiAG,ZAP,ZAT,
293     &                    r13,r14,r15,r16,r17,r18)
294#else
295         Call nwxc_vs98ss(tol_rho,PA,GAA,TA,FA,FPA,FGA,FTA,EUA,ZA,
296     &                    ChiA,EUPA,ChiAP,ChiAG,ZAP,ZAT,
297     &                    r13,r14,r15,r16,r17,r18)
298#endif
299#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV)
300         Call nwxc_vs98ss_d2(tol_rho,PA,GAA,TA,FA,FPA,FGA,FTA,EUA,ZA,
301     &                    ChiA,EUPA,ChiAP,ChiAG,ZAP,ZAT,
302     &                    r13,r14,r15,r16,r17,r18)
303#else
304         Call nwxc_vs98ss_d3(tol_rho,PA,GAA,TA,FA,FPA,FGA,FTA,EUA,ZA,
305     &                    ChiA,EUPA,ChiAP,ChiAG,ZAP,ZAT,
306     &                    r13,r14,r15,r16,r17,r18)
307#endif
308         func(n)        = func(n)        + FA*wght
309c        Amat(n,D1_RA)  = Amat(n,D1_RA)  + FPA*wght
310c        Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + FGA*wght
311c        2*0.5=1.0 for Mmat
312c        Mmat(n,D1_TA)  = Mmat(n,D1_TA)  + FTA*wght
313
314c
315c        In the vs98ss subroutine, we use 2*TA as the tau,
316c
317c
318c       Beta
319c
320 25      continue
321         PB   = 0.0d0
322         GBB  = 0.0d0
323         TB   = 0.0d0
324         EUB  = 0.0d0
325         ChiB = 0.0d0
326         ZB   = 0.0d0
327
328         if(rho(n,R_B).le.0.5d0*DTol) go to 30
329
330         PB = rho(n,R_B)
331
332         if(tau(n,T_B).le.0.5d0*DTol) go to 30
333         GBB = rgamma(n,G_BB)
334         TB = tau(n,T_B)*2.0d0
335
336#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
337#if defined(NWAD_PRINT)
338         Call nwxc_vs98ss_p(tol_rho,PB,GBB,TB,FB,FPB,FGB,FTB,EUB,ZB,
339     &                    ChiB,EUPB,ChiBP,ChiBG,ZBP,ZBT,
340     &                    r13,r14,r15,r16,r17,r18)
341#else
342         Call nwxc_vs98ss(tol_rho,PB,GBB,TB,FB,FPB,FGB,FTB,EUB,ZB,
343     &                    ChiB,EUPB,ChiBP,ChiBG,ZBP,ZBT,
344     &                    r13,r14,r15,r16,r17,r18)
345#endif
346#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV)
347         Call nwxc_vs98ss_d2(tol_rho,PB,GBB,TB,FB,FPB,FGB,FTB,EUB,ZB,
348     &                    ChiB,EUPB,ChiBP,ChiBG,ZBP,ZBT,
349     &                    r13,r14,r15,r16,r17,r18)
350#else
351         Call nwxc_vs98ss_d3(tol_rho,PB,GBB,TB,FB,FPB,FGB,FTB,EUB,ZB,
352     &                    ChiB,EUPB,ChiBP,ChiBG,ZBP,ZBT,
353     &                    r13,r14,r15,r16,r17,r18)
354#endif
355         func(n)        = func(n)        + FB*wght
356c        Amat(n,D1_RB)  = Amat(n,D1_RB)  + FPB*wght
357c        Cmat(n,D1_GBB) = Cmat(n,D1_GBB) + FGB*wght
358c        Mmat(n,D1_TB)  = Mmat(n,D1_TB)  + FTB*wght
359
360#if 0
361      write (0,'(A,3F20.6)') "BAmat Cmat Mmat",FPB,FGB,FTB
362#endif
363      endif
364 30   continue
365      P = PA + PB
366      If(PA.gt.DTol.or.PB.gt.DTol) then
367          RS = (Pi34/P) ** F13
368c         RSP = -RS/(F3*P)
369          Zeta = (PA-PB)/P
370c         dZdA = (F1-Zeta)/P
371c         dZdB = (-F1-Zeta)/P
372#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
373#if defined(NWAD_PRINT)
374          Call nwxc_c_lsda_p(tol_rho,RS,Zeta,PotLC,
375     &               dLdS,dLdZ,d2LdSS,d2LdSZ,d2LdZZ)
376#else
377          Call nwxc_c_lsda(tol_rho,RS,Zeta,PotLC,
378     &               dLdS,dLdZ,d2LdSS,d2LdSZ,d2LdZZ)
379#endif
380#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV)
381          Call nwxc_c_lsda_d2(tol_rho,RS,Zeta,PotLC,
382     &               dLdS,dLdZ,d2LdSS,d2LdSZ,d2LdZZ)
383#else
384          Call nwxc_c_lsda_d3(tol_rho,RS,Zeta,PotLC,
385     &               dLdS,dLdZ,d2LdSS,d2LdSZ,d2LdZZ)
386#endif
387          EUEG = P*PotLC - EUA - EUB
388          ZAB = ZA + ZB
389          XAB = ChiA+ChiB
390          kab = F1 + gab*(XAB+ZAB)
391          xk = XAB/kab
392          zk = ZAB/kab
393#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
394#if defined(NWAD_PRINT)
395          call nwxc_gvt4_p(gcab,dgdx,dgdz,xk,zk,kab,gab,gab,
396     &                   r7,r8,r9,r10,r11,r12)
397#else
398          call nwxc_gvt4(gcab,dgdx,dgdz,xk,zk,kab,gab,gab,
399     &                   r7,r8,r9,r10,r11,r12)
400#endif
401#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV)
402          call nwxc_gvt4_d2(gcab,dgdx,dgdz,xk,zk,kab,gab,gab,
403     &                   r7,r8,r9,r10,r11,r12)
404#else
405          call nwxc_gvt4_d3(gcab,dgdx,dgdz,xk,zk,kab,gab,gab,
406     &                   r7,r8,r9,r10,r11,r12)
407#endif
408          func(n) = func(n) + gcab*EUEG*wght
409c         dgdPA = dgdx*ChiAP + dgdz*ZAP
410c         dgdGA = dgdx*ChiAG
411c         dgdTA = dgdz*ZAT
412c         dgdPB = dgdx*ChiBP + dgdz*ZBP
413c         dgdGB = dgdx*ChiBG
414c         dgdTB = dgdz*ZBT
415c         EUEGPA = PotLC + P*dLdS*RSP + P*dLdZ*dZdA - EUPA
416c         EUEGPB = PotLC + P*dLdS*RSP + P*dLdZ*dZdB - EUPB
417c         if (ipol.eq.1) then
418c           Amat(n,D1_RA)  = Amat(n,D1_RA)
419c    &                     + (EUEGPA*gcab   + EUEG*dgdPA)*wght
420c           Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + EUEG*dgdGA*wght
421c           Mmat(n,D1_TA)  = Mmat(n,D1_TA)  + EUEG*dgdTA*wght
422c         else
423c           Amat(n,D1_RA)  = Amat(n,D1_RA)
424c    &                     + (EUEGPA*gcab   + EUEG*dgdPA)*wght
425c           Amat(n,D1_RB)  = Amat(n,D1_RB)
426c    &                     + (EUEGPB*gcab   + EUEG*dgdPB)*wght
427c           Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + EUEG*dgdGA*wght
428c           Cmat(n,D1_GBB) = Cmat(n,D1_GBB) + EUEG*dgdGB*wght
429c           Mmat(n,D1_TA)  = Mmat(n,D1_TA)  + EUEG*dgdTA*wght
430c           Mmat(n,D1_TB)  = Mmat(n,D1_TB)  + EUEG*dgdTB*wght
431c         endif
432      endIf
433c      write (*,*) "Amat(n,1),Cmat(n,1),Mmat(n,1)",Amat(n,1),Cmat(n,1)
434c     & ,Mmat(n,1),ipol
435c      stop
43620    continue
437      end
438C>
439C> \brief Compute the same-spin part of the VS98 correlation functional
440C>
441C> This routine evaluates the same-spin part of the VS98 functional for
442C> 1 grid point and 1 spin-case.
443C>
444#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
445#if defined(NWAD_PRINT)
446      Subroutine nwxc_vs98ss_p(tol_rho,PX,GX,TX,F,FP,FG,FT,EUEG,Z,Chi,
447     +                         EUEGP,ChiP,ChiG,ZP,ZT,
448     +                         r13,r14,r15,r16,r17,r18)
449#else
450      Subroutine nwxc_vs98ss(tol_rho,PX,GX,TX,F,FP,FG,FT,EUEG,Z,Chi,
451     +                       EUEGP,ChiP,ChiG,ZP,ZT,
452     +                       r13,r14,r15,r16,r17,r18)
453#endif
454#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV)
455      Subroutine nwxc_vs98ss_d2(tol_rho,PX,GX,TX,F,FP,FG,FT,EUEG,Z,Chi,
456     +                       EUEGP,ChiP,ChiG,ZP,ZT,
457     +                       r13,r14,r15,r16,r17,r18)
458#else
459      Subroutine nwxc_vs98ss_d3(tol_rho,PX,GX,TX,F,FP,FG,FT,EUEG,Z,Chi,
460     +                       EUEGP,ChiP,ChiG,ZP,ZT,
461     +                       r13,r14,r15,r16,r17,r18)
462#endif
463#include "nwad.fh"
464      Implicit none
465c
466#include "intf_nwxc_c_lsda.fh"
467#include "intf_nwxc_gvt4.fh"
468C
469C     Compute the same-spin part of the vs98 correlation functional for one grid
470C     point and one spin-case.
471C
472
473      double precision tol_rho
474#if defined(NWAD_PRINT)
475#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
476      type(nwad_dble)::r13, r14, r15, r16, r17, r18
477#else
478      double precision r13, r14, r15, r16, r17, r18
479#endif
480#else
481      double precision r13, r14, r15, r16, r17, r18
482#endif
483      type(nwad_dble)::PX, GX, TX, F, EUEG, Chi, RS, Z, P, E, PotLC
484      type(nwad_dble)::rhoo, rho53, rho83
485      type(nwad_dble)::kc,xk,zk,d,gc, Zeta
486      double precision FP, FG, FT, DTol, ZP, ZT
487      double precision EUEGP, ChiP, ChiG, cf, gcc
488      double precision Zero, Pt25, F1, F2, F3, F4, F5, F6, F8, F11
489      double precision Pi, Pi34, F13, F23, F43, F53, F83, F113
490      double precision RSP, DX, DZ, dgdP, dgdG, dgdT
491      double precision DP, DG, DT
492      double precision F4o3, dgdx, dgdz
493      double precision d2LdSS, d2LdSZ, d2LdZZ, dLdS, dLdZ
494
495      Data Zero/0.0d0/, Pt25/0.25d0/, F1/1.0d0/, F2/2.0d0/, F3/3.0d0/,
496     $  F4/4.0d0/, F5/5.0d0/, F6/6.0d0/, F8/8.0d0/, F11/11.0d0/,
497     $  gcc/0.00515088d0/,cf/9.115599720d0/
498
499
500      F4o3 = 4.0d0/3.0d0
501      dtol=tol_rho
502      If(PX.le.DTol) then
503        EUEG = Zero
504        Chi = Zero
505        EUEGP = Zero
506        ChiP = Zero
507        ChiG = Zero
508        PX = Zero
509        GX = Zero
510        TX = Zero
511        F  = Zero
512        FP = Zero
513        FG = Zero
514        FT = Zero
515        Z  = Zero
516        ZP = Zero
517        ZT = Zero
518      else
519        Pi = F4*ATan(F1)
520        Pi34 = F3 / (F4*Pi)
521        F13 = F1 / F3
522        F23 = F2 / F3
523        F43 = F2 * F23
524        F53 = F5 / F3
525        F83 = F8 / F3
526        F113 = F11 / F3
527        rhoo = PX
528c       rrho = 1.0d0/rhoo
529c       rho43 = rhoo**F4o3
530c       rho13 = rho43*rrho
531        rho53 = rhoo**F53
532        rho83 = rho53*rhoo
533
534        RS = (Pi34/PX) ** F13
535        Zeta = F1
536#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
537#if defined(NWAD_PRINT)
538        Call nwxc_c_lsda_p(tol_rho,
539     A       RS,Zeta,PotLC,dLdS,dLdZ,d2LdSS,d2LdSZ,d2LdZZ)
540#else
541        Call nwxc_c_lsda(tol_rho,
542     A       RS,Zeta,PotLC,dLdS,dLdZ,d2LdSS,d2LdSZ,d2LdZZ)
543#endif
544#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV)
545        Call nwxc_c_lsda_d2(tol_rho,
546     A       RS,Zeta,PotLC,dLdS,dLdZ,d2LdSS,d2LdSZ,d2LdZZ)
547#else
548        Call nwxc_c_lsda_d3(tol_rho,
549     A       RS,Zeta,PotLC,dLdS,dLdZ,d2LdSS,d2LdSZ,d2LdZZ)
550#endif
551        EUEG = PX*PotLC
552        Chi = GX/rho83
553        Z = (TX/rho53) - cf
554        kc = F1 + gcc*(Chi + Z)
555        xk = Chi/kc
556        zk = Z/kc
557        D = F1 - Chi/(F4*(Z + cf))
558#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
559#if defined(NWAD_PRINT)
560        call nwxc_gvt4_p(gc,dgdx,dgdz,xk,zk,kc,gcc,gcc,
561     &                 r13,r14,r15,r16,r17,r18)
562#else
563        call nwxc_gvt4(gc,dgdx,dgdz,xk,zk,kc,gcc,gcc,
564     &                 r13,r14,r15,r16,r17,r18)
565#endif
566#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV)
567        call nwxc_gvt4_d2(gc,dgdx,dgdz,xk,zk,kc,gcc,gcc,
568     &                 r13,r14,r15,r16,r17,r18)
569#else
570        call nwxc_gvt4_d3(gc,dgdx,dgdz,xk,zk,kc,gcc,gcc,
571     &                 r13,r14,r15,r16,r17,r18)
572#endif
573        E = D*EUEG*gc
574c         write (*,*) "Chi, Z, gc", CHi, Z, gc
575        F = E
576c
577c       RSP = -RS/(F3*Px)
578c       ChiG = F1/PX**F83
579c       ChiP = -F83*Chi/PX
580c       ZP = -F53 * TX/rho83
581c       ZT =  F1/rho53
582c       DZ = Chi/(F4*(Z + cf)*(Z + cf))
583c       DX = -F1/(F4*(Z + cf))
584c       DP = DZ*ZP + DX*ChiP
585c       DG = DX*ChiG
586c       DT = DZ*ZT
587c       dgdP = dgdx*ChiP + dgdz*ZP
588c       dgdG = dgdx*ChiG
589c       dgdT = dgdz*ZT
590c       EUEGP = PotLC + PX*dLdS*RSP
591c       FP = DP*EUEG*gc + D*EUEGP*gc + D*EUEG*dgdP
592c       FG = DG*EUEG*gc + D*EUEG*dgdG
593c       FT = DT*EUEG*gc + D*EUEG*dgdT
594       Endif
595       Return
596       End
597#ifndef NWAD_PRINT
598#define NWAD_PRINT
599c
600c     Compile source again for Maxima
601c
602#include "nwxc_c_vs98.F"
603#endif
604#ifndef SECOND_DERIV
605#define SECOND_DERIV
606c
607c     Compile source again for the 2nd derivative case
608c
609#include "nwxc_c_vs98.F"
610#endif
611#ifndef THIRD_DERIV
612#define THIRD_DERIV
613c
614c     Compile source again for the 3rd derivative case
615c
616#include "nwxc_c_vs98.F"
617#endif
618#undef NWAD_PRINT
619C> @}
620