1c
2c$Id$
3c
4#include "dft2drv.fh"
5c
6c    Tao,Perdew, Staroverov, Scuseria exchange functional
7c           META GGA
8C         utilizes ingredients:
9c                              rho   -  density
10c                              delrho - gradient of density
11c                              tau (tauN)- K.S kinetic energy density
12c                              tauW - von Weiszacker kinetic energy density
13c                              tauU - uniform-gas KE density
14
15
16      Subroutine xc_ctpss03(tol_rho, cfac, lcfac, nlcfac, rho, delrho,
17     &                     nq, ipol, Ec, qwght, ldew, func,
18     &                     tau, Amat, Cmat, Mmat)
19c     References:
20c     [a] J. Tao, J.P. Perdew, V.N.Staroverov, G. Scuseria
21c         PRL 91, 146401 (2003).
22c     [b] J. Tao, J.P. Perdew, V.N.Staroverov, G. Scuseria
23c         JCP 120, 6898 (2004).
24
25      Implicit none
26c
27c
28c     Input and other parameters
29c
30      integer ipol, nq
31
32      double precision cfac
33      logical  ldew,lcfac,nlcfac
34      double precision func(*)  ! value of the functional [output]
35
36      double precision fac
37      double precision tol_rho
38c
39c     Correlation energy
40c
41      double precision Ec
42c
43c     Charge Density
44c
45      double precision rho(nq,ipol*(ipol+1)/2)
46
47c
48c     Charge Density Gradient
49c
50      double precision delrho(nq,3,ipol), gammaval, gam12
51
52c
53c     Kinetic Energy Density
54c
55      double precision tau(nq,ipol), tauN
56
57c
58c     Quadrature Weights
59c
60      double precision qwght(nq)
61c
62c     Sampling Matrices for the XC Potential
63c
64      double precision Amat(nq,ipol), Cmat(nq,*)
65      double precision Mmat(nq,*)
66
67      integer n
68      double precision rhoval
69
70c    call to the cPBE subroutine
71
72      double precision  neGGA, dneGGAdn(2), dneGGAdg(3)
73      double precision rho_t(3), delrho_t(3,2)
74      double precision  tauNA,tauNB
75      double precision  neFSP, dneFSPdn(2), dneFSPdg(3)
76      double precision delrho_A(3,2), rho_A(3)
77
78c     TPSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSs
79
80      double precision THRD, F43, F73
81      double precision zeta, ccc
82      double precision tauw
83      double precision xx2,xx3
84      double precision pbe, en
85      double precision dd
86      double precision rhoa,rhob
87      double precision pbeup, delc, erevc, erevsic,
88     A     delca,delcb,pbedown
89
90c     derivsssssssssssssssssssssssss
91
92      double precision  decggadn, dtwt2dn, dtwt3dn
93      double precision  dpbeupdn, drevdn
94      double precision  atermn, btermn
95      double precision  finaln
96      double precision  drevdg, dpbeupdg
97      double precision  dtwt2dg, dtwt3dg, decggadg
98      double precision  atermg, btermg
99      double precision  apartg, cpartg(2), finalg
100      double precision  finalgaa,finalgbb,finalgab
101      double precision drevdt, apartt, bpartt(2),finalt
102      double precision drevdta,drevdtb
103
104      double precision dzetadna, dzetadnb
105      double precision dcccdna, dcccdnb
106      double precision dcccdgaa,dcccdgbb,dcccdgab
107      double precision drevdna, drevdnb
108
109      double precision drevdgaa
110      double precision drevdgbb
111      double precision drevdgab
112      double precision dtwt3dt
113      double precision etildea,etildeb,
114     N     detiladna,detiladnb,detilbdna,detilbdnb,
115     D     detiladgaa,detiladgbb,detilbdgaa,detilbdgbb
116      double precision detiladgab,detilbdgab
117      double precision pi,fabup,fabdown
118      double precision czeta0,czeta1
119      double precision gaa,gbb,gab
120      double precision xi2,delzeta2,dencxi2zeta,
121     D     ddez2dna,ddez2dnb,dxi2dna,dxi2dnb,
122     D     ddencxi2dna,ddencxi2dnb,
123     D     ddencxi2dxi2,ddencxi2dzeta
124      double precision dxi2dgaa
125      double precision dxi2dgbb
126      double precision dxi2dgab
127      double precision denccc,onemzeta,onepzeta
128      double precision rhoval2,rhoval3
129      double precision denxi2,denxi22
130      double precision term1,term2
131      double precision ddelzeta2dna,ddelzeta2dnb
132      double precision ddelzeta2dgaa,ddelzeta2dgbb,ddelzeta2dgab
133      double precision ddencccdna,ddencccdnb
134c
135      parameter (dd = 2.8d0)
136      parameter (THRD = 1d0/3d0)
137      parameter (F43 = 4d0/3d0)
138      parameter (F73 = 7d0/3d0)
139c
140      czeta0(zeta)=
141     &  0.53d0+0.87d0*zeta**2+0.5d0*zeta**4+2.26d0*zeta**6
142      czeta1(zeta)=
143     &  2.d0*0.87d0*zeta+4d0*0.5d0*zeta**3+6d0*2.26d0*zeta**5
144      denccc(zeta,xi2)=
145     &  1.d0+0.5d0*xi2*((1.d0+zeta)**(-F43)+(1.d0-zeta)**(-F43))
146c
147      pi=acos(-1d0)
148      fac = cfac
149c
150      if (ipol.eq.1 )then
151c        ======> SPIN-RESTRICTED <======
152
153         do 22  n = 1, nq
154            if (rho(n,1).lt.tol_rho) goto 22
155
156         rhoval = rho(n,1)
157
158C   set up values to call PBE subroutine
159         rho_t(1) = rho(n,1)
160c do delrho
161         delrho_t(1,1) = delrho(n,1,1)
162         delrho_t(2,1) = delrho(n,2,1)
163         delrho_t(3,1) = delrho(n,3,1)
164         gammaval = delrho(n,1,1)*delrho(n,1,1) +
165     &        delrho(n,2,1)*delrho(n,2,1) +
166     &        delrho(n,3,1)*delrho(n,3,1)
167         gam12=dsqrt(gammaval)
168c
169c     get E_GGA[rho,gamma]
170c
171         neGGA = 0.0d0  !Ec in PBE
172         dneGGAdn(1) = 0.0d0   !Amat in PBE
173         dneGGAdg(1) = 0.0d0  !Cmat in PBE
174         dneGGAdg(2) = 0.0d0  !Cmat in PBE
175
176         call xc_cMpbe96(tol_rho,
177     &        rho_t, delrho_t,
178     &        dneGGAdn,dneGGAdg,
179     &        1, ipol, neGGA)
180         pbe = neGGA
181
182c
183c        epGGA = n*(epsilon_c^GGA) / n =cor. energy per electron
184c        epGGA= ec^LDA +H  = pbe
185
186         tauN = tau(n,1)
187
188         ccc = 0.53d0           !since zeta=0
189
190       if(sqrt(gammaval).gt.tol_rho) then
191         tauw = 0.125d0*gammaval/rhoval
192         xx2 = (tauw/tauN)**2.d0
193         xx3 = (tauw/tauN)**3.d0
194      else
195         tauw=0d0
196         xx2=0d0
197         xx3=0d0
198      endif
199         en = pbe*(1.d0 + ccc*xx2)
200c
201c    set up values to call PBE subroutine as
202c     Fully SpinPolarized system
203c
204
205         rho_A(1) = (0.5d0)*rho(n,1)   ! total   equals (1/2)n_tot
206         rho_A(2) = (0.5d0)*rho(n,1)   ! alpha   equals (1/2)n_tot
207         rho_A(3) = 0.d0               ! beta  equals zero
208         delrho_A(1,1) = (0.5d0)*delrho_t(1,1) ! nabla n_up x
209         delrho_A(2,1) = (0.5d0)*delrho_t(2,1) ! nabla n_up y
210         delrho_A(3,1) = (0.5d0)*delrho_t(3,1) ! nabla n_up z
211
212         delrho_A(1,2) = 0.d0   ! set beta gradient to zero
213         delrho_A(2,2) = 0.d0   ! set beta gradient to zero
214         delrho_A(3,2) = 0.d0   ! set beta gradient to zero
215
216         neFSP = 0.0d0  !Ec in PBE
217         dneFSPdn(1) = 0.0d0   !Amat in PBE
218         dneFSPdn(2) = 0.0d0   !Amat in PBE
219
220         dneFSPdg(1) = 0.0d0  !Cmat in PBE
221         dneFSPdg(2) = 0.0d0  !Cmat in PBE
222         dneFSPdg(3) = 0.0d0  !Cmat in PBE
223
224c
225c     get E_GGA[rho_alpha,0,gamma_alpha,0]
226c
227         call xc_cMpbe96(tol_rho, rho_A, delrho_A,
228     &        dneFSPdn,dneFSPdg,
229     &        1, 2, neFSP)
230
231         pbeup = neFSP
232
233c        functional deriv info below fffffffffffff
234
235         dtwt2dn = -2.d0*xx2/rhoval
236         dtwt3dn = -3.d0*xx3/rhoval
237         decggadn= dneGGAdn(1)
238         if(sqrt(gammaval).gt.tol_rho) then
239            dtwt2dg = 2.d0*0.125d0*tauw/(rhoval*tauN**2)
240         else
241            dtwt2dg = 0d0
242         endif
243         if(abs(taun).gt.tol_rho) then
244            dtwt3dg = 3.d0*xx2*0.125d0/(rhoval*tauN)
245         else
246            dtwt3dg = 0d0
247         endif
248         decggadg= dneGGAdg(1)
249
250       if (pbeup.lt.pbe) then
251          delc= xx2*pbe
252C  eps-tilda is eps_c
253c     erev = egga*(1-xx2)
254C  functional deriv info below fffffffffffffffff
255
256          drevdn= -pbe*dtwt2dn+decggadn*(1.d0 - xx2)
257          drevdg= -pbe*dtwt2dg+decggadg*(1.d0 - xx2)
258          drevdt= 2.d0*pbe*xx2/tauN
259       else
260          delc= xx2*pbeup
261
262C  eps-tilda is eps^FSP
263C  functional deriv info below fffffffffffffffff
264
265          dpbeupdn = 0.5d0*dneFSPdn(1)
266c  above note the .5's.  you are taking deriv wrt total density n
267c                        not deriv wrt n_up
268          dpbeupdg = 0.25d0*dneFSPdg(1)
269c  note .25 above is because you want gamma=deln_tot*deln_tot
270c
271          atermn=pbe*ccc*dtwt2dn+(1.d0+ccc*xx2)*decggadn
272          btermn=(1.d0+ccc)*(xx2*dpbeupdn+pbeup*dtwt2dn)
273          drevdn=atermn-btermn
274c
275          atermg=pbe*ccc*dtwt2dg+(1.d0+ccc*xx2)*decggadg
276          btermg=(1.d0+ccc)*(xx2*dpbeupdg+pbeup*dtwt2dg)
277          drevdg=atermg-btermg
278c
279       if(abs(taun).gt.tol_rho) then
280          drevdt=(ccc*pbe-(1.d0+ccc)*pbeup)*xx2*(-2.d0/tauN)
281       else
282          drevdt=0d0
283       endif
284       endif
285
286       delc = -(1.d0 + ccc)*delc
287       erevc = en + delc
288       erevsic = erevc*(1.d0 + dd*erevc*xx3)
289
290       if(ldew) func(n) = func(n) + rhoval*erevsic*fac
291       Ec = Ec + rhoval*erevsic*qwght(n)*fac
292
293c     derivs wrt n
294       finaln= rhoval*drevdn + erevc +
295     &      dd*(erevc*erevc*xx3 +
296     +      rhoval*(xx3*2.d0*erevc*drevdn +
297     +      erevc*erevc*dtwt3dn))
298         Amat(n,1)=Amat(n,1)+(finaln)*fac
299
300c     derivs wrt g
301         apartg=rhoval*drevdg
302         cpartg(1)=erevc*erevc*dtwt3dg
303         cpartg(2)=xx3*2.d0*erevc*drevdg
304         finalg=apartg+rhoval*dd*( cpartg(1)+cpartg(2) )
305         Cmat(n,D1_GAA)=Cmat(n,D1_GAA)+ 2d0*finalg*fac
306
307c     derivs wrt t
308         apartt=rhoval*drevdt
309       if(abs(taun).gt.tol_rho) then
310          bpartt(1)=-erevc*erevc*xx3*3.d0/tauN
311       else
312          bpartt(1)=0d0
313       endif
314         bpartt(2)=xx3*2.d0*erevc*drevdt
315         finalt=apartt+dd*rhoval*(bpartt(1)+bpartt(2))
316         Mmat(n,1)=Mmat(n,1)+0.5d0*finalt*fac
317
31822     continue
319
320      else !ipol=2
321c
322c        ======> SPIN-UNRESTRICTED <======
323c
324         do 20 n = 1, nq
325c
326         if (rho(n,1).lt.tol_rho) goto 20
327
328         rhoval = rho(n,1)
329
330         rho_t(1) = rho(n,1)
331         rho_t(2) = rho(n,2)
332         rho_t(3) = rho(n,3)
333         delrho_t(1,1) = delrho(n,1,1)
334         delrho_t(2,1) = delrho(n,2,1)
335         delrho_t(3,1) = delrho(n,3,1)
336         delrho_t(1,2) = delrho(n,1,2)
337         delrho_t(2,2) = delrho(n,2,2)
338         delrho_t(3,2) = delrho(n,3,2)
339c
340         neGGA = 0.0d0  !Ec in PBE
341         dneGGAdn(1) = 0.0d0   !Amat in PBE (n,1)
342         dneGGAdn(2) = 0.0d0   !Amat in PBE (n,2)
343         dneGGAdg(1) = 0.0d0  !Cmat in PBE--aa
344         dneGGAdg(2) = 0.0d0  !Cmat in PBE--ab
345         dneGGAdg(3) = 0.0d0  !Cmat in PBE--bb
346c
347c     get E_GGA[rho,gamma]
348c
349         call xc_cMpbe96(tol_rho,
350     &        rho_t, delrho_t,
351     &        dneGGAdn,dneGGAdg,
352     &        1, ipol, neGGA)
353         pbe = neGGA
354c
355c        epGGA = (epsilon_c^GGA)  =cor. energy per electron
356c        epGGA= ec^LDA +H  = pbe
357c
358         gammaval = delrho(n,1,1)*delrho(n,1,1) +
359     &        delrho(n,1,2)*delrho(n,1,2) +
360     &        delrho(n,2,1)*delrho(n,2,1) +
361     &        delrho(n,2,2)*delrho(n,2,2) +
362     &        delrho(n,3,1)*delrho(n,3,1) +
363     &        delrho(n,3,2)*delrho(n,3,2) +
364     &        2.d0*(delrho(n,1,1)*delrho(n,1,2) +
365     &        delrho(n,2,1)*delrho(n,2,2) +
366     &        delrho(n,3,1)*delrho(n,3,2))
367         gam12=dsqrt(gammaval)
368         tauNa = tau(n,1)
369         tauNb = tau(n,2)
370         taun = tauna+taunb
371         rhoa=rho(n,2)
372         rhob=rho(n,3)
373c
374c        Check for small densities (H atom case as well)
375c
376         if ((rhoa.lt.tol_rho).or.
377     &          (rhob.lt.tol_rho)) goto 20
378c
379         dcccdna=0.d0
380         dcccdnb=0.d0
381         dcccdgaa=0.d0
382         dcccdgbb=0.d0
383         dcccdgab=0.d0
384         ccc=0d0
385c
386         if(rhoval.gt.tol_rho) then
387c
388            zeta=(rhoa-rhob)/rhoval
389c
390            onepzeta = 1.d0+zeta
391            onemzeta = 1.d0-zeta
392c
393            gaa = delrho(n,1,1)*delrho(n,1,1) +
394     &            delrho(n,2,1)*delrho(n,2,1) +
395     &            delrho(n,3,1)*delrho(n,3,1)
396            gbb = delrho(n,1,2)*delrho(n,1,2) +
397     &            delrho(n,2,2)*delrho(n,2,2) +
398     &            delrho(n,3,2)*delrho(n,3,2)
399            gab = delrho(n,1,1)*delrho(n,1,2) +
400     &            delrho(n,2,1)*delrho(n,2,2) +
401     &            delrho(n,3,1)*delrho(n,3,2)
402c
403            rhoval2 = rhoval*rhoval
404            rhoval3 = rhoval*rhoval*rhoval
405            denxi2 = 2.d0*(3.d0*pi*pi*rhoval)**(1.d0/3.d0)
406            denxi22 = denxi2*denxi2
407            delzeta2 = (gaa*onemzeta*onemzeta + gbb*onepzeta*onepzeta
408     &                  -2.d0*gab*onemzeta*onepzeta)/rhoval2
409c
410            dzetadna =  onemzeta/rhoval
411            dzetadnb = -onepzeta/rhoval
412c
413            term1=-2.d0*gaa*onemzeta+2.d0*gbb*onepzeta+4.d0*gab*zeta
414            term1=(term1/rhoval2)*dzetadna
415            term2=-2.d0*delzeta2/rhoval
416            ddelzeta2dna = term1 + term2
417c
418            term1=-2.d0*gaa*onemzeta+2.d0*gbb*onepzeta+4.d0*gab*zeta
419            term1=(term1/rhoval2)*dzetadnb
420            term2=-2.d0*delzeta2/rhoval
421            ddelzeta2dnb = term1 + term2
422c
423            ddelzeta2dgaa =  onemzeta*onemzeta/rhoval2
424            ddelzeta2dgbb =  onepzeta*onepzeta/rhoval2
425            ddelzeta2dgab =  -2.d0*onepzeta*onemzeta/rhoval2
426c
427            xi2=delzeta2/denxi22
428            dxi2dna=(ddelzeta2dna -(2.d0/3d0)*delzeta2/rhoval)/denxi22
429            dxi2dnb=(ddelzeta2dnb -(2.d0/3d0)*delzeta2/rhoval)/denxi22
430c
431            dxi2dgaa=onemzeta*onemzeta/rhoval2/denxi22
432	    dxi2dgbb=onepzeta*onepzeta/rhoval2/denxi22
433            dxi2dgab=-2.d0*onepzeta*onemzeta/rhoval2/denxi22
434c
435            ccc=czeta0(zeta)/(denccc(zeta,xi2)**4)
436            ddencccdna=2.d0*(denccc(zeta,xi2)**3)*
437     &       (dxi2dna*(onepzeta**(-F43) + onemzeta**(-F43))
438     &        + xi2*F43*(onemzeta**(-F73) - onepzeta**(-F73))*dzetadna)
439            ddencccdnb=2.d0*(denccc(zeta,xi2)**3)*
440     &       (dxi2dnb*(onepzeta**(-F43) + onemzeta**(-F43))
441     &        + xi2*F43*(onemzeta**(-F73) - onepzeta**(-F73))*dzetadnb)
442c
443            dcccdna=(czeta1(zeta)*dzetadna/denccc(zeta,xi2)**4) -
444     &              (czeta0(zeta)*ddencccdna/(denccc(zeta,xi2)**8))
445            dcccdnb=(czeta1(zeta)*dzetadnb/denccc(zeta,xi2)**4) -
446     &              (czeta0(zeta)*ddencccdnb/(denccc(zeta,xi2)**8))
447c
448            dcccdgaa=-(czeta0(zeta)/(denccc(zeta,xi2)**8))*
449     &       2.d0*(denccc(zeta,xi2)**3)*
450     &       dxi2dgaa*(onepzeta**(-F43) + onemzeta**(-F43))
451            dcccdgbb=-(czeta0(zeta)/(denccc(zeta,xi2)**8))*
452     &       2.d0*(denccc(zeta,xi2)**3)*
453     &       dxi2dgbb*(onepzeta**(-F43) + onemzeta**(-F43))
454            dcccdgab=-(czeta0(zeta)/(denccc(zeta,xi2)**8))*
455     &       2.d0*(denccc(zeta,xi2)**3)*
456     &       dxi2dgab*(onepzeta**(-F43) + onemzeta**(-F43))
457         endif
458c
459         tauw = 0.125d0*gammaval/rhoval
460         if(abs(tauw).gt.tol_rho) then
461c
462         xx2 = (tauw/tauN)**2.d0
463         xx3 = (tauw/tauN)**3.d0
464         dtwt2dn = -2.d0*xx2/rhoval
465         dtwt3dn = -3.d0*xx3/rhoval
466         dtwt3dt = -3.d0*xx3/taun
467         dtwt2dg = 2.d0*0.125d0*tauw/(rhoval*tauN**2)
468         dtwt3dg = 3.d0*xx2*0.125d0/(rhoval*tauN)
469      else
470         tauw=0d0
471         xx2=0d0
472         xx3=0d0
473         dtwt2dn = 0d0
474         dtwt3dn = 0d0
475         dtwt3dt = 0d0
476         dtwt2dg = 0d0
477         dtwt3dg = 0d0
478
479      endif
480c
481         en = pbe*(1.d0 + ccc*xx2)
482c
483c     Alpha bit
484c     set up values to call PBE subroutine as
485c     Fully SpinPolarized system for alpha spin
486c     to get E_GGA[rho_alpha,0,gamma_alpha,0]
487c
488         rho_A(1) = rhoa
489         rho_A(2) = rhoa
490         rho_A(3) = 0.d0               ! beta  equals zero
491         delrho_A(1,1) = delrho_t(1,1) ! nabla n_up x
492         delrho_A(2,1) = delrho_t(2,1) ! nabla n_up y
493         delrho_A(3,1) = delrho_t(3,1) ! nabla n_up z
494         delrho_A(1,2) = 0.d0   ! set beta gradient to zero
495         delrho_A(2,2) = 0.d0   ! set beta gradient to zero
496         delrho_A(3,2) = 0.d0   ! set beta gradient to zero
497
498         neFSP = 0.0d0  !Ec in PBE
499         dneFSPdn(1) = 0.0d0   !Amat in PBE
500         dneFSPdn(2) = 0.0d0   !Amat in PBE
501
502         dneFSPdg(1) = 0.0d0  !Cmat in PBE
503         dneFSPdg(2) = 0.0d0  !Cmat in PBE
504         dneFSPdg(3) = 0.0d0  !Cmat in PBE
505c
506         call xc_cMpbe96(tol_rho, rho_A, delrho_A,
507     &        dneFSPdn, dneFSPdg,  1, 2, neFSP)
508
509         pbeup = neFSP
510
511c        functional deriv info below fffffffffffff
512       if (pbeup.lt.pbe) then
513          etildea    = pbe
514          detiladna  = dneggadn(1)
515          detiladnb  = dneggadn(2)
516          detiladgaa = dneggadg(D1_GAA)
517          detiladgbb = dneggadg(D1_GBB)
518          detiladgab = dneggadg(D1_GAB)
519       else
520          etildea    = pbeup
521          detiladna  = dneFSPdn(1)
522          detiladnb  = 0.d0
523          detiladgaa = dneFSPdg(D1_GAA)
524          detiladgbb = 0.d0
525          detiladgab = 0.d0
526       endif
527
528c     n_sigma/n_total factor
529       fabup=rhoa/rhoval
530
531       delc= xx2*etildea
532       delca = -(1.d0 + ccc)*fabup*delc
533       erevc = en + delca
534c
535c     Beta bit
536c     set up values to call PBE subroutine as
537c     Fully SpinPolarized system for beta spin
538c     to get E_GGA[rho_beta,0,gamma_beta,0]
539c
540         rho_A(1) = rhob
541         rho_A(2) = rhob
542         rho_A(3) = 0.d0               ! beta  equals zero
543         delrho_A(1,1) = delrho_t(1,2) ! nabla n_up x
544         delrho_A(2,1) = delrho_t(2,2) ! nabla n_up y
545         delrho_A(3,1) = delrho_t(3,2) ! nabla n_up z
546         delrho_A(1,2) = 0.d0   ! set beta gradient to zero
547         delrho_A(2,2) = 0.d0   ! set beta gradient to zero
548         delrho_A(3,2) = 0.d0   ! set beta gradient to zero
549
550         neFSP = 0.0d0  !Ec in PBE
551         dneFSPdn(1) = 0.0d0   !Amat in PBE
552         dneFSPdn(2) = 0.0d0   !Amat in PBE
553
554         dneFSPdg(1) = 0.0d0  !Cmat in PBE
555         dneFSPdg(2) = 0.0d0  !Cmat in PBE
556         dneFSPdg(3) = 0.0d0  !Cmat in PBE
557c
558         call xc_cMpbe96(tol_rho, rho_A, delrho_A,
559     &        dneFSPdn,dneFSPdg, 1, 2, neFSP)
560
561         pbedown = neFSP
562
563c        functional deriv info below fffffffffffff
564
565       if (pbedown.lt.pbe) then
566          etildeb=pbe
567          detilbdna  = dneggadn(1)
568          detilbdnb  = dneggadn(2)
569          detilbdgaa = dneggadg(D1_GAA)
570          detilbdgbb = dneggadg(D1_GBB)
571          detilbdgab = dneggadg(D1_GAB)
572       else
573          etildeb    = pbedown
574          detilbdna  = 0.d0
575          detilbdnb  = dneFSPdn(1)
576          detilbdgaa = 0.d0
577          detilbdgbb = dneFSPdg(D1_GAA)
578          detilbdgab = 0.d0
579       endif
580c
581c     n_sigma/n_total factor
582       fabdown=rhob/rhoval
583
584       delc= xx2*etildeb
585       delcb = -(1.d0 + ccc)*fabdown*delc
586       erevc = erevc + delcb
587
588       erevsic = erevc*(1.d0 + dd*erevc*xx3)
589
590       if(ldew) func(n) = func(n) + rhoval*erevsic*fac
591       Ec = Ec + rhoval*erevsic*qwght(n)*fac
592c
593c na
594       atermn=pbe*ccc*dtwt2dn+(1.d0+ccc*xx2)*dneggadn(1)+
595     &      pbe*xx2*dcccdna
596       btermn=(1.d0+ccc)*(
597     &      dtwt2dn*(fabup*etildea+fabdown*etildeb) +
598     +      xx2*( (etildea - etildeb)*rhob/(rhoval*rhoval) +
599     +      fabup*detiladna + fabdown*detilbdna) ) +
600     Z      xx2*(fabup*etildea+fabdown*etildeb)*dcccdna
601       drevdna=atermn-btermn
602c
603c nb
604       atermn=pbe*ccc*dtwt2dn+(1.d0+ccc*xx2)*dneggadn(2)+
605     &      pbe*xx2*dcccdnb
606       btermn=(1.d0+ccc)*(
607     &      dtwt2dn*(fabup*etildea+fabdown*etildeb) +
608     &      xx2*( (etildeb-etildea)*rhoa/(rhoval*rhoval) +
609     &      fabup*detilbdna+fabdown*detilbdnb) )+
610     &      xx2*(fabup*etildea+fabdown*etildeb)*dcccdnb
611       drevdnb=atermn-btermn
612c
613c gaa
614       atermg=pbe*ccc*dtwt2dg+(1.d0+ccc*xx2)*dneggadg(D1_GAA)+
615     &      pbe*xx2*dcccdgaa
616       btermg=(1.d0+ccc)*
617     &      (xx2*(fabup*detiladgaa+fabdown*detilbdgaa)+
618     &      (etildea*fabup+etildeb*fabdown)*dtwt2dg)+
619     &      xx2*(etildea*fabup+etildeb*fabdown)*dcccdgaa
620       drevdgaa=atermg-btermg
621c
622c gbb
623       atermg=pbe*ccc*dtwt2dg+(1.d0+ccc*xx2)*dneggadg(D1_GBB)+
624     &      pbe*xx2*dcccdgbb
625       btermg=(1.d0+ccc)*
626     &      (xx2*(fabup*detiladgbb+fabdown*detilbdgbb)+
627     &      (etildea*fabup+etildeb*fabdown)*dtwt2dg)+
628     &      xx2*(etildea*fabup+etildeb*fabdown)*dcccdgbb
629       drevdgbb=atermg-btermg
630c
631c gab
632       atermg=pbe*ccc*2.d0*dtwt2dg+(1.d0+ccc*xx2)*dneggadg(D1_GAB)+
633     &      pbe*xx2*dcccdgab
634       btermg=(1.d0+ccc)*
635     &      (xx2*(fabup*detiladgab+fabdown*detilbdgab)+
636     &      (etildea*fabup+etildeb*fabdown)*2.d0*dtwt2dg)+
637     &      xx2*(etildea*fabup+etildeb*fabdown)*dcccdgab
638       drevdgab=atermg-btermg
639c
640c     ta,tb
641       if(abs(taun).gt.tol_rho) then
642       drevdta=-2.d0*xx2/tauN*
643     *(ccc*pbe-(1.d0+ccc)*(fabup*etildea+fabdown*etildeb))
644       drevdtb=-2d0*xx2/tauN*
645     *      (ccc*pbe-(1.d0+ccc)*(fabup*etildea+fabdown*etildeb))
646      else
647         drevdta=0d0
648         drevdtb=0d0
649       endif
650c
651c derivs wrt na
652         finaln= rhoval*drevdna + erevc +
653     &      dd*(erevc*erevc*xx3 + rhoval*(xx3*2.d0*erevc*drevdna
654     &     + erevc*erevc*dtwt3dn))
655         Amat(n,1)=Amat(n,1)+finaln*fac
656c
657c derivs wrt nb
658         finaln= rhoval*drevdnb + erevc +
659     &      dd*(erevc*erevc*xx3 + rhoval*(xx3*2.d0*erevc*drevdnb
660     &     + erevc*erevc*dtwt3dn))
661         Amat(n,2)=Amat(n,2)+finaln*fac
662c
663c derivs wrt gaa
664         apartg=rhoval*drevdgaa
665         cpartg(1)=erevc*erevc*dtwt3dg
666         cpartg(2)=xx3*2.d0*erevc*drevdgaa
667         finalgaa=apartg+rhoval*dd*(cpartg(1)+cpartg(2))
668         Cmat(n,D1_GAA)=Cmat(n,D1_GAA)+finalgaa*fac
669c
670c derivs wrt gbb
671         apartg=rhoval*drevdgbb
672         cpartg(1)=erevc*erevc*dtwt3dg
673         cpartg(2)=xx3*2.d0*erevc*drevdgbb
674         finalgbb=apartg+rhoval*dd*(cpartg(1)+cpartg(2))
675         Cmat(n,D1_GBB)=Cmat(n,D1_GBB)+finalgbb*fac
676c
677c derivs wrt gab
678         apartg=rhoval*drevdgab
679         cpartg(1)=erevc*erevc*2.d0*dtwt3dg
680         cpartg(2)=xx3*2.d0*erevc*drevdgab
681         finalgab=apartg+rhoval*dd*(cpartg(1)+cpartg(2))
682         Cmat(n,D1_GAB)=Cmat(n,D1_GAB)+finalgab*fac
683c
684c derivs wrt ta
685         apartt=rhoval*drevdta
686         bpartt(1)=erevc*erevc*dtwt3dt
687         bpartt(2)=xx3*2.d0*erevc*drevdta
688         finalt=apartt+dd*rhoval*(bpartt(1)+bpartt(2))
689         Mmat(n,1)=Mmat(n,1)+0.5d0*finalt*fac
690c
691c derivs wrt tb
692         apartt=rhoval*drevdtb
693         bpartt(1)=erevc*erevc*dtwt3dt
694         bpartt(2)=xx3*2.d0*erevc*drevdtb
695         finalt=apartt+dd*rhoval*(bpartt(1)+bpartt(2))
696         Mmat(n,2)=Mmat(n,2)+0.5d0*finalt*fac
697c
69820     continue
699
700      endif  !end ipol=2 case
701
702      return
703      end
704c
705      Subroutine xc_ctpss03_d2()
706         call errquit(' ctpss03: d2 not coded ',0,0)
707      return
708      end
709