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