1*     ************************************************
2*     *                                              *
3*     *              nwpw_tpss03_x                   *
4*     *                                              *
5*     ************************************************
6      subroutine nwpw_tpss03_x(pi,thrd,twthrd,frthrd,fvthrd,etthrd,
7     >                         C920,C1,C2,C3,kappa,mu,b,c,e,
8     >                         Cx,P23,es,
9     >                         n,agr,tau,
10     >                         xe,dfdnx,dfdagrx,dfdtaux)
11      implicit none
12*     ***** input *****
13      real*8 pi,thrd,twthrd,frthrd,fvthrd,etthrd,C920,C1,C2,C3
14      real*8 kappa,mu,b,c,e
15      real*8 Cx,P23,es
16      real*8 n,agr,tau
17*     ***** output *****
18      real*8 xe,dfdnx,dfdagrx,dfdtaux
19*     ***** local declarations *****
20      real*8 n_13,n_53,n_83,inv_n,agr2,tauW,tauU
21      real*8 p,z,fz,alpha,qb,p2,p3,z2,z3,qb2,fa,thresA,dalpha_dfa
22      real*8 x1a,x1,x2,x3a,x3b,x3,x4,x5,x6,x7,xnum,x
23      real*8 dp_dn,dz_dn,dalpha_dn,qba,qbb,dqba_dn,dqbb_dn,dqb_dn
24      real*8 dx1_dn,dx2_dn,dx3_dn,dx4_dn,dx5_dn,dx6_dn,dx7_dn
25      real*8 dx3a_dn,dx3b_dn,dxnum_dn,dx_dn
26      real*8 dp_dagr,dz_dagr,dalpha_dagr,dqba_dagr,dqbb_dagr,dqb_dagr
27      real*8 dx1_dagr,dx2_dagr,dx3_dagr,dx4_dagr,dx5_dagr,dx6_dagr
28      real*8 dx7_dagr,dx3a_dagr,dx3b_dagr,dxnum_dagr,dx_dagr
29      real*8 dz_dtau,dalpha_dtau,dqb_dtau
30      real*8 dx1_dtau,dx2_dtau,dx3_dtau,dx5_dtau,dxnum_dtau,dx_dtau
31      real*8 Fx,ex0,nex0,dFx_dx,dFx_dn,dFx_dagr,dFx_dtau
32
33
34      n_13  = n**thrd
35      n_53  = n_13*n_13*n
36      n_83  = n_53*n
37      inv_n = 1.0d0/n
38      agr2  = agr*agr
39
40      p       =  agr2/(4.0d0*P23*n_83)
41      dp_dn   = -etthrd*p*inv_n
42      dp_dagr =  2.0d0*p/agr
43c     dp_dtau =  0.0d0
44
45      p2 = p*p
46      p3 = p2*p
47
48      tauW  = 0.125d0*agr2*inv_n
49      tauU  = 0.3d0*P23*n_53
50      fz    = tauW/tau
51
52      if (fz .gt. 1.0d0) then
53         z       = 1.0d0
54         dz_dn   = 0.0d0
55         dz_dagr = 0.0d0
56         dz_dtau = 0.0d0
57      else
58         z       =  fz
59         dz_dn   = -z*inv_n
60         dz_dagr =  2.0d0*z/agr
61         dz_dtau = -z/tau
62      end if
63
64      z2 = z*z
65      z3 = z2*z
66
67      alpha = fvthrd*p*(1.0d0/z - 1.0d0)
68c     alpha = (tau - tauW)/tauU
69
70      if (alpha .le. 0.0d0) then
71        alpha       = 0.0d0
72        dalpha_dn   = 0.0d0
73        dalpha_dagr = 0.0d0
74        dalpha_dtau = 0.0d0
75      else
76        dalpha_dn   = fvthrd*(-p*dz_dn/z2 + dp_dn*(1.0d0/z - 1.0d0))
77        dalpha_dagr = (alpha/p)*dp_dagr - fvthrd*(p/z2)*dz_dagr
78        dalpha_dtau = 1.0d0/tauU
79c       dalpha_dtau = fvthrd*p*(-1.0d0/z2)*dz_dtau
80      end if
81
82      qb = C920*(alpha - 1.0d0)/dsqrt(1.0d0+b*alpha*(alpha - 1.0d0))
83     &   + twthrd*p
84
85      qb2 = qb*qb
86
87      x1a = (c*z2)/((1.0d0 + z2)**2.0d0)
88      x1  = (C1 + x1a)*p
89      x2  = C2*qb2
90      x3a = C3*qb
91      x3b = dsqrt(0.5d0*(0.36d0*z2 + p2))
92      x3  = x3a*x3b
93      x4  = C1*C1*p2/kappa
94      x5  = 2.0d0*es*C1*0.36d0*z2
95      x6  = e*mu*p3
96      x7  = 1.0d0/((1.0d0 + es*p)**2.0d0)
97
98      xnum = x1 + x2 + x3 + x4 + x5 + x6
99      x    = xnum*x7
100
101*     ***** dxdn *****
102      qba     = alpha - 1.0d0
103      qbb     = 1.0d0/dsqrt(1.0d0 + b*alpha*(alpha - 1.0d0))
104      dqba_dn = dalpha_dn
105      dqbb_dn = b*dalpha_dn*(2.0d0*alpha - 1.0d0)*
106     &          (-0.5d0)*(qbb**3.0d0)
107      dqb_dn  = C920*(qbb*dqba_dn + qba*dqbb_dn) + twthrd*dp_dn
108
109      dx1_dn  = (x1/p)*dp_dn
110     &        + 2.0d0*c*p*z*dz_dn/((1.0d0 + z2)**3.0d0)*(1.0d0 - z2)
111      dx2_dn  = 2.0d0*C2*qb*dqb_dn
112
113      dx3a_dn = C3*dqb_dn
114      dx3b_dn = 0.5d0/x3b*(0.36d0*z*dz_dn + p*dp_dn)
115      dx3_dn  = x3b*dx3a_dn+x3a*dx3b_dn
116
117      dx4_dn  = (2.0d0*x4/p)*dp_dn
118      dx5_dn  = (2.0d0*x5/z)*dz_dn
119      dx6_dn  = (3.0d0*x6/p)*dp_dn
120      dx7_dn  = -2.0d0*es*dp_dn/(1.0d0 + es*p)**3.0d0
121
122      dxnum_dn = dx1_dn + dx2_dn + dx3_dn + dx4_dn + dx5_dn + dx6_dn
123      dx_dn    = x7*dxnum_dn + xnum*dx7_dn
124
125*     ***** dxdagr *****
126      dqba_dagr = dalpha_dagr
127      dqbb_dagr = -0.5d0*dalpha_dagr*b*(2.0d0*alpha - 1.0d0)
128     &            *qbb**3.0d0
129      dqb_dagr  = C920*(qba*dqbb_dagr + qbb*dqba_dagr)
130     &          + twthrd*dp_dagr
131      dx1_dagr  = (x1/p)*dp_dagr
132     &          + 2.0d0*c*p*z*dz_dagr/((1.0d0 + z2)**3.0d0)*(1.0d0 - z2)
133
134      dx2_dagr  = C2*2.0d0*qb*dqb_dagr
135
136      dx3a_dagr = C3*dqb_dagr
137      dx3b_dagr = 0.5d0/x3b*( 0.36d0*z*dz_dagr + p*dp_dagr)
138      dx3_dagr  = x3b*dx3a_dagr + x3a*dx3b_dagr
139
140      dx4_dagr  = (2.0d0*x4/p)*dp_dagr
141      dx5_dagr  = (2.0d0*x5/z)*dz_dagr
142      dx6_dagr  = (3.0d0*x6/p)*dp_dagr
143
144      dx7_dagr  = -2.0d0*es*dp_dagr/(1.0d0 + es*p)**3.0d0
145
146      dxnum_dagr = dx1_dagr + dx2_dagr + dx3_dagr
147     &           + dx4_dagr + dx5_dagr + dx6_dagr
148      dx_dagr    = x7*dxnum_dagr + xnum*dx7_dagr
149
150*     ***** dxdtau *****
151
152      dqb_dtau = C920*dalpha_dtau*qbb*(1.0d0
153     &         - 0.5d0*b*qba*qbb*qbb*(2.0d0*alpha - 1.0d0))
154
155      dx1_dtau = c*p*dz_dtau*2.0d0*z*(1.0d0 - z2)/
156     &           ((1.0d0 + z2)**3.0d0)
157      dx2_dtau = 2.0d0*c2*qb*dqb_dtau
158      dx3_dtau = x3*(dqb_dtau/qb
159     &         + 0.5d0*0.36d0*z*dz_dtau/(x3b*x3b))
160      dx5_dtau = 2.0d0*(x5/z)*dz_dtau
161
162      dxnum_dtau = dx1_dtau + dx2_dtau + dx3_dtau + dx5_dtau
163      dx_dtau    = x7*dxnum_dtau
164
165      Fx = 1.0d0 + kappa - kappa/(1.0d0 + x/kappa)
166
167      ex0  = Cx*n_13
168      xe   = ex0*Fx
169      nex0 = n*ex0
170
171      dFx_dx = 1.0d0/(1.0d0 + x/kappa)**2.0d0
172      dFx_dn = dFx_dx*dx_dn
173      dfdnx  = frthrd*xe + nex0*dFx_dn
174
175      dFx_dagr = dFx_dx*dx_dagr
176      dfdagrx  = nex0*dFx_dagr
177
178      dFx_dtau = dFx_dx*dx_dtau
179      dfdtaux  = nex0*dFx_dtau
180
181      return
182      end
183
184*     ************************************************
185*     *                                              *
186*     *              gen_TPSS03_restricted           *
187*     *                                              *
188*     ************************************************
189
190*    This function returns the TPSS03 exchange-correlation
191*  energy density, xce, and its derivatives with respect
192*  to n, |grad n|, tau.
193
194*
195*   Entry - n2ft3d   : number of grid points
196*           rho_in(*) :  density (nup+ndn)
197*           agr_in(*): |grad rho_in|
198*           tau_in(*): tau
199*           x_parameter: scale parameter for exchange
200*           c_parameter: scale parameter for correlation
201*
202*     Exit  - xce(n2ft3d) : TPSS03 exchange correlation energy density
203*             fn(n2ft3d)  : d(n*xce)/dn
204*             fdn(n2ft3d) : d(n*xce)/d|grad n|
205*             fdtau(n2ft3d) : d(n*xce)/dtau
206*
207
208      subroutine gen_TPSS03_restricted(n2ft3d,rho_in,agr_in,tau_in,
209     >                               x_parameter,c_parameter,
210     >                               xce,fn,fdn,fdtau)
211      implicit none
212*     ***** input *****
213      integer n2ft3d
214      real*8 rho_in(*),agr_in(*),tau_in(*)
215      real*8 x_parameter,c_parameter
216*     ***** output *****
217      real*8 xce(*),fn(*),fdn(*),fdtau(*)
218*     ***** local declarations *****
219      integer i
220      real*8 n,agr,tau
221      real*8 eup0c,fnup0c,fdnup0c
222      real*8 e0c,fn0c,fdn0c
223      real*8 Cx,P23,es
224      real*8 ex,fnx,fdnx,fdtaux
225      real*8 z,z2,z3,tmp1,tmp2,agr2,c00
226      real*8 tauW,dz2_dn,dz2_dagr,dz2_dtau,dz3_dn,dz3_dagr,dz3_dtau
227      real*8 PKZB1,dPKZB1_dn,dPKZB1_dagr,dPKZB1_dtau
228      real*8 pbec,dpbec_dn,dpbec_dagr
229      real*8 pbeupc,dpbeupc_dn,dpbeupc_dagr
230      real*8 etil,detil_dn,detil_dagr
231      real*8 PKZB2,dPKZB2_dn,dPKZB2_dagr,dPKZB2_dtau
232      real*8 revPKZB,drev_dn,drev_dagr,drev_dtau
233      real*8 revz3,drevz3_dn,drevz3_dagr,drevz3_dtau
234      real*8 ec,fnc,fdnc,fdtauc
235*     ***** constants *****
236      real*8 pi,thrd,twthrd,frthrd,fvthrd,etthrd
237      parameter (pi     =  3.14159265358979311599d0)
238      parameter (thrd   =  1.0d0/3.0d0)
239      parameter (twthrd =  2.0d0/3.0d0)
240      parameter (frthrd =  4.0d0/3.0d0)
241      parameter (fvthrd =  5.0d0/3.0d0)
242      parameter (etthrd =  8.0d0/3.0d0)
243*     ***** density cutoff parameters *****
244      real*8 tol,ETA
245      parameter (tol = 1.0d-18)
246      parameter (ETA            =      1.0d-20)
247*     ***** TPSS03 constants *****
248      real*8 kappa,mu,b,c,d,e,C920,C1,C2,C3
249      parameter (kappa =  0.8040d0)
250      parameter (mu    =  0.21951d0)
251      parameter (b     =  0.40d0)
252      parameter (c     =  1.59096d0)
253      parameter (d     =  2.8d0)
254      parameter (e     =  1.537d0)
255      parameter (C920  =  9.0d0/20.0d0)
256      parameter (C1    =  10.0d0/81.0d0)
257      parameter (C2    =  146.0d0/2025.0d0)
258      parameter (C3    = -73.0d0/405.0d0)
259
260      Cx  = (-0.75d0)*(3.0d0/pi)**thrd
261      P23 = (3.0d0*pi**2.0d0)**twthrd
262      es  = dsqrt(e)
263
264      do i=1,n2ft3d
265        n       =       rho_in(i) + ETA
266        agr     =       agr_in(i) + ETA
267        tau     = 2.0d0*tau_in(i) + ETA
268
269*       ***** TPSS03 Exchange *****
270        call nwpw_tpss03_x(pi,thrd,twthrd,frthrd,fvthrd,etthrd,
271     >                     C920,C1,C2,C3,kappa,mu,b,c,e,
272     >                     Cx,P23,es,
273     >                     n,agr,tau,
274     >                     ex,fnx,fdnx,fdtaux)
275
276*       ***** TPSS03 Correlation *****
277        call gen_PBE96_c_restricted(rho_in(i),agr_in(i),
278     >                            e0c,fn0c,fdn0c)
279
280        pbec       = e0c
281        dpbec_dn   = (fn0c - pbec)/n
282        dpbec_dagr = fdn0c/n
283
284        call gen_PBE96_c_unrestricted(0.50d0*rho_in(i),0.50d0*agr_in(i),
285     >                              eup0c,fnup0c,fdnup0c)
286
287        pbeupc       = eup0c
288        dpbeupc_dn   = (fnup0c - pbeupc)/n
289        dpbeupc_dagr = fdnup0c/n
290
291        c00  = 0.53d0
292
293        agr2 = agr*agr
294
295        tauW = 0.125d0*agr2/n
296        z    = tauW/tau
297
298        if (z .gt. 1.0d0) then
299          z2       = 1.0d0
300          dz2_dn   = 0.0d0
301          dz2_dagr = 0.0d0
302          dz2_dtau = 0.0d0
303
304          z3       = 1.0d0
305          dz3_dn   = 0.0d0
306          dz3_dagr = 0.0d0
307          dz3_dtau = 0.0d0
308        else
309          z2       =  z*z
310          dz2_dn   = -2.0d0*z2/n
311          dz2_dagr =  4.0d0*z2/agr
312          dz2_dtau = -2.0d0*z2/tau
313
314          z3       =  z2*z
315          dz3_dn   = -3.0d0*z3/n
316          dz3_dagr =  6.0d0*z3/agr
317          dz3_dtau = -3.0d0*z3/tau
318        end if
319
320        tmp1        = 1.0d0 + c00*z2
321        PKZB1       = pbec*tmp1
322        dPKZB1_dn   = dpbec_dn*tmp1     + pbec*c00*dz2_dn
323        dPKZB1_dagr = dpbec_dagr*tmp1   + pbec*c00*dz2_dagr
324        dPKZB1_dtau = pbec*c00*dz2_dtau
325
326        if (pbeupc .lt. pbec) then
327          etil       = pbec
328          detil_dn   = dpbec_dn
329          detil_dagr = dpbec_dagr
330        else
331          etil       = pbeupc
332          detil_dn   = dpbeupc_dn
333          detil_dagr = dpbeupc_dagr
334        endif
335
336        tmp1        = -(1.0d0 + c00)
337        tmp2        =  etil
338        PKZB2       =  tmp1*z2*tmp2
339        dPKZB2_dn   =  tmp1*(dz2_dn*tmp2   + z2*detil_dn)
340        dPKZB2_dagr =  tmp1*(dz2_dagr*tmp2 + z2*detil_dagr)
341        dPKZB2_dtau =  tmp1*dz2_dtau*tmp2
342
343        revPKZB   = PKZB1       + PKZB2
344        drev_dn   = dPKZB1_dn   + dPKZB2_dn
345        drev_dagr = dPKZB1_dagr + dPKZB2_dagr
346        drev_dtau = dPKZB1_dtau + dPKZB2_dtau
347
348        revz3       = 1.0d0           + d*revPKZB*z3
349        drevz3_dn   = d*(drev_dn*z3   + revPKZB*dz3_dn)
350        drevz3_dagr = d*(drev_dagr*z3 + revPKZB*dz3_dagr)
351        drevz3_dtau = d*(drev_dtau*z3 + revPKZB*dz3_dtau)
352
353        ec     = revPKZB*revz3
354        fnc    = n*(drev_dn*revz3   + revPKZB*drevz3_dn)    + ec
355        fdnc   = n*(drev_dagr*revz3 + revPKZB*drevz3_dagr)
356        fdtauc = n*(drev_dtau*revz3 + revPKZB*drevz3_dtau)
357
358        xce(i)   = x_parameter*ex     + c_parameter*ec
359        fn(i)    = x_parameter*fnx    + c_parameter*fnc
360        fdn(i)   = x_parameter*fdnx   + c_parameter*fdnc
361        fdtau(i) = x_parameter*fdtaux + c_parameter*fdtauc
362
363      end do
364      return
365      end
366
367*     ************************************************
368*     *                                              *
369*     *            gen_TPSS03_unrestricted           *
370*     *                                              *
371*     ************************************************
372
373*    This function returns the TPSS03 exchange-correlation
374*  energy density, xce, and its derivatives with respect
375*  to nup, ndn, |grad nup|, |grad ndn|, |grad n|, tauup, taudn.
376
377*
378*   Entry - n2ft3d   : number of grid points
379*           rho_in(*,2) :  density (nup and ndn)
380*           agr_in(*,3): |grad rho_in| (nup, ndn and n)
381*           tau_in(*,2): tau (nup and ndn)
382*           x_parameter: scale parameter for exchange
383*           c_parameter: scale parameter for correlation
384*
385*     Exit  - xce(n2ft3d) : TPSS03 exchange correlation energy density
386*             fn(n2ft3d,2)  : d(n*xce)/dnup, d(n*xce)/dndn
387*             fdn(n2ft3d,3) : d(n*xce)/d|grad nup|, d(n*xce)/d|grad ndn|, d(n*xce)/d|grad n|
388*             fdtau(n2ft3d,2) : d(n*xce)/dtauup, d(n*xce)/dtaudn
389*
390
391      subroutine gen_TPSS03_unrestricted(n2ft3d,rho_in,agr_in,tau_in,
392     >                               x_parameter,c_parameter,
393     >                               xce,fn,fdn,fdtau)
394      implicit none
395*     ***** input *****
396      integer n2ft3d
397      real*8 rho_in(n2ft3d,2),agr_in(n2ft3d,3),tau_in(n2ft3d,2)
398      real*8 x_parameter,c_parameter
399*     ***** output *****
400      real*8 xce(n2ft3d),fn(n2ft3d,2),fdn(n2ft3d,3),fdtau(n2ft3d,2)
401*     ***** local declarations *****
402      integer i
403      real*8 n,agr,tau
404      real*8 nup,agrup,tauup
405      real*8 ndn,agrdn,taudn
406      real*8 Cx,P23,es
407      real*8 z,z2,z3
408      real*8 ex,eupx,fnupx,fdnupx,fdtauupx,ednx,fndnx,fdndnx,fdtaudnx
409      real*8 e0c,fn0c1,fn0c2,fdn0c1,fdn0c2,fdn0c3
410      real*8 eup0c,fnup0c,fdnup0c
411      real*8 edn0c,fndn0c,fdndn0c
412      real*8 zeta,dzeta_dnup,dzeta_dndn
413      real*8 onepzeta,onemzeta,onepzeta2,onemzeta2
414      real*8 zeta2,one2mzeta2
415      real*8 n2,agr2,agraa,agrbb,agrab
416      real*8 tmp1,tmp2,denxi,denxi2,gzeta2
417      real*8 dgzeta2_dnup,dgzeta2_dndn
418      real*8 dgzeta2_dagrup,dgzeta2_dagrdn,dgzeta2_dagr
419      real*8 xi2,dxi2_dnup,dxi2_dndn
420      real*8 dxi2_dagrup,dxi2_dagrdn,dxi2_dagr
421      real*8 cz0,dcz0_dzeta,onezetap43,onezetam73
422      real*8 denczx,ddenczx_dnup,ddenczx_dndn,denczx4,denczx5,denczx8
423      real*8 czx,dczx_dnup,dczx_dndn
424      real*8 dczx_dagrup,dczx_dagrdn,dczx_dagr
425      real*8 tauW,dz2_dn,dz2_dagr,dz2_dtau,dz3_dn,dz3_dagr,dz3_dtau
426      real*8 PKZB1,dPKZB1_dnup,dPKZB1_dndn
427      real*8 dPKZB1_dagrup,dPKZB1_dagrdn,dPKZB1_dagr
428      real*8 dPKZB1_dtauup,dPKZB1_dtaudn
429      real*8 pbec,dpbec_dnup,dpbec_dndn
430      real*8 dpbec_dagrup,dpbec_dagrdn,dpbec_dagr
431      real*8 pbeupc,dpbeupc_dnup,dpbeupc_dndn
432      real*8 dpbeupc_dagrup,dpbeupc_dagrdn,dpbeupc_dagr
433      real*8 pbednc,dpbednc_dnup,dpbednc_dndn
434      real*8 dpbednc_dagrup,dpbednc_dagrdn,dpbednc_dagr
435      real*8 etilup,detilup_dnup,detilup_dndn
436      real*8 detilup_dagrup,detilup_dagrdn,detilup_dagr
437      real*8 etildn,detildn_dnup,detildn_dndn
438      real*8 detildn_dagrup,detildn_dagrdn,detildn_dagr
439      real*8 fa,fb
440      real*8 PKZB2,dPKZB2_dnup,dPKZB2_dndn
441      real*8 dPKZB2_dagrup,dPKZB2_dagrdn,dPKZB2_dagr
442      real*8 dPKZB2_dtauup,dPKZB2_dtaudn
443      real*8 revPKZB,drev_dnup,drev_dndn
444      real*8 drev_dagrup,drev_dagrdn,drev_dagr
445      real*8 drev_dtauup,drev_dtaudn
446      real*8 revz3,drevz3_dnup,drevz3_dndn
447      real*8 drevz3_dagrup,drevz3_dagrdn,drevz3_dagr
448      real*8 drevz3_dtauup,drevz3_dtaudn
449      real*8 ec,fnupc,fndnc,fdnupc,fdndnc,fdnc,fdtauupc,fdtaudnc
450*     ***** constants *****
451      real*8 pi,thrd,twthrd,frthrd,fvthrd,svthrd,etthrd
452      parameter (pi     = 3.14159265358979311599d0)
453      parameter (thrd   = 1.0d0/3.0d0)
454      parameter (twthrd = 2.0d0/3.0d0)
455      parameter (frthrd = 4.0d0/3.0d0)
456      parameter (fvthrd = 5.0d0/3.0d0)
457      parameter (svthrd = 7.0d0/3.0d0)
458      parameter (etthrd = 8.0d0/3.0d0)
459*     ***** density cutoff parameters *****
460      real*8 tol,ETA
461      parameter (tol = 1.0d-18)
462      parameter (ETA            =      1.0d-20)
463*     ***** TPSS03 constants *****
464      real*8 kappa,mu,b,c,d,e,C920,C1,C2,C3
465      parameter (kappa =  0.8040d0)
466      parameter (mu    =  0.21951d0)
467      parameter (b     =  0.40d0)
468      parameter (c     =  1.59096d0)
469      parameter (d     =  2.8d0)
470      parameter (e     =  1.537d0)
471      parameter (C920  =  9.0d0/20.0d0)
472      parameter (C1    =  10.0d0/81.0d0)
473      parameter (C2    =  146.0d0/2025.0d0)
474      parameter (C3    = -73.0d0/405.0d0)
475
476      Cx  = (-0.75d0)*(3.0d0/pi)**thrd
477      P23 = (3.0d0*pi**2.0d0)**twthrd
478      es  = dsqrt(e)
479
480      do i=1,n2ft3d
481        nup       = rho_in(i,1) + ETA
482        agrup     = agr_in(i,1) + ETA
483        tauup     = tau_in(i,1) + ETA
484        ndn       = rho_in(i,2) + ETA
485        agrdn     = agr_in(i,2) + ETA
486        taudn     = tau_in(i,2) + ETA
487
488
489*       ***** TPSS03 Exchange *****
490*       ***** UP *****
491        n   = 2.0d0*nup
492        agr = 2.0d0*agrup
493        tau = 2.0d0*tauup
494
495        call nwpw_tpss03_x(pi,thrd,twthrd,frthrd,fvthrd,etthrd,
496     >                     C920,C1,C2,C3,kappa,mu,b,c,e,
497     >                     Cx,P23,es,
498     >                     n,agr,tau,eupx,fnupx,fdnupx,fdtauupx)
499
500*       ***** DOWN *****
501        n   = 2.0d0*ndn
502        agr = 2.0d0*agrdn
503        tau = 2.0d0*taudn
504
505        call nwpw_tpss03_x(pi,thrd,twthrd,frthrd,fvthrd,etthrd,
506     >                     C920,C1,C2,C3,kappa,mu,b,c,e,
507     >                     Cx,P23,es,
508     >                     n,agr,tau,ednx,fndnx,fdndnx,fdtaudnx)
509
510        n  = nup + ndn
511
512        ex = (eupx*nup + ednx*ndn)/n
513
514*       ***** TPSS03 Correlation *****
515        agr       = agr_in(i,3) + ETA
516
517        tau  = tauup + taudn
518        n2   = n*n
519        agr2 = agr*agr
520
521        call gen_PBE96_c_full_unrestricted(rho_in(i,1),rho_in(i,2),
522     >                            agr_in(i,1),agr_in(i,2),agr_in(i,3),
523     >                            e0c,fn0c1,fn0c2,fdn0c1,fdn0c2,fdn0c3)
524
525        pbec         = e0c
526        dpbec_dnup   = (fn0c1 - pbec)/n
527        dpbec_dndn   = (fn0c2 - pbec)/n
528        dpbec_dagrup = fdn0c1/n
529        dpbec_dagrdn = fdn0c2/n
530        dpbec_dagr   = fdn0c3/n
531
532        call gen_PBE96_c_unrestricted(rho_in(i,1),agr_in(i,1),
533     >                            eup0c,fnup0c,fdnup0c)
534
535        pbeupc         = eup0c
536        dpbeupc_dnup   = (fnup0c - pbeupc)/nup
537        dpbeupc_dndn   = 0.0d0
538        dpbeupc_dagrup = fdnup0c/nup
539        dpbeupc_dagrdn = 0.0d0
540        dpbeupc_dagr   = 0.0d0
541
542        call gen_PBE96_c_unrestricted(rho_in(i,2),agr_in(i,2),
543     >                            edn0c,fndn0c,fdndn0c)
544
545        pbednc         = edn0c
546        dpbednc_dnup   = 0.0d0
547        dpbednc_dndn   = (fndn0c - pbednc)/ndn
548        dpbednc_dagrup = 0.0d0
549        dpbednc_dagrdn = fdndn0c/ndn
550        dpbednc_dagr   = 0.0d0
551
552c       if (dabs(nup - ndn) .lt. 1.0d-7) then
553c         czx         = 0.53d0
554c         dczx_dnup   = 0.0d0
555c         dczx_dndn   = 0.0d0
556c         dczx_dagrup = 0.0d0
557c         dczx_dagrdn = 0.0d0
558c         dczx_dagr   = 0.0d0
559c       else
560        zeta = (nup - ndn)/n
561
562        onepzeta  = 1.0d0 + zeta
563        onemzeta  = 1.0d0 - zeta
564        onepzeta2 = onepzeta*onepzeta
565        onemzeta2 = onemzeta*onemzeta
566
567        zeta2      = zeta*zeta
568        one2mzeta2 = 1.0d0 - zeta2
569
570        dzeta_dnup =  onemzeta/n
571        dzeta_dndn = -onepzeta/n
572
573        agraa = agrup*agrup
574        agrbb = agrdn*agrdn
575        agrab = agraa + agrbb - agr2
576
577        denxi  = 2.0d0*((3.0d0*pi*pi*n)**thrd)
578        denxi2 = denxi*denxi
579
580        gzeta2       = (agraa*onemzeta2 + agrbb*onepzeta2
581     &               + agrab*one2mzeta2)/n2
582
583        tmp1         = (-agraa*onemzeta + agrbb*onepzeta
584     &               - agrab*zeta)*2.0d0*dzeta_dnup
585        tmp1         = tmp1/n2
586        tmp2         = -2.0d0*gzeta2/n
587        dgzeta2_dnup = tmp1 + tmp2
588
589        tmp1         = (-agraa*onemzeta + agrbb*onepzeta
590     &               - agrab*zeta)*2.0d0*dzeta_dndn
591        tmp1         = tmp1/n2
592        tmp2         = -2.0d0*gzeta2/n
593        dgzeta2_dndn = tmp1 + tmp2
594
595        dgzeta2_dagrup =  2.0d0*agrup*(onemzeta2 + one2mzeta2)/n2
596        dgzeta2_dagrdn =  2.0d0*agrdn*(onepzeta2 + one2mzeta2)/n2
597        dgzeta2_dagr   = -2.0d0*agr*one2mzeta2/n2
598
599        xi2 = gzeta2/denxi2
600
601        dxi2_dnup = (dgzeta2_dnup - twthrd*gzeta2/n)/denxi2
602        dxi2_dndn = (dgzeta2_dndn - twthrd*gzeta2/n)/denxi2
603
604        dxi2_dagrup = dgzeta2_dagrup/denxi2
605        dxi2_dagrdn = dgzeta2_dagrdn/denxi2
606        dxi2_dagr   = dgzeta2_dagr/denxi2
607
608        cz0        = 0.53d0 + 0.87d0*zeta**2.0d0
609     &             + 0.50d0*zeta**4.0d0 + 2.26*zeta**6.0d0
610        dcz0_dzeta = 1.74d0*zeta + 2.0d0*zeta**3.0d0
611     &             + 13.56d0*zeta**5.0d0
612        onezetap43 = onepzeta**(-frthrd) + onemzeta**(-frthrd)
613        onezetam73 = onemzeta**(-svthrd) - onepzeta**(-svthrd)
614
615        denczx       = 1.0d0 + 0.50d0*xi2*onezetap43
616        ddenczx_dnup = 2.0d0*(denczx**3.0d0)*(dxi2_dnup*onezetap43
617     &               + xi2*frthrd*onezetam73*dzeta_dnup)
618        ddenczx_dndn = 2.0d0*(denczx**3.0d0)*(dxi2_dndn*onezetap43
619     &               + xi2*frthrd*onezetam73*dzeta_dndn)
620
621        denczx4 = denczx**4.0d0
622        denczx5 = denczx4*denczx
623        denczx8 = denczx4*denczx4
624
625        czx         = cz0/denczx4
626        dczx_dnup   = dcz0_dzeta*dzeta_dnup/denczx4
627     &              - cz0*ddenczx_dnup/denczx8
628        dczx_dndn   = dcz0_dzeta*dzeta_dndn/denczx4
629     &              - cz0*ddenczx_dndn/denczx8
630        dczx_dagrup = -2.0d0*cz0/denczx5*onezetap43*dxi2_dagrup
631        dczx_dagrdn = -2.0d0*cz0/denczx5*onezetap43*dxi2_dagrdn
632        dczx_dagr   = -2.0d0*cz0/denczx5*onezetap43*dxi2_dagr
633c       end if
634
635        tauW = 0.125d0*agr2/n
636        z    = tauW/tau
637
638        if (z .gt. 1.0d0) then
639          z2       = 1.0d0
640          dz2_dn   = 0.0d0
641          dz2_dagr = 0.0d0
642          dz2_dtau = 0.0d0
643
644          z3       = 1.0d0
645          dz3_dn   = 0.0d0
646          dz3_dagr = 0.0d0
647          dz3_dtau = 0.0d0
648        else
649          z2       =  z*z
650          dz2_dn   = -2.0d0*z2/n
651c         dz2_dnup == dz2_dndn == dz2_dn
652          dz2_dagr =  4.0d0*z2/agr
653          dz2_dtau = -2.0d0*z2/tau
654c         dz2_dtauup == dz2_dtaudn == dz2_dtau
655
656          z3       =  z2*z
657          dz3_dn   = -3.0d0*z3/n
658c         dz3_dnup == dz3_dndn == dz3_dn
659          dz3_dagr =  6.0d0*z3/agr
660          dz3_dtau = -3.0d0*z3/tau
661c         dz3_dtauup == dz3_dtaudn == dz3_dtau
662        end if
663
664        tmp1          = 1.0d0 + czx*z2
665        PKZB1         = pbec*tmp1
666        dPKZB1_dnup   = dpbec_dnup*tmp1
667     &                + pbec*(dczx_dnup*z2 + czx*dz2_dn)
668        dPKZB1_dndn   = dpbec_dndn*tmp1
669     &                + pbec*(dczx_dndn*z2 + czx*dz2_dn)
670        dPKZB1_dagrup = dpbec_dagrup*tmp1
671     &                + pbec*dczx_dagrup*z2
672        dPKZB1_dagrdn = dpbec_dagrdn*tmp1
673     &                + pbec*dczx_dagrdn*z2
674        dPKZB1_dagr   = dpbec_dagr*tmp1
675     &                + pbec*(dczx_dagr*z2 + czx*dz2_dagr)
676        dPKZB1_dtauup = pbec*czx*dz2_dtau
677        dPKZB1_dtaudn = dPKZB1_dtauup
678
679        if (pbeupc .lt. pbec) then
680          etilup         = pbec
681          detilup_dnup   = dpbec_dnup
682          detilup_dndn   = dpbec_dndn
683          detilup_dagrup = dpbec_dagrup
684          detilup_dagrdn = dpbec_dagrdn
685          detilup_dagr   = dpbec_dagr
686        else
687          etilup         = pbeupc
688          detilup_dnup   = dpbeupc_dnup
689          detilup_dndn   = dpbeupc_dndn
690          detilup_dagrup = dpbeupc_dagrup
691          detilup_dagrdn = dpbeupc_dagrdn
692          detilup_dagr   = dpbeupc_dagr
693        endif
694
695        if (pbednc .lt. pbec) then
696          etildn         = pbec
697          detildn_dnup   = dpbec_dnup
698          detildn_dndn   = dpbec_dndn
699          detildn_dagrup = dpbec_dagrup
700          detildn_dagrdn = dpbec_dagrdn
701          detildn_dagr   = dpbec_dagr
702        else
703          etildn         = pbednc
704          detildn_dnup   = dpbednc_dnup
705          detildn_dndn   = dpbednc_dndn
706          detildn_dagrup = dpbednc_dagrup
707          detildn_dagrdn = dpbednc_dagrdn
708          detildn_dagr   = dpbednc_dagr
709        endif
710
711        fa = nup/n
712        fb = ndn/n
713
714        tmp1          = -(1.0d0 + czx)
715        tmp2          =  fa*etilup + fb*etildn
716        PKZB2         =  tmp1*z2*tmp2
717        dPKZB2_dnup   = -dczx_dnup*z2*tmp2 + tmp1*(dz2_dn*tmp2
718     &                +  z2*(ndn/n2*(etilup - etildn)
719     &                +  fa*detilup_dnup + fb*detildn_dnup))
720        dPKZB2_dndn   = -dczx_dndn*z2*tmp2 + tmp1*(dz2_dn*tmp2
721     &                +  z2*(nup/n2*(etildn - etilup)
722     &                +  fb*detildn_dndn + fa*detilup_dndn))
723        dPKZB2_dagrup = -dczx_dagrup*z2*tmp2
724     &                +  tmp1*z2*(fa*detilup_dagrup + fb*detildn_dagrup)
725        dPKZB2_dagrdn = -dczx_dagrdn*z2*tmp2
726     &                +  tmp1*z2*(fb*detildn_dagrdn + fa*detilup_dagrdn)
727        dPKZB2_dagr   = -dczx_dagr*z2*tmp2 + tmp1*(dz2_dagr*tmp2
728     &                +  z2*(fa*detilup_dagr + fb*detildn_dagr))
729        dPKZB2_dtauup =  tmp1*dz2_dtau*tmp2
730        dPKZB2_dtaudn =  dPKZB2_dtauup
731
732        revPKZB     = PKZB1         + PKZB2
733        drev_dnup   = dPKZB1_dnup   + dPKZB2_dnup
734        drev_dndn   = dPKZB1_dndn   + dPKZB2_dndn
735        drev_dagrup = dPKZB1_dagrup + dPKZB2_dagrup
736        drev_dagrdn = dPKZB1_dagrdn + dPKZB2_dagrdn
737        drev_dagr   = dPKZB1_dagr   + dPKZB2_dagr
738        drev_dtauup = dPKZB1_dtauup + dPKZB2_dtauup
739        drev_dtaudn = dPKZB1_dtaudn + dPKZB2_dtaudn
740
741        revz3         = 1.0d0             + d*revPKZB*z3
742        drevz3_dnup   = d*(drev_dnup*z3   + revPKZB*dz3_dn)
743        drevz3_dndn   = d*(drev_dndn*z3   + revPKZB*dz3_dn)
744        drevz3_dagrup = d*drev_dagrup*z3
745        drevz3_dagrdn = d*drev_dagrdn*z3
746        drevz3_dagr   = d*(drev_dagr*z3   + revPKZB*dz3_dagr)
747        drevz3_dtauup = d*(drev_dtauup*z3 + revPKZB*dz3_dtau)
748        drevz3_dtaudn = d*(drev_dtaudn*z3 + revPKZB*dz3_dtau)
749
750        ec       = revPKZB*revz3
751        fnupc    = n*(drev_dnup*revz3   + revPKZB*drevz3_dnup)   + ec
752        fndnc    = n*(drev_dndn*revz3   + revPKZB*drevz3_dndn)   + ec
753        fdnupc   = n*(drev_dagrup*revz3 + revPKZB*drevz3_dagrup)
754        fdndnc   = n*(drev_dagrdn*revz3 + revPKZB*drevz3_dagrdn)
755        fdnc     = n*(drev_dagr*revz3   + revPKZB*drevz3_dagr)
756        fdtauupc = n*(drev_dtauup*revz3 + revPKZB*drevz3_dtauup)
757        fdtaudnc = n*(drev_dtaudn*revz3 + revPKZB*drevz3_dtaudn)
758
759
760        xce(i)     = x_parameter*ex       + c_parameter*ec
761        fn(i,1)    = x_parameter*fnupx    + c_parameter*fnupc
762        fn(i,2)    = x_parameter*fndnx    + c_parameter*fndnc
763        fdn(i,1)   = x_parameter*fdnupx   + c_parameter*fdnupc
764        fdn(i,2)   = x_parameter*fdndnx   + c_parameter*fdndnc
765        fdn(i,3)   = c_parameter*fdnc
766        fdtau(i,1) = x_parameter*fdtauupx + c_parameter*fdtauupc
767        fdtau(i,2) = x_parameter*fdtaudnx + c_parameter*fdtaudnc
768
769      end do
770
771      return
772      end
773