1*     ************************************************
2*     *                                              *
3*     *                   nwpw_GVT4                  *
4*     *                                              *
5*     ************************************************
6      subroutine nwpw_GVT4(ag,bg,cg,dg,eg,fg,alphag,alphagt,
7     >                     xg,zg,gammag,G,dGdx,dGdz)
8      implicit none
9      real*8 ag,bg,cg,dg,eg,fg,alphag,alphagt,xg,zg,gammag
10      real*8 G,dGdx,dGdz
11      real*8 xg2,zg2,gammag2
12
13      xg2     = xg*xg
14      zg2     = zg*zg
15      gammag2 = gammag*gammag
16
17      G = (ag + bg*xg + cg*zg
18     &     + dg*xg2 + eg*xg*zg + fg*zg2)/gammag
19
20      dGdx = ( -ag*alphag
21     &        + bg*(1.0d0 - 2.0d0*alphag*xg)
22     &        - 2.0d0*cg*alphag*zg
23     &        + dg*(2.0d0*xg - 3.0d0*alphag*xg2)
24     &        + eg*(zg - 3.0d0*alphag*zg*xg)
25     &        - 3.0d0*fg*alphag*zg2)/gammag2
26
27      dGdz = ( -ag*alphagt
28     &        - 2.0d0*bg*alphagt*xg
29     &        + cg*(1.0d0 - 2.0d0*alphagt*zg)
30     &        - 3.0d0*dg*alphagt*xg2
31     &        + eg*(xg - 3.0d0*alphagt*xg*zg)
32     &        + fg*(2.0d0*zg - 3.0d0*alphagt*zg2))/gammag2
33
34      return
35      end
36*     ************************************************
37*     *                                              *
38*     *              gen_VS98_restricted             *
39*     *                                              *
40*     ************************************************
41
42*    This function returns the VS98 exchange-correlation
43*  energy density, xce, and its derivatives with respect
44*  to n, |grad n|, tau.
45
46*
47*   Entry - n2ft3d   : number of grid points
48*           rho_in(*) :  density (nup+ndn)
49*           agr_in(*): |grad rho_in|
50*           tau_in(*): tau
51*           x_parameter: scale parameter for exchange
52*           c_parameter: scale parameter for correlation
53*
54*     Exit  - xce(n2ft3d) : VS98 exchange correlation energy density
55*             fn(n2ft3d)  : d(n*xce)/dn
56*             fdn(n2ft3d) : d(n*xce)/d|grad n|
57*             fdtau(n2ft3d) : d(n*xce)/dtau
58*
59      subroutine gen_VS98_restricted(n2ft3d,rho_in,agr_in,tau_in,
60     >                               x_parameter,c_parameter,
61     >                               xce,fn,fdn,fdtau)
62      implicit none
63*     ***** input *****
64      integer n2ft3d
65      real*8 rho_in(*),agr_in(*),tau_in(*)
66      real*8 x_parameter,c_parameter
67*     ***** output *****
68      real*8 xce(*),fn(*),fdn(*),fdtau(*)
69*     ***** local declarations *****
70      integer i
71      real*8 n,agr,tau
72      real*8 xr,zr,gamma
73      real*8 inv_n,n_13,n_43,n_53,n_83,agr2
74      real*8 x,dx_dn,dx_dagr,z,dz_dn,dz_dtau
75      real*8 GG,dGdx,dGdz,dG_dn,dG_dagr,dG_dtau
76      real*8 ex,fnx,fdnx,fdtaux
77      real*8 ess0c,dess0c_drs,dess0c_dn
78      real*8 eud0c,deud0c_drs
79      real*8 eudc,eud1c,fud1c,dfud1c_dn
80      real*8 dfudc_dn,dfudc_dagr,dfudc_dtau
81      real*8 rs,drs_dn,dummy
82      real*8 D,dD_dx,dD_dz,dD_dn,dD_dagr,dD_dtau
83      real*8 Dt,dD_dDt,dDt_dx,dDt_dz
84      real*8 ec,fnc,fdnc,fdtauc
85*     ***** constants *****
86      real*8 pi,thrd,twthrd,frthrd,fvthrd,etthrd
87      parameter (pi     = 3.14159265358979311599d0)
88      parameter (thrd   = 1.0d0/3.0d0)
89      parameter (twthrd = 2.0d0/3.0d0)
90      parameter (frthrd = 4.0d0/3.0d0)
91      parameter (fvthrd = 5.0d0/3.0d0)
92      parameter (etthrd = 8.0d0/3.0d0)
93*     ***** density cutoff parametersi *****
94      real*8 tol,ETA,thresd
95      parameter (tol = 1.0d-10)
96      parameter (ETA            =      1.0d-20)
97      parameter (thresd  =   1.0d-2)
98*     ***** VS98 constants *****
99      real*8 aax,bbx,ccx,ddx,eex,ffx,alphax
100      real*8 aass,bbss,ccss,ddss,eess,ffss,alphass
101      real*8 aaopp,bbopp,ccopp,ddopp,eeopp,ffopp,alphaopp
102      real*8 cf,clda
103      parameter (cf       =  9.115599720d0)
104c     cf = 3.0d0/5.0d0*(6.0d0*pi*pi)**(2.0d0/3.0d0))
105      parameter (clda     =  0.9305257363491d0)
106      parameter (aax      = -9.800683d-1)
107      parameter (bbx      = -3.556788d-3)
108      parameter (ccx      =  6.250326d-3)
109      parameter (ddx      = -2.354518d-5)
110      parameter (eex      = -1.282732d-4)
111      parameter (ffx      =  3.574822d-4)
112      parameter (alphax   =  0.00186726d0)
113      parameter (aass     =  3.270912d-1)
114      parameter (bbss     = -3.228915d-2)
115      parameter (ccss     = -2.942406d-2)
116      parameter (ddss     =  2.134222d-3)
117      parameter (eess     = -5.451559d-3)
118      parameter (ffss     =  1.577575d-2)
119      parameter (alphass  =  0.00515088d0)
120      parameter (aaopp    =  7.035010d-1)
121      parameter (bbopp    =  7.694574d-3)
122      parameter (ccopp    =  5.152765d-2)
123      parameter (ddopp    =  3.394308d-5)
124      parameter (eeopp    = -1.269420d-3)
125      parameter (ffopp    =  1.296118d-3)
126      parameter (alphaopp =  0.00304966d0)
127
128      do i=1,n2ft3d
129        n        = rho_in(i) + ETA
130        agr      = agr_in(i) + ETA
131        tau      = 2.0d0*tau_in(i) + ETA
132
133        n   = 0.50d0*n
134        agr = 0.50d0*agr
135        tau = 0.50d0*tau
136
137        agr2  = agr*agr
138        inv_n = 1.0d0/n
139        n_13  = n**thrd
140        n_43  = n_13*n
141        n_53  = n_43*n_13
142        n_83  = n_53*n
143
144        x       =  agr2/n_83
145        dx_dn   = -etthrd*x*inv_n
146        dx_dagr =  2.0d0*agr/n_83
147        z       =  tau/n_53 - cf
148        dz_dn   = -fvthrd*tau/n_83
149        dz_dtau =  1.0d0/n_53
150
151*       ***** VS98 Exchange *****
152        gamma = 1.0d0 + alphax*(x + z)
153        xr    = x/gamma
154        zr    = z/gamma
155
156        call nwpw_GVT4(aax,bbx,ccx,ddx,eex,ffx,alphax,alphax,
157     >              xr,zr,gamma,GG,dGdx,dGdz)
158
159        dG_dn   = dGdx*dx_dn + dGdz*dz_dn
160        dG_dagr = dGdx*dx_dagr
161        dG_dtau = dGdz*dz_dtau
162
163        ex     = clda*n_13*GG
164        fnx    = frthrd*ex + clda*n_43*dG_dn
165        fdnx   = clda*n_43*dG_dagr
166        fdtaux = clda*n_43*dG_dtau
167
168*       ***** VS98 Correlation *****
169*       ***** Same-Spin (alpha-alpha/beta-beta) *****
170        rs     = (0.75d0/(pi*n))**thrd
171        drs_dn = -thrd*rs/n
172
173        call gen_PW91_c_rz(tol,rs,1.0d0,ess0c,dess0c_drs,dummy)
174        dess0c_dn = dess0c_drs*drs_dn
175
176        Dt      = 1.0d0 - 0.25d0*x/(z + cf)
177        D       = 0.0d0
178        dD_dn   = 0.0d0
179        dD_dagr = 0.0d0
180        dD_dtau = 0.0d0
181
182        if (Dt .le. 0.0d0) then
183           D       = 0.0d0
184           dD_dn   = 0.0d0
185           dD_dagr = 0.0d0
186           dD_dtau = 0.0d0
187        else if (Dt .gt. 0.0d0 .and. Dt .lt. thresd) then
188           D       =  2.0d0*Dt*Dt/thresd - Dt*Dt*Dt/(thresd*thresd)
189           dD_dDt  =  4.0d0*Dt/thresd - 3.0d0*Dt*Dt/(thresd*thresd)
190           dDt_dx  = -0.25d0/(z + cf)
191           dDt_dz  =  0.25d0*x/((z + cf)*(z + cf))
192           dD_dn   =  dD_dDt*(dDt_dx*dx_dn + dDt_dz*dz_dn)
193           dD_dagr =  dD_dDt*dDt_dx*dx_dagr
194           dD_dtau =  dD_dDt*dDt_dz*dz_dtau
195        else if(Dt .ge. thresd) then
196           D       =  Dt
197           dD_dx   = -0.25d0/(z + cf)
198           dD_dz   =  0.25d0*x/((z + cf)*(z + cf))
199           dD_dn   =  dD_dx*dx_dn + dD_dz*dz_dn
200           dD_dagr =  dD_dx*dx_dagr
201           dD_dtau =  dD_dz*dz_dtau
202        end if
203
204        gamma = 1.0d0 + alphass*(x + z)
205        xr = x/gamma
206        zr = z/gamma
207
208        call nwpw_GVT4(aass,bbss,ccss,ddss,eess,ffss,alphass,alphass,
209     >                 xr,zr,gamma,GG,dGdx,dGdz)
210
211        dG_dn   = dGdx*dx_dn + dGdz*dz_dn
212        dG_dagr = dGdx*dx_dagr
213        dG_dtau = dGdz*dz_dtau
214
215        ec     = ess0c*GG*D
216        fnc    = ec + n*dess0c_drs*drs_dn*GG*D
217     &         + n*ess0c*dG_dn*D
218     &         + n*ess0c*GG*dD_dn
219        fdnc   = n*ess0c*(dG_dagr*D + GG*dD_dagr)
220        fdtauc = n*ess0c*(dG_dtau*D + GG*dD_dtau)
221
222*       ***** Opposite-Spin (alpha-beta) *****
223        n = 2.0d0*n
224
225        rs         =  (0.75d0/(pi*n))**thrd
226        drs_dn     = -thrd*rs/n
227
228        call gen_PW91_c_rz(tol,rs,0.0d0,eud0c,deud0c_drs,dummy)
229        eud1c      = eud0c - ess0c
230        fud1c      = n*eud1c
231        dfud1c_dn  = eud0c + n*deud0c_drs*drs_dn
232     &             - ess0c - 0.50d0*n*dess0c_dn
233
234        x   = 2.0d0*x
235        z   = 2.0d0*z
236
237        gamma = 1.0d0 + alphaopp*(x + z)
238        xr    = x/gamma
239        zr    = z/gamma
240
241        call nwpw_GVT4(aaopp,bbopp,ccopp,ddopp,eeopp,ffopp,
242     >                 alphaopp,alphaopp,
243     >                 xr,zr,gamma,GG,dGdx,dGdz)
244
245        dG_dn   = dGdx*dx_dn + dGdz*dz_dn
246        dG_dagr = dGdx*dx_dagr
247        dG_dtau = dGdz*dz_dtau
248
249        eudc       = eud1c*GG
250        dfudc_dn   = fud1c*dG_dn    + GG*dfud1c_dn
251        dfudc_dagr = fud1c*dG_dagr
252        dfudc_dtau = fud1c*dG_dtau
253
254        ec     = ec     + eudc
255        fnc    = fnc    + dfudc_dn
256        fdnc   = fdnc   + dfudc_dagr
257        fdtauc = fdtauc + dfudc_dtau
258
259        xce(i)   = x_parameter*ex     + c_parameter*ec
260        fn(i)    = x_parameter*fnx    + c_parameter*fnc
261        fdn(i)   = x_parameter*fdnx   + c_parameter*fdnc
262        fdtau(i) = x_parameter*fdtaux + c_parameter*fdtauc
263      end do
264
265      return
266      end
267*     ************************************************
268*     *                                              *
269*     *              gen_VS98_unrestricted             *
270*     *                                              *
271*     ************************************************
272
273*    This function returns the VS98 exchange-correlation
274*  energy density, xce, and its derivatives with respect
275*  to nup, ndn, |grad nup|, |grad ndn|, tauup, taudn.
276
277*
278*   Entry - n2ft3d   : number of grid points
279*           rho_in(*,2) :  density (nup and ndn)
280*           agr_in(*,3): |grad rho_in| (nup, ndn and n)
281*           tau_in(*,2): tau (nup and ndn)
282*           x_parameter: scale parameter for exchange
283*           c_parameter: scale parameter for correlation
284*
285*     Exit  - xce(n2ft3d) : VS98 exchange correlation energy density
286*             fn(n2ft3d,2)  : d(n*xce)/dnup, d(n*xce)/dndn
287*             fdn(n2ft3d,3) : d(n*xce)/d|grad nup|, d(n*xce)/d|grad ndn|, d(n*xce)/d|grad n|
288*             fdtau(n2ft3d,2) : d(n*xce)/dtauup, d(n*xce)/dtaudn
289*
290      subroutine gen_VS98_unrestricted(n2ft3d,rho_in,agr_in,tau_in,
291     >                               x_parameter,c_parameter,
292     >                               xce,fn,fdn,fdtau)
293      implicit none
294*     ***** input *****
295      integer n2ft3d
296      real*8 rho_in(n2ft3d,2),agr_in(n2ft3d,3),tau_in(n2ft3d,2)
297      real*8 x_parameter,c_parameter
298*     ***** output *****
299      real*8 xce(n2ft3d),fn(n2ft3d,2),fdn(n2ft3d,3),fdtau(n2ft3d,2)
300*     ***** local declarations *****
301      integer i
302      real*8 nup,agrup,tauup
303      real*8 ndn,agrdn,taudn
304      real*8 n,agr,tau
305      real*8 x,dxdn,dxdagr,z,dzdn,dzdtau
306      real*8 xr,zr,gamma
307      real*8 inv_nup,nup_13,nup_43,nup_53,nup_83,agrup2
308      real*8 inv_ndn,ndn_13,ndn_43,ndn_53,ndn_83,agrdn2
309      real*8 xup,dxup_dnup,dxup_dagrup,zup,dzup_dnup,dzup_dtauup
310      real*8 xdn,dxdn_dndn,dxdn_dagrdn,zdn,dzdn_dndn,dzdn_dtaudn
311      real*8 GG,dGdx,dGdz
312      real*8 eupx,fnupx,fdnupx,fdtauupx
313      real*8 ednx,fndnx,fdndnx,fdtaudnx
314      real*8 euu0c,deuu0c_dnup
315      real*8 euuc,dfuuc_dnup,dfuuc_dagrup,dfuuc_dtauup
316      real*8 edd0c,dedd0c_dndn
317      real*8 eddc,dfddc_dndn,dfddc_dagrdn,dfddc_dtaudn
318      real*8 eud0c,deud0c_drs,deud0c_dzeta
319      real*8 eudc,eud1c,fud1c,dfud1c_dnup,dfud1c_dndn
320      real*8 dfudc_dnup,dfudc_dagrup,dfudc_dtauup
321      real*8 dfudc_dndn,dfudc_dagrdn,dfudc_dtaudn
322      real*8 rs,drs_dn,zeta,dzeta_dnup,dzeta_dndn,deuu0c_drs,dedd0c_drs
323      real*8 D,dD_dx,dD_dz,dD_dn,dD_dagr,dD_dtau,dummy
324      real*8 Dt,dD_dDt,dDt_dx,dDt_dz
325      real*8 eupc,fnupc,fdnupc,fdtauupc
326      real*8 ednc,fndnc,fdndnc,fdtaudnc
327      real*8 ex,ec
328*     ***** constants *****
329      real*8 pi,thrd,twthrd,frthrd,fvthrd,etthrd
330      parameter (pi     = 3.14159265358979311599d0)
331      parameter (thrd   = 1.0d0/3.0d0)
332      parameter (twthrd = 2.0d0/3.0d0)
333      parameter (frthrd = 4.0d0/3.0d0)
334      parameter (fvthrd = 5.0d0/3.0d0)
335      parameter (etthrd = 8.0d0/3.0d0)
336*     ***** density cutoff parameters *****
337      real*8 tol,ETA,thresd
338      parameter (tol = 1.0d-10)
339      parameter (ETA            =      1.0d-20)
340      parameter (thresd  =   1.0d-2)
341*     ***** VS98 constants *****
342      real*8 aax,bbx,ccx,ddx,eex,ffx,alphax
343      real*8 aass,bbss,ccss,ddss,eess,ffss,alphass
344      real*8 aaopp,bbopp,ccopp,ddopp,eeopp,ffopp,alphaopp
345      real*8 cf,clda
346      parameter (cf       =  9.115599720d0)
347c     cf = 3.0d0/5.0d0*(6.0d0*pi*pi)**(2.0d0/3.0d0))
348      parameter (clda     =  0.9305257363491d0)
349      parameter (aax      = -9.800683d-1)
350      parameter (bbx      = -3.556788d-3)
351      parameter (ccx      =  6.250326d-3)
352      parameter (ddx      = -2.354518d-5)
353      parameter (eex      = -1.282732d-4)
354      parameter (ffx      =  3.574822d-4)
355      parameter (alphax   =  0.00186726d0)
356      parameter (aass     =  3.270912d-1)
357      parameter (bbss     = -3.228915d-2)
358      parameter (ccss     = -2.942406d-2)
359      parameter (ddss     =  2.134222d-3)
360      parameter (eess     = -5.451559d-3)
361      parameter (ffss     =  1.577575d-2)
362      parameter (alphass  =  0.00515088d0)
363      parameter (aaopp    =  7.035010d-1)
364      parameter (bbopp    =  7.694574d-3)
365      parameter (ccopp    =  5.152765d-2)
366      parameter (ddopp    =  3.394308d-5)
367      parameter (eeopp    = -1.269420d-3)
368      parameter (ffopp    =  1.296118d-3)
369      parameter (alphaopp =  0.00304966d0)
370
371      do i=1,n2ft3d
372        nup        = rho_in(i,1) + ETA
373        agrup      = agr_in(i,1) + ETA
374        tauup      = tau_in(i,1) + ETA
375        ndn        = rho_in(i,2) + ETA
376        agrdn      = agr_in(i,2) + ETA
377        taudn      = tau_in(i,2) + ETA
378
379        n = nup + ndn
380
381*       ***** VS98 Exchange *****
382*       ***** Up *****
383        agrup2  = agrup*agrup
384        inv_nup = 1.0d0/nup
385        nup_13  = nup**thrd
386        nup_43  = nup_13*nup
387        nup_53  = nup_43*nup_13
388        nup_83  = nup_53*nup
389
390        xup         =  agrup2/nup_83
391        dxup_dnup   = -etthrd*xup*inv_nup
392        dxup_dagrup =  2.0d0*agrup/nup_83
393        zup         =  tauup/nup_53 - cf
394        dzup_dnup   = -fvthrd*tauup/nup_83
395        dzup_dtauup =  1.0d0/nup_53
396
397        gamma = 1.0d0 + alphax*(xup + zup)
398        xr    = xup/gamma
399        zr    = zup/gamma
400
401        call nwpw_GVT4(aax,bbx,ccx,ddx,eex,ffx,alphax,alphax,
402     >              xr,zr,gamma,GG,dGdx,dGdz)
403
404        eupx     = clda*nup_13*GG
405        fnupx    = frthrd*eupx + clda*nup_43*(dGdx*dxup_dnup
406     &           + dGdz*dzup_dnup)
407        fdnupx   = clda*nup_43*(dGdx*dxup_dagrup)
408        fdtauupx = clda*nup_43*(dGdz*dzup_dtauup)
409
410*       ***** Down *****
411        agrdn2  = agrdn*agrdn
412        inv_ndn = 1.0d0/ndn
413        ndn_13  = ndn**thrd
414        ndn_43  = ndn_13*ndn
415        ndn_53  = ndn_43*ndn_13
416        ndn_83  = ndn_53*ndn
417
418        xdn         =  agrdn2/ndn_83
419        dxdn_dndn   = -etthrd*xdn*inv_ndn
420        dxdn_dagrdn =  2.0d0*agrdn/ndn_83
421        zdn         =  taudn/ndn_53 - cf
422        dzdn_dndn   = -fvthrd*taudn/ndn_83
423        dzdn_dtaudn =  1.0d0/ndn_53
424
425        gamma = 1.0d0 + alphax*(xdn + zdn)
426        xr    = xdn/gamma
427        zr    = zdn/gamma
428
429        call nwpw_GVT4(aax,bbx,ccx,ddx,eex,ffx,alphax,alphax,
430     >              xr,zr,gamma,GG,dGdx,dGdz)
431
432        ednx     = clda*ndn_13*GG
433        fndnx    = frthrd*ednx + clda*ndn_43*(dGdx*dxdn_dndn
434     &           + dGdz*dzdn_dndn)
435        fdndnx   = clda*ndn_43*(dGdx*dxdn_dagrdn)
436        fdtaudnx = clda*ndn_43*(dGdz*dzdn_dtaudn)
437
438        ex = (eupx*nup + ednx*ndn)/n
439
440*       ***** VS98 Correlation *****
441*       ***** Same-Spin *****
442*       ***** alpha-alpha *****
443        rs     = (0.75d0/(pi*nup))**thrd
444        drs_dn = -thrd*rs/nup
445
446        call gen_PW91_c_rz(tol,rs,1.0d0,euu0c,deuu0c_drs,dummy)
447        deuu0c_dnup = deuu0c_drs*drs_dn
448
449        Dt      = 1.0d0 - 0.25d0*xup/(zup + cf)
450        D       = 0.0d0
451        dD_dn   = 0.0d0
452        dD_dagr = 0.0d0
453        dD_dtau = 0.0d0
454
455        if (Dt .le. 0.0d0) then
456           D       = 0.0d0
457           dD_dn   = 0.0d0
458           dD_dagr = 0.0d0
459           dD_dtau = 0.0d0
460        else if (Dt .gt. 0.0d0 .and. Dt .lt. thresd) then
461           D       =  2.0d0*Dt*Dt/thresd - Dt*Dt*Dt/(thresd*thresd)
462           dD_dDt  =  4.0d0*Dt/thresd - 3.0d0*Dt*Dt/(thresd*thresd)
463           dDt_dx  = -0.25d0/(zup + cf)
464           dDt_dz  =  0.25d0*xup/((zup + cf)*(zup + cf))
465           dD_dn   =  dD_dDt*(dDt_dx*dxup_dnup + dDt_dz*dzup_dnup)
466           dD_dagr =  dD_dDt*dDt_dx*dxup_dagrup
467           dD_dtau =  dD_dDt*dDt_dz*dzup_dtauup
468        else if(Dt .ge. thresd) then
469           D       =  Dt
470           dD_dx   = -0.25d0/(zup + cf)
471           dD_dz   =  0.25d0*xup/((zup + cf)*(zup + cf))
472           dD_dn   =  dD_dx*dxup_dnup + dD_dz*dzup_dnup
473           dD_dagr =  dD_dx*dxup_dagrup
474           dD_dtau =  dD_dz*dzup_dtauup
475        end if
476
477        gamma = 1.0d0 + alphass*(xup + zup)
478        xr    = xup/gamma
479        zr    = zup/gamma
480
481        call nwpw_GVT4(aass,bbss,ccss,ddss,eess,ffss,alphass,alphass,
482     >                 xr,zr,gamma,GG,dGdx,dGdz)
483
484        eupc     = euu0c*GG*D
485        fnupc    = eupc + nup*deuu0c_drs*drs_dn*GG*D
486     &           + nup*euu0c*(dGdx*dxup_dnup+dGdz*dzup_dnup)*D
487     &           + nup*euu0c*GG*dD_dn
488        fdnupc   = nup*euu0c*(dGdx*dxup_dagrup*D + GG*dD_dagr)
489        fdtauupc = nup*euu0c*(dGdz*dzup_dtauup*D + GG*dD_dtau)
490
491*       ***** beta-beta *****
492        rs     = (0.75d0/(pi*ndn))**thrd
493        drs_dn = -thrd*rs/ndn
494
495        call gen_PW91_c_rz(tol,rs,1.0d0,edd0c,dedd0c_drs,dummy)
496        dedd0c_dndn = dedd0c_drs*drs_dn
497
498        Dt      = 1.0d0 - 0.25d0*xdn/(zdn + cf)
499        D       = 0.0d0
500        dD_dn   = 0.0d0
501        dD_dagr = 0.0d0
502        dD_dtau = 0.0d0
503
504        if (Dt .le. 0.0d0) then
505           D       = 0.0d0
506           dD_dn   = 0.0d0
507           dD_dagr = 0.0d0
508           dD_dtau = 0.0d0
509        else if (Dt .gt. 0.0d0 .and. Dt .lt. thresd) then
510           D       =  2.0d0*Dt*Dt/thresd - Dt*Dt*Dt/(thresd*thresd)
511           dD_dDt  =  4.0d0*Dt/thresd - 3.0d0*Dt*Dt/(thresd*thresd)
512           dDt_dx  = -0.25d0/(zdn + cf)
513           dDt_dz  =  0.25d0*xdn/((zdn + cf)*(zdn + cf))
514           dD_dn   =  dD_dDt*(dDt_dx*dxdn_dndn + dDt_dz*dzdn_dndn)
515           dD_dagr =  dD_dDt*dDt_dx*dxdn_dagrdn
516           dD_dtau =  dD_dDt*dDt_dz*dzdn_dtaudn
517        else if(Dt .ge. thresd) then
518           D       =  Dt
519           dD_dx   = -0.25d0/(zdn + cf)
520           dD_dz   =  0.25d0*xdn/((zdn + cf)*(zdn + cf))
521           dD_dn   =  dD_dx*dxdn_dndn + dD_dz*dzdn_dndn
522           dD_dagr =  dD_dx*dxdn_dagrdn
523           dD_dtau =  dD_dz*dzdn_dtaudn
524        end if
525
526        gamma = 1.0d0 + alphass*(xdn + zdn)
527        xr    = xdn/gamma
528        zr    = zdn/gamma
529
530        call nwpw_GVT4(aass,bbss,ccss,ddss,eess,ffss,alphass,alphass,
531     >                 xr,zr,gamma,GG,dGdx,dGdz)
532
533        ednc     = edd0c*GG*D
534        fndnc    = ednc + ndn*dedd0c_drs*drs_dn*GG*D
535     &           + ndn*edd0c*(dGdx*dxdn_dndn+dGdz*dzdn_dndn)*D
536     &           + ndn*edd0c*GG*dD_dn
537        fdndnc   = ndn*edd0c*(dGdx*dxdn_dagrdn*D + GG*dD_dagr)
538        fdtaudnc = ndn*edd0c*(dGdz*dzdn_dtaudn*D + GG*dD_dtau)
539
540        ec = (eupc*nup + ednc*ndn)/n
541
542*       ***** Opposite-Spin (alpha-beta) *****
543        rs         =  (0.75d0/(pi*n))**thrd
544        drs_dn     = -thrd*rs/n
545        zeta       =  (nup - ndn)/n
546        dzeta_dnup =  ( 1.0d0 - zeta)/n
547        dzeta_dndn =  (-1.0d0 - zeta)/n
548
549        call gen_PW91_c_rz(tol,rs,zeta,eud0c,deud0c_drs,deud0c_dzeta)
550        eud1c       = eud0c - (euu0c*nup + edd0c*ndn)/n
551        fud1c       = n*eud1c
552        dfud1c_dnup = eud0c + n*(deud0c_drs*drs_dn
553     &              + deud0c_dzeta*dzeta_dnup) - euu0c
554     &              - nup*deuu0c_dnup
555        dfud1c_dndn = eud0c + n*(deud0c_drs*drs_dn
556     &              + deud0c_dzeta*dzeta_dndn) - edd0c
557     &              - ndn*dedd0c_dndn
558
559        x   = xup + xdn
560        z   = zup + zdn
561
562        gamma = 1.0d0 + alphaopp*(x + z)
563        xr    = x/gamma
564        zr    = z/gamma
565
566        call nwpw_GVT4(aaopp,bbopp,ccopp,ddopp,eeopp,ffopp,
567     >                 alphaopp,alphaopp,
568     >                 xr,zr,gamma,GG,dGdx,dGdz)
569
570        eudc         = eud1c*GG
571        dfudc_dnup   = dfud1c_dnup*GG
572     &               + fud1c*(dGdx*dxup_dnup + dGdz*dzup_dnup)
573        dfudc_dndn   = dfud1c_dndn*GG +
574     &               fud1c*(dGdx*dxdn_dndn + dGdz*dzdn_dndn)
575        dfudc_dagrup = fud1c*dGdx*dxup_dagrup
576        dfudc_dagrdn = fud1c*dGdx*dxdn_dagrdn
577        dfudc_dtauup = fud1c*dGdz*dzup_dtauup
578        dfudc_dtaudn = fud1c*dGdz*dzdn_dtaudn
579
580        ec       = ec       + eudc
581        fnupc    = fnupc    + dfudc_dnup
582        fndnc    = fndnc    + dfudc_dndn
583        fdnupc   = fdnupc   + dfudc_dagrup
584        fdndnc   = fdndnc   + dfudc_dagrdn
585        fdtauupc = fdtauupc + dfudc_dtauup
586        fdtaudnc = fdtaudnc + dfudc_dtaudn
587
588        xce(i)     = x_parameter*ex       + c_parameter*ec
589        fn(i,1)    = x_parameter*fnupx    + c_parameter*fnupc
590        fn(i,2)    = x_parameter*fndnx    + c_parameter*fndnc
591        fdn(i,1)   = x_parameter*fdnupx   + c_parameter*fdnupc
592        fdn(i,2)   = x_parameter*fdndnx   + c_parameter*fdndnc
593        fdn(i,3)   = 0.0d0
594        fdtau(i,1) = x_parameter*fdtauupx + c_parameter*fdtauupc
595        fdtau(i,2) = x_parameter*fdtaudnx + c_parameter*fdtaudnc
596      end do
597
598      return
599      end
600
601