1c note that cfac, lcfac, nlcfac are single numbers
2c in the original cpbe96 file, they are arrays
3
4      Subroutine xc_cMpbe96(tol_rho, rho, delrho,
5     &                    Amat, Cmat, nq, ipol, Ec)
6
7
8c
9c$Id$
10c
11      Implicit none
12#include "dft2drv.fh"
13c
14c
15c     Input and other parameters
16c
17      integer lnq
18      integer ipol, nq
19      double precision dummy(1)
20      logical lfac, nlfac, lfacl, nlfacl
21c note that cfac, lcfac, nlcfac are single numbers
22c in the original cpbe96 file, they are arrays
23      double precision fac, facl
24      double precision lqwght
25      double precision tol_rho
26c
27c     Constants in PBE functional
28c
29      double precision GAMMA, BETA, PI
30      parameter (GAMMA = 0.03109069086965489503494086371273d0)
31      parameter (BETA = 0.06672455060314922d0)
32      parameter (PI = 3.1415926535897932385d0)
33c
34c     Threshold parameters
35c
36      double precision TOLL, EXPTOL
37      double precision EPS
38      parameter (TOLL = 1.0D-40, EXPTOL = 40.0d0)
39      parameter (EPS = 1.0e-8)
40      double precision cfac
41      parameter(cfac=1d0)
42c
43c     Correlation energy
44c
45      double precision Ec
46c
47c     Charge Density
48c
49      double precision rho(nq,ipol*(ipol+1)/2)
50      double precision rho_t(3)
51c
52c     Charge Density Gradient
53c
54      double precision delrho(nq,3,ipol)
55      double precision dsqgamma
56c
57c     Quadrature Weights
58c
59c     Sampling Matrices for the XC Potential
60c
61      double precision Amat(nq,ipol), Cmat(nq,*)
62c
63c     Intermediate derivative results, etc.
64c
65      integer n
66      double precision rhoval, gammaval
67      double precision nepsc, dnepscdn(2)
68      double precision epsc, depscdna, depscdnb
69      double precision H0, dH0dna, dH0dnb, dH0dg
70      double precision phi, dphidna, dphidnb, dphidzeta
71      double precision zeta, dzetadna, dzetadnb
72      double precision arglog, darglogdna, darglogdnb, darglogdg
73      double precision fAt, dfAtdt, dfAtdA
74      double precision fAtnum, dfAtnumdt, dfAtnumdA
75      double precision fAtden, dfAtdendt, dfAtdendA
76      double precision dfAtdna, dfAtdnb, dfAtdg
77      double precision A, dAdna, dAdnb
78      double precision t, dtdna, dtdnb, dtdg
79      double precision ks, dksdna, dksdnb
80      double precision argexp, dargexpdna, dargexpdnb
81      double precision expinA
82c
83c References:
84c [a] J. P. Perdew, K. Burke, and M. Ernzerhof,
85c     {\it Generalized gradient approximation made simple},
86c     Phys.\ Rev.\ Lett. {\bf 77,} 3865 (1996).
87c [b] J. P. Perdew, K. Burke, and Y. Wang, {\it Real-space cutoff
88c     construction of a generalized gradient approximation: The PW91
89c     density functional}, submitted to Phys.\ Rev.\ B, Feb. 1996.
90c [c] J. P. Perdew and Y. Wang, Phys.\ Rev.\ B {\bf 45}, 13244 (1992).
91c
92c  E_c(PBE) = Int n (epsilon_c + H0) dxdydz
93c
94c  n*epsilon_c                <=== supplied by another subroutine
95c  d(n*epsilon_c)/d(na)       <=== supplied by another subroutine
96c  d2(n*epsilon_c)/d(na)d(na) <=== supplied by another subroutine
97c  d2(n*epsilon_c)/d(na)d(nb) <=== supplied by another subroutine
98c  d2(n*epsilon_c)/d(nb)d(nb) <=== supplied by another subroutine
99c
100c  H0 = GAMMA * phi**3 * log{ 1 + BETA/GAMMA * t**2 * [ ... ]}
101c
102c  phi = (1/2)[(1+zeta)**(2/3)+(1-zeta)**(2/3)]
103c
104c  zeta = (na - nb)/n
105c
106c  [ ... ] = (1 + A * t**2)/(1 + A * t**2 + A**2 * t**4)
107c
108c  A = BETA/GAMMA [exp{-epsilon_c/(GAMMA*phi**3)}-1]**(-1)
109c
110c  t = |Nabla n|/(2*phi*ks*n)
111c
112c  ks = 2 * (3 * PI**2 * n)**(1/6) / sqrt(PI)
113c
114c  |Nabla n| = sqrt(g_aa + g_bb + 2*g_ab)
115c
116c  Names of variables
117c
118c  E_c(PBE)                  : Ec
119c  n (alpha+beta density)    : rhoval
120c  na, nb                    : rho(*,2), rho(*,3)
121c  epsilon_c                 : epsc
122c  H0                        : H0
123c  n*epsilon_c               : nepsc
124c  phi                       : phi
125c  zeta                      : zeta
126c  { ... }                   : arglog
127c  [ ... ]                   : fAt
128c  (1 + A * t**2)            : fAtnum
129c  (1 + A * t**2 + A**2 * t**4) : fAtden
130c  A                         : A
131c  t                         : t
132c  |Nabla n|                 : gammaval
133c  ks                        : ks
134c  {-epsilon_c ... }         : argexp
135c  g_aa, g_bb, g_ab          : g
136c
137c  Derivatives of these are named like d...dna, d2...dnadnb,
138c  d2...dna2, etc.
139c
140
141c      write(0,*) 'upon arrival in cpbe  Ec=',Ec
142
143      fac = cfac
144      lfac = .false.
145      nlfac = .true.
146
147c
148c     ======> BOTH SPIN-RESTRICETED AND UNRESTRICTED <======
149c
150      do 20 n = 1, nq
151c
152c        n and zeta = (na - nb)/n
153c
154         rhoval = rho(n,1)
155         rho_t(1) = rho(n,1)
156         if (ipol.eq.2) then
157            rho_t(2) = rho(n,2)
158            rho_t(3) = rho(n,3)
159         endif
160         if (rhoval.le.tol_rho) goto 20
161         if (ipol.eq.1) then
162            gammaval = delrho(n,1,1)*delrho(n,1,1) +
163     &                 delrho(n,2,1)*delrho(n,2,1) +
164     &                 delrho(n,3,1)*delrho(n,3,1)
165         else
166            gammaval = delrho(n,1,1)*delrho(n,1,1) +
167     &                 delrho(n,1,2)*delrho(n,1,2) +
168     &                 delrho(n,2,1)*delrho(n,2,1) +
169     &                 delrho(n,2,2)*delrho(n,2,2) +
170     &                 delrho(n,3,1)*delrho(n,3,1) +
171     &                 delrho(n,3,2)*delrho(n,3,2) +
172     &           2.d0*(delrho(n,1,1)*delrho(n,1,2) +
173     &                 delrho(n,2,1)*delrho(n,2,2) +
174     &                 delrho(n,3,1)*delrho(n,3,2))
175         endif
176         dsqgamma = max(dsqrt(gammaval),tol_rho)
177         nepsc = 0.0d0
178         dnepscdn(1) = 0.0d0
179         if (ipol.eq.2) dnepscdn(2) = 0.0d0
180c
181c        call for LDA bit
182c        this implementation temporarily assigns the pw91LDA for
183c        use in the metaGGA local part
184c
185
186            call xc_pw91lda(tol_rho,1d0,.true.,.false.,rho_t,
187     &         dnepscdn,1,ipol,nepsc,1d0,
188     &         .false.,dummy)
189
190c
191c        ==================
192c        PBE non-local part
193c        ==================
194         if(abs(nepsc).lt.tol_rho*tol_rho) goto 20
195c
196c        epsilon_c = n*epsilon_c / n
197c
198         epsc = nepsc/rhoval
199         if (ipol.eq.1) then
200            depscdna = dnepscdn(1)/rhoval-nepsc/(rhoval**2)
201            depscdnb = depscdna
202         else
203            depscdna = dnepscdn(1)/rhoval-nepsc/(rhoval**2)
204            depscdnb = dnepscdn(2)/rhoval-nepsc/(rhoval**2)
205         endif
206c
207c        ks = 2*(3*PI**2*n)**(1/6)/sqrt(PI) and its derivs
208c
209         ks = 2.0d0*(3.0d0*PI*PI*rhoval)**(1.0d0/6.0d0)/dsqrt(PI)
210         dksdna = (1.0d0/6.0d0)*ks/rhoval
211         dksdnb = dksdna
212c
213c        zeta = (na-nb)/n and its derivs
214c
215         if (ipol.eq.1) then
216            zeta = 0.0d0
217         else
218            zeta = (rho(n,2)-rho(n,3))/rhoval
219         endif
220         if(zeta.lt.-1.0d0) zeta=-1.0d0
221         if(zeta.gt. 1.0d0) zeta= 1.0d0
222         if (ipol.eq.1) then
223            dzetadna = 1.0d0/rhoval
224            dzetadnb = -1.0d0/rhoval
225         else
226            dzetadna =  2.0d0*rho(n,3)/(rhoval**2)
227            dzetadnb = -2.0d0*rho(n,2)/(rhoval**2)
228         endif
229c
230c        phi = (1/2)[(1+zeta)**(2/3)+(1-zeta)**(2/3)] and its derivs
231c
232         phi = 0.5d0*((1.0d0+zeta)**(2.0d0/3.0d0)
233     &               +(1.0d0-zeta)**(2.0d0/3.0d0))
234         if ((1.0d0-zeta).lt.tol_rho) then
235            dphidzeta = 0.5d0*(2.0d0/3.0d0)*(
236     &             (1.0d0+zeta)**(2.0d0/3.0d0)/(1.0d0+zeta))
237         else if ((1.0d0+zeta).lt.tol_rho) then
238            dphidzeta = 0.5d0*(2.0d0/3.0d0)*(
239     &            -(1.0d0-zeta)**(2.0d0/3.0d0)/(1.0d0-zeta))
240         else
241            dphidzeta = 0.5d0*(2.0d0/3.0d0)*(
242     &         (1.0d0+zeta)**(2.0d0/3.0d0)/(1.0d0+zeta)
243     &        -(1.0d0-zeta)**(2.0d0/3.0d0)/(1.0d0-zeta))
244         endif
245         dphidna = dphidzeta*dzetadna
246         dphidnb = dphidzeta*dzetadnb
247c
248c        t = |Nabla n|/(2*phi*ks*n) and its derivs
249c
250         t = dsqgamma/(2.0d0*phi*ks*rhoval)
251         dtdna = -t/rhoval-t/phi*dphidna-t/ks*dksdna
252         dtdnb = -t/rhoval-t/phi*dphidnb-t/ks*dksdnb
253c
254c        { ... } in A (see below) and its derivs
255c
256         argexp = -epsc/GAMMA/(phi**3)
257         dargexpdna = -depscdna/GAMMA/(phi**3)
258     &                +3.0d0*epsc/GAMMA/(phi**4)*dphidna
259         dargexpdnb = -depscdnb/GAMMA/(phi**3)
260     &                +3.0d0*epsc/GAMMA/(phi**4)*dphidnb
261c
262c        A = BETA/GAMMA [exp{-epsilon_c/(GAMMA*phi**3)}-1]**(-1)
263c
264         if (dabs(argexp).lt.EXPTOL) then
265            expinA=dexp(argexp)
266         else
267            expinA=0.0d0
268         endif
269         A = BETA/GAMMA/(expinA-1.0d0)
270         dAdna = -BETA/GAMMA*dargexpdna*expinA/(expinA-1.0d0)**2
271         dAdnb = -BETA/GAMMA*dargexpdnb*expinA/(expinA-1.0d0)**2
272c
273c        fAt = (1 + A * t**2)/(1 + A * t**2 + A**2 * t**4) and its derivs
274c
275         fAtnum = 1.0d0+A*t**2
276         fAtden = 1.0d0+A*t**2+A**2*t**4
277         fAt = fAtnum/fAtden
278         dfAtnumdt = 2.0d0*A*t
279         dfAtnumdA = t**2
280         dfAtdendt = 2.0d0*A*t+4.0d0*A**2*t**3
281         dfAtdendA = t**2+2.0d0*A*t**4
282         dfAtdt = (dfAtnumdt*fAtden-fAtnum*dfAtdendt)/(fAtden**2)
283         dfAtdA = (dfAtnumdA*fAtden-fAtnum*dfAtdendA)/(fAtden**2)
284         dfAtdna = dfAtdt * dtdna + dfAtdA * dAdna
285         dfAtdnb = dfAtdt * dtdnb + dfAtdA * dAdnb
286c
287c        arglog = 1 + BETA/GAMMA * t**2 * fAt and its derivs
288c
289         arglog = 1.0d0 + BETA/GAMMA*t**2*fAt
290         darglogdna = BETA/GAMMA*(2.0d0*t*dtdna*fAt
291     &                            +t*t*dfAtdna)
292         darglogdnb = BETA/GAMMA*(2.0d0*t*dtdnb*fAt
293     &                            +t*t*dfAtdnb)
294c
295c        H0 = GAMMA * phi**3 * log{arglog} and its derivs
296c
297         H0 = GAMMA*(phi**3)*dlog(arglog)
298         dH0dna = GAMMA*(3.0d0*(phi**2)*dphidna*dlog(arglog)
299     &                  +(phi**3)*darglogdna/arglog)
300         dH0dnb = GAMMA*(3.0d0*(phi**2)*dphidnb*dlog(arglog)
301     &                  +(phi**3)*darglogdnb/arglog)
302c
303c        Now we update Ec, Amat, and Amat2
304c
305
306c          NOTE:  this PBE does the LDA part of Ec in house
307            Ec = Ec+epsc*fac
308            Ec = Ec+H0*fac
309            Amat(n,1) = Amat(n,1) + depscdna*fac
310            if (ipol.eq.2) Amat(n,2) = Amat(n,2) + depscdnb*fac
311
312            Amat(n,1) = Amat(n,1) +  dH0dna*fac
313            if (ipol.eq.2) Amat(n,2) = Amat(n,2) + dH0dnb*fac
314c
315c        Now we go into gradient-correction parts
316c        Note that the functional depends on |Nabla n| through "t" only
317c
318         if (dsqgamma.gt.TOLL)then
319            dtdg = 0.25d0/(phi*ks*rhoval)/dsqgamma
320            dfAtdg = dfAtdt*dtdg
321            darglogdg = BETA/GAMMA*(2.0d0*t*dtdg*fAt+t*t*dfAtdg)
322            dH0dg = GAMMA*(phi**3)*darglogdg/arglog
323
324            Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + dH0dg*fac
325            Cmat(n,D1_GAB) = Cmat(n,D1_GAB) + dH0dg*fac*2.0d0
326            if (ipol.eq.2) Cmat(n,D1_GBB) = Cmat(n,D1_GBB) + dH0dg*fac
327         endif
328   20 continue
329c
330
331      return
332      end
333c
334