1#include "dft2drv.fh"
2!#define PRINTA 1
3      Subroutine xc_cpkzb99(tol_rho, cfac, lcfac, nlcfac, rho, delrho,
4     &                     nq, ipol, Ec, qwght, ldew, func,
5     &                     tau, Amat, Cmat, Mmat)
6
7
8c
9c$Id$
10c
11
12c     References:
13c     [a] J.P. Perdew, S. Kurth, A. Zupan and P. Blaha,
14c         PRL 82, 2544 (1999).
15
16      Implicit none
17c
18c
19c     Input and other parameters
20c
21
22      integer ipol, nq
23
24      double precision cfac
25      logical lcfac, nlcfac
26      logical  ldew
27      double precision func(*)
28
29      double precision fac
30      double precision tol_rho
31c
32c     Correlation energy
33c
34      double precision Ec
35c
36c     Charge Density
37c
38      double precision rho(nq,ipol*(ipol+1)/2)
39
40c
41c     Charge Density Gradient
42c
43      double precision delrho(nq,3,ipol), gammaval
44
45c
46c     Kinetic Energy Density
47c
48      double precision tau(nq,ipol)
49
50c
51c     Quadrature Weights
52c
53      double precision qwght(nq)
54c
55c     Sampling Matrices for the XC Potential
56c
57      double precision Amat(nq,ipol), Cmat(nq,*)
58      double precision Mmat(nq,*)
59
60      integer n
61      double precision rhoval,rhoa,rhob
62
63c    first sigma term
64      double precision  taun
65      double precision  ccc
66      parameter (ccc = 0.53d0) !cpkzb empirical parameter
67
68c   Second call to the cPBE subroutine
69
70      double precision  neGGA, dneGGAdn(2), dneGGAdg(3)
71      double precision rho_t(3), delrho_t(3,2)
72      double precision  tauNA,tauNB
73c
74      double precision gam12,pbe,tauw,xx2,en,
75     ,     tauwa,tauwb,xx2a,xx2b,dtwat2dg,dtwat2dn,
76     ,     dtwbt2dg,dtwbt2dn
77      double precision pbeup,dtwt2dn,decggadn,dtwt2dg,
78     ,     delc,decggadg,drevdn,drevdg,drevdt,
79     ,     dpbeupdn,dpbeupdg,atermn,btermn,atermg,btermg,
80     ,     erevc,finaln,apartg,finalg,apartt,finalt
81c
82      double precision  neFSP, dneFSPdn(2), dneFSPdg(3)
83c
84      double precision drevdna,drevdnb,drevdgaa,drevdgbb,
85     A     drevdta,drevdtb,finalgbb
86      double precision delca,delcb,
87     A     detiladga,detiladgb,detilbdga,detilbdgb,
88     A     detiladna,detiladnb,detilbdna,detilbdnb
89      double precision etildea,etildeb,gaa,gbb,gab
90      double precision fabup,fabdown
91      double precision delrho_A(3,2), rho_A(3)
92c
93      double precision xx1,xx1a,xx1b,pbedown
94      double precision tauwplus,taunplus,rhoval2
95      double precision dxx1dna,dxx1dnb
96      double precision dxx1adna,dxx1bdnb
97      double precision dxx1dgaa,dxx1dgbb
98      double precision dxx1adgaa,dxx1bdgbb
99      double precision drevdgab
100      double precision dxx1dta,dxx1dtb
101      double precision dxx1adta,dxx1bdtb
102      double precision finalna,finalnb
103      double precision finalgaa,finalgab
104      double precision rhoa2,rhob2
105      double precision detiladgaa,detiladgbb
106      double precision detilbdgaa,detilbdgbb
107c
108      fac = cfac
109      if (ipol.eq.1 )then
110c        ======> SPIN-RESTRICTED <======
111         do 12  n = 1, nq
112         if (rho(n,1).lt.tol_rho) goto 12
113
114         rhoval = rho(n,1)
115
116C   set up values to call PBE subroutine
117         rho_t(1) = rho(n,1)
118c do delrho
119         delrho_t(1,1) = delrho(n,1,1)
120         delrho_t(2,1) = delrho(n,2,1)
121         delrho_t(3,1) = delrho(n,3,1)
122         gammaval = delrho(n,1,1)*delrho(n,1,1) +
123     &              delrho(n,2,1)*delrho(n,2,1) +
124     &              delrho(n,3,1)*delrho(n,3,1)
125         gam12=dsqrt(gammaval)
126c
127c     get E_GGA[rho,gamma]
128c
129         neGGA = 0.0d0  !Ec in PBE
130         dneGGAdn(1) = 0.0d0   !Amat in PBE
131         dneGGAdg(1) = 0.0d0  !Cmat in PBE
132         dneGGAdg(2) = 0.0d0  !Cmat in PBE
133
134         call xc_cMpbe96(tol_rho,
135     &        rho_t, delrho_t,
136     &        dneGGAdn,dneGGAdg,
137     &        1, ipol, neGGA)
138         pbe = neGGA
139
140         tauN = tau(n,1)
141         tauw = 0.125d0*gammaval/rhoval
142         xx2 = (tauw/tauN)**2.d0
143         en = pbe*(1.d0 + ccc*xx2)
144c
145c    set up values to call PBE subroutine as
146c     Fully SpinPolarized system
147c
148
149         rho_A(1) = (0.5d0)*rho(n,1)   ! total   equals (1/2)n_tot
150         rho_A(2) = (0.5d0)*rho(n,1)   ! alpha   equals (1/2)n_tot
151         rho_A(3) = 0.d0               ! beta  equals zero
152         delrho_A(1,1) = (0.5d0)*delrho_t(1,1) ! nabla n_up x
153         delrho_A(2,1) = (0.5d0)*delrho_t(2,1) ! nabla n_up y
154         delrho_A(3,1) = (0.5d0)*delrho_t(3,1) ! nabla n_up z
155
156         delrho_A(1,2) = 0.d0   ! set beta gradient to zero
157         delrho_A(2,2) = 0.d0   ! set beta gradient to zero
158         delrho_A(3,2) = 0.d0   ! set beta gradient to zero
159
160         neFSP = 0.0d0  !Ec in PBE
161         dneFSPdn(1) = 0.0d0   !Amat in PBE
162         dneFSPdn(2) = 0.0d0   !Amat in PBE
163
164         dneFSPdg(1) = 0.0d0  !Cmat in PBE
165         dneFSPdg(2) = 0.0d0  !Cmat in PBE
166         dneFSPdg(3) = 0.0d0  !Cmat in PBE
167
168c
169c     get E_GGA[rho_alpha,0,gamma_alpha,0]
170c
171         call xc_cMpbe96(tol_rho, rho_A, delrho_A,
172     &        dneFSPdn,dneFSPdg, 1, 2, neFSP)
173         pbeup = neFSP
174
175c        functional deriv info below fffffffffffff
176         dtwt2dn = -2.d0*xx2/rhoval
177         decggadn= dneGGAdn(1)
178         dtwt2dg = 2.d0*0.125d0*tauw/(rhoval*tauN**2)
179         decggadg= dneGGAdg(1)
180         delc= xx2*pbeup
181
182C  eps-tilda is eps^FSP
183C  functional deriv info below fffffffffffffffff
184
185        dpbeupdn = 0.5d0*dneFSPdn(1)
186c  above note the .5's.  you are taking deriv wrt total density n
187c                        not deriv wrt n_up
188        dpbeupdg = 0.25d0*dneFSPdg(1)
189c  note .25 above is because you want gamma=deln_tot*deln_tot
190
191
192        atermn=pbe*ccc*dtwt2dn+(1.d0+ccc*xx2)*decggadn
193        btermn=(1.d0+ccc)*(xx2*dpbeupdn + pbeup*dtwt2dn)
194        drevdn=atermn - btermn
195
196        atermg=pbe*ccc*dtwt2dg+(1.d0+ccc*xx2)*decggadg
197        btermg=(1.d0+ccc)*(xx2*dpbeupdg+pbeup*dtwt2dg)
198        drevdg=atermg-btermg
199        drevdt=(ccc*pbe-(1.d0+ccc)*pbeup)*xx2*(-2.d0/tauN)
200
201       delc = -(1.d0 + ccc)*delc
202       erevc = en + delc
203
204       if(ldew) func(n) = func(n) + rhoval*erevc*fac
205       Ec = Ec + rhoval*erevc*qwght(n)*fac
206
207c     derivs wrt n
208         finaln= rhoval*drevdn + erevc
209         Amat(n,1)=Amat(n,1)+(finaln)*fac
210
211c     derivs wrt g
212         apartg=rhoval*drevdg
213         finalg=apartg
214         Cmat(n,D1_GAA)=Cmat(n,D1_GAA)+ 2.d0*finalg*fac
215
216c     derivs wrt t
217         apartt=rhoval*drevdt
218         finalt=apartt
219         Mmat(n,1)=Mmat(n,1)+0.5d0*finalt*fac
220
22112     continue
222c
223c     open-shell
224c
225       else   !ipol=2 and do alpha beta cases
226         do 20 n = 1, nq
227c
228         if (rho(n,1).lt.tol_rho) goto 20
229c
230         rhoval = rho(n,1)
231         rhoval2 = rhoval*rhoval
232c
233         rho_t(1) = rho(n,1)
234         rho_t(2) = rho(n,2)
235         rho_t(3) = rho(n,3)
236         delrho_t(1,1) = delrho(n,1,1)
237         delrho_t(2,1) = delrho(n,2,1)
238         delrho_t(3,1) = delrho(n,3,1)
239         delrho_t(1,2) = delrho(n,1,2)
240         delrho_t(2,2) = delrho(n,2,2)
241         delrho_t(3,2) = delrho(n,3,2)
242
243         neGGA = 0.0d0  !Ec in PBE
244         dneGGAdn(1) = 0.0d0   !Amat in PBE (n,1)
245         dneGGAdn(2) = 0.0d0   !Amat in PBE (n,2)
246         dneGGAdg(1) = 0.0d0  !Cmat in PBE--aa
247         dneGGAdg(2) = 0.0d0  !Cmat in PBE--ab
248         dneGGAdg(3) = 0.0d0  !Cmat in PBE--bb
249c
250c     get E_GGA[rho,gamma]
251c
252         call xc_cMpbe96(tol_rho,
253     &        rho_t, delrho_t,
254     &        dneGGAdn,dneGGAdg,
255     &        1, ipol, neGGA)
256         pbe = neGGA
257c
258c        epGGA = (epsilon_c^GGA)  =cor. energy per electron
259c        epGGA= ec^LDA +H  = pbe
260c
261         gaa = delrho(n,1,1)*delrho(n,1,1) +
262     &         delrho(n,2,1)*delrho(n,2,1) +
263     &         delrho(n,3,1)*delrho(n,3,1)
264         gbb = delrho(n,1,2)*delrho(n,1,2) +
265     &         delrho(n,2,2)*delrho(n,2,2) +
266     &         delrho(n,3,2)*delrho(n,3,2)
267         gab = delrho(n,1,1)*delrho(n,1,2) +
268     &         delrho(n,2,1)*delrho(n,2,2) +
269     &         delrho(n,3,1)*delrho(n,3,2)
270c
271         rhoa=rho(n,2)
272         rhoa2 = rhoa*rhoa
273         rhob=rho(n,3)
274         rhob2 = rhob*rhob
275c
276c        Check for small densities (H atom case as well)
277c
278         if ((rhoa.lt.tol_rho).or.
279     &          (rhob.lt.tol_rho)) goto 20
280c
281         tauwa = 0.125d0*gaa/rhoa
282         tauwb = 0.125d0*gbb/rhob
283c
284         tauna = tau(n,1)
285         taunb = tau(n,2)
286c
287         tauw = tauwa+tauwb
288         taun = tauna+taunb
289c
290         xx1 = tauw/taun
291         xx2 = xx1*xx1
292c
293         xx1a = tauwa/tauna
294         xx2a = xx1a*xx1a
295c
296         xx1b = tauwb/taunb
297         xx2b = xx1b*xx1b
298c
299         en = pbe*(1.d0 + ccc*xx2)
300c
301c     Alpha bit
302c    set up values to call PBE subroutine as
303c     Fully SpinPolarized system for alpha spin
304c     to get E_GGA[rho_alpha,0,gamma_alpha,0]
305c
306         rho_A(1) = rhoa
307         rho_A(2) = rhoa
308         rho_A(3) = 0.d0               ! beta  equals zero
309         delrho_A(1,1) = delrho_t(1,1) ! nabla n_up x
310         delrho_A(2,1) = delrho_t(2,1) ! nabla n_up y
311         delrho_A(3,1) = delrho_t(3,1) ! nabla n_up z
312         delrho_A(1,2) = 0.d0   ! set beta gradient to zero
313         delrho_A(2,2) = 0.d0   ! set beta gradient to zero
314         delrho_A(3,2) = 0.d0   ! set beta gradient to zero
315
316         neFSP = 0.0d0  !Ec in PBE
317         dneFSPdn(1) = 0.0d0   !Amat in PBE
318         dneFSPdn(2) = 0.0d0   !Amat in PBE
319
320         dneFSPdg(1) = 0.0d0  !Cmat in PBE
321         dneFSPdg(2) = 0.0d0  !Cmat in PBE
322         dneFSPdg(3) = 0.0d0  !Cmat in PBE
323c
324         call xc_cMpbe96(tol_rho, rho_A, delrho_A,
325     &        dneFSPdn,dneFSPdg, 1, 2, neFSP)
326         pbeup = neFSP
327c
328c        functional deriv info below fffffffffffff
329         etildea= pbeup
330         detiladna = dneFSPdn(1)
331         detiladnb = 0d0
332         detiladgaa = dneFSPdg(D1_GAA)
333         detiladgbb = 0d0
334c
335c     n_sigma/n_total factor
336       fabup=rhoa/rhoval
337       delca = -(1.d0 + ccc)*fabup*xx2a*etildea
338       erevc = en + delca
339c
340c     Beta bit
341c    set up values to call PBE subroutine as
342c     Fully SpinPolarized system for beta spin
343c     to get E_GGA[rho_beta,0,gamma_beta,0]
344c
345       rho_A(1) = rhob
346       rho_A(2) = rhob
347       rho_A(3) = 0.d0          ! beta  equals zero
348       delrho_A(1,1) = delrho_t(1,2) ! nabla n_up x
349       delrho_A(2,1) = delrho_t(2,2) ! nabla n_up y
350       delrho_A(3,1) = delrho_t(3,2) ! nabla n_up z
351       delrho_A(1,2) = 0.d0     ! set beta gradient to zero
352       delrho_A(2,2) = 0.d0     ! set beta gradient to zero
353       delrho_A(3,2) = 0.d0     ! set beta gradient to zero
354
355       neFSP = 0.0d0            !Ec in PBE
356       dneFSPdn(1) = 0.0d0      !Amat in PBE
357       dneFSPdn(2) = 0.0d0      !Amat in PBE
358       dneFSPdg(1) = 0.0d0      !Cmat in PBE
359       dneFSPdg(2) = 0.0d0      !Cmat in PBE
360       dneFSPdg(3) = 0.0d0      !Cmat in PBE
361c
362       call xc_cMpbe96(tol_rho, rho_A, delrho_A,
363     &        dneFSPdn,dneFSPdg, 1, 2, neFSP)
364       pbedown = neFSP
365c
366c      functional deriv info below fffffffffffff
367       etildeb= pbedown
368       detilbdna=0d0
369       detilbdnb = dneFSPdn(1)
370       detilbdgaa=0d0
371       detilbdgbb = dneFSPdg(D1_GAA)
372c
373c     n_sigma/n_total factor
374       fabdown=rhob/rhoval
375       delcb = -(1.d0 + ccc)*fabdown*xx2b*etildeb
376       erevc = erevc + delcb
377c
378       if(ldew) func(n) = func(n) + rhoval*erevc*fac
379       Ec = Ec + rhoval*erevc*qwght(n)*fac
380c
381c na
382       dxx1dna = -0.125d0*gaa/(taun*rhoa2)
383       dxx1adna = -0.125d0*gaa/(tauna*rhoa2)
384       atermn=pbe*ccc*2.d0*xx1*dxx1dna + (1.d0+ccc*xx2)*dneggadn(1)
385       btermn= (1.d0+ccc)*(2.d0*xx1a*dxx1adna*fabup*etildea +
386     &                     xx2a*etildea*fabdown/rhoval +
387     &                     xx2a*fabup*detiladna -
388     &                     xx2b*etildeb*fabdown/rhoval)
389       drevdna = atermn - btermn
390c
391c nb
392       dxx1dnb = -0.125d0*gbb/(taun*rhob2)
393       dxx1bdnb = -0.125d0*gbb/(taunb*rhob2)
394       atermn=pbe*ccc*2.d0*xx1*dxx1dnb + (1.d0+ccc*xx2)*dneggadn(2)
395       btermn= (1.d0+ccc)*(2.d0*xx1b*dxx1bdnb*fabdown*etildeb +
396     &                     xx2b*etildeb*fabup/rhoval +
397     &                     xx2b*fabdown*detilbdnb -
398     &                     xx2a*etildea*fabup/rhoval)
399       drevdnb = atermn - btermn
400c
401c gaa
402       dxx1dgaa = 0.125d0/(taun*rhoa)
403       dxx1adgaa = 0.125d0/(tauna*rhoa)
404       atermg=(1.d0+ccc*xx2)*dneggadg(D1_GAA)+ pbe*ccc*2.d0*xx1*dxx1dgaa
405       btermg=(1.d0+ccc)*(2.d0*xx1a*dxx1adgaa*fabup*etildea +
406     &    xx2a*fabup*detiladgaa)
407       drevdgaa = atermg - btermg
408c
409c gbb
410       dxx1dgbb = 0.125d0/(taun*rhob)
411       dxx1bdgbb = 0.125d0/(taunb*rhob)
412       atermg=(1.d0+ccc*xx2)*dneggadg(D1_GBB)+ pbe*ccc*2.d0*xx1*dxx1dgbb
413       btermg=(1.d0+ccc)*(2.d0*xx1b*dxx1bdgbb*fabdown*etildeb +
414     &   xx2b*fabdown*detilbdgbb)
415       drevdgbb = atermg - btermg
416c
417c gab
418       atermg=(1.d0+ccc*xx2)*dneggadg(D1_GAB)
419       drevdgab = atermg
420c
421c ta
422       dxx1dta=-xx1/taun
423       dxx1adta=-xx1a/tauna
424       drevdta=pbe*2.d0*ccc*xx1*dxx1dta
425     &        -(1.d0+ccc)*2.d0*xx1a*dxx1adta*fabup*etildea
426c
427c tb
428       dxx1dtb=-xx1/taun
429       dxx1bdtb=-xx1b/taunb
430       drevdtb=pbe*2.d0*ccc*xx1*dxx1dtb
431     &        -(1.d0+ccc)*2.d0*xx1b*dxx1bdtb*fabdown*etildeb
432c
433c derivs wrt na,nb
434       finalna= rhoval*drevdna + erevc
435       Amat(n,1)=Amat(n,1)+finalna*fac
436
437       finalnb= rhoval*drevdnb + erevc
438       Amat(n,2)=Amat(n,2)+finalnb*fac
439c
440c     derivs wrt gaa
441       finalgaa=rhoval*drevdgaa
442       Cmat(n,D1_GAA)=Cmat(n,D1_GAA)+ finalgaa*fac
443c
444c     derivs wrt gbb
445       finalgbb=rhoval*drevdgbb
446       Cmat(n,D1_GBB)=Cmat(n,D1_GBB)+ finalgbb*fac
447c
448c     derivs wrt gab
449       finalgab=rhoval*drevdgab
450       Cmat(n,D1_GAB)=Cmat(n,D1_GAB)+ finalgab*fac
451c
452c     derivs wrt ta,tb
453       apartt=rhoval*drevdta
454       finalt=apartt
455       Mmat(n,1)=Mmat(n,1)+0.5d0*finalt*fac
456
457       apartt=rhoval*drevdtb
458       finalt=apartt
459       Mmat(n,2)=Mmat(n,2)+0.5d0*finalt*fac
460
46120     continue
462
463      endif
464
465      return
466      end
467
468c
469      Subroutine xc_cpkzb99_d2()
470      call errquit(' cpkzb99: d2 not coded ',0,0)
471      return
472      end
473