1c Perdew-Kurth-Zupan-Blaha '99 exchange functional 2c META GGA 3C utilizes ingredients: 4c rho - density 5c delrho - gradient of density 6c tau - K.S kinetic energy density 7c cor only: tauW - von Weiszacker kinetic energy density 8c References: 9c [a] J.P. Perdew, S. Kurth, A. Zupan and P. Blaha, 10c PRL 82, 2544 (1999). 11 12 13 Subroutine xc_xpkzb99(tol_rho, fac, lfac, nlfac, rho, delrho, 14 & Amat, Cmat, nq, ipol, Ex, 15 & qwght, ldew, func, tau,Mmat) 16 17 18c 19c$Id$ 20c 21 implicit none 22c 23c 24 double precision fac, Ex 25 integer nq, ipol 26 logical lfac, nlfac,ldew 27 double precision func(*) ! value of the functional [output] 28c 29c Charge Density & Its Cube Root 30c 31 double precision rho(nq,ipol*(ipol+1)/2) 32c 33c Charge Density Gradient 34c 35 double precision delrho(nq,3,ipol) 36c 37c Quadrature Weights 38c 39 double precision qwght(nq) 40c 41c Sampling Matrices for the XC Potential & Energy 42c 43 double precision Amat(nq,ipol), Cmat(nq,*) 44 double precision tol_rho, pi 45c 46 integer n 47 double precision rrho, rho43, rho13, gamma 48c 49c kinetic energy density or tau 50c 51 double precision tau(nq,ipol), Mmat(nq,*) 52 double precision tauN 53 54 double precision p, qtil, x 55 double precision rho53, rho83, mt 56 double precision F83, F23, F53, F43, F13 57 double precision G920 58 double precision bigD, uk, ruk 59 double precision CX1, CX2, CX3, CX4 60 double precision P32, Ax 61c functional derivatives below 62 double precision d1p(3), d1x(3) 63 double precision rg2, d1mt(3) 64 double precision d1qtil(3), T1 65c functional derivatives above 66 67 parameter(uk=0.8040d0, bigD=0.113d0, ruk=1.d0/uk) 68 parameter (F43=4.d0/3.d0, F13=1.d0/3.d0) 69 parameter (F83=8.d0/3.d0, F23=2.d0/3.d0, F53=5.d0/3.d0) 70 parameter (G920 =9.d0/20.d0 ) 71 parameter (CX1 = 10.d0/81.d0, 72 & CX2 = 146.d0/2025.d0, 73 & CX3 = -73.d0/405.d0 ) 74 parameter (CX4= (bigD + ruk*CX1**2) ) 75c 76 pi=acos(-1d0) 77 Ax = (-0.75d0)*(3d0/pi)**F13 78 P32 = 1.d0/( 2.d0*(3.d0*pi**2)**(F23) ) 79 80c 81 if (ipol.eq.1 )then 82c 83c ======> SPIN-RESTRICTED <====== 84c 85 86 do 10 n = 1, nq 87 if (rho(n,1).lt.tol_rho) goto 10 88 89c rho43= n*e_x^unif=exchange energy per electron for uniform electron gas 90c = n* Ax*n^(1/3) or n*C*n^(1/3) 91 92 rho43 = Ax*rho(n,1)**F43 ! Ax*n^4/3 93 rrho = 1d0/rho(n,1) ! reciprocal of rho 94 rho13 = F43*rho43*rrho !functional deriv of rho43 95 96 rho53 = rho(n,1)**F53 97 rho83 = rho(n,1)**F83 98 99 100C Below we just sum up the LDA contribution to the functional 101 if (lfac) then 102 Ex = Ex + rho43*qwght(n)*fac 103 Amat(n,1) = Amat(n,1) + rho13*fac 104 if(ldew)func(n) = func(n) + rho43*fac 105 endif 106c 107 gamma = delrho(n,1,1)*delrho(n,1,1) + 108 & delrho(n,2,1)*delrho(n,2,1) + 109 & delrho(n,3,1)*delrho(n,3,1) 110c gam12 = dsqrt(gamma) 111c if (.not.(nlfac.and.gam12.gt.tol_rho)) goto 10 112 tauN = tau(n,1) 113 114 p=P32*gamma/(rho83*2.d0) 115 qtil=(3.d0*tauN*P32/rho53)-G920-(p/12.d0) 116c 117c Evaluate the GC part of Fx, i.e. mt = Fx(p,qtil) - 1 118c 119 120 x= CX1*p + CX2*qtil*qtil + CX3*qtil*p+ CX4*p*p 121 122 if (.not.(nlfac.and.x.gt.tol_rho)) goto 10 123 124 125 mt = uk - uk/(1.d0 + x*ruk) 126 127 128C functional derivatives 129 130 rg2=1.d0/( (1.d0 + x*ruk)*(1.d0 + x*ruk) ) 131 132c deriv wrt n, the density (for Amat) 133 d1p(1) = -F83*rrho*p 134 135 T1=3.d0*P32*tauN*(1.d0/rho53) 136 d1qtil(1) = -F53*T1*rrho - d1p(1)/12.d0 137 138 d1x(1) = CX1*d1p(1) + CX2*2.d0*qtil*d1qtil(1) + 139 & CX3*(qtil*d1p(1) + p*d1qtil(1)) + CX4*2d0*p*d1p(1) 140 141 d1mt(1) = rg2*d1x(1) 142 143 144c deriv wrt gamma, the gradient (for Cmat) 145 d1p(2) = 0.5d0*P32/rho83 146 d1qtil(2) = -d1p(2)/12.d0 147 148 d1x(2) = CX1*d1p(2) + CX2*2d0*qtil*d1qtil(2) + 149 & CX3*(qtil*d1p(2) + p*d1qtil(2)) + CX4*2d0*p*d1p(2) 150 151 d1mt(2) = rg2*d1x(2) 152 153 154c deriv wrt tau, the Kinetic Energy Density (for Mmat) 155 156 d1p(3) = 0.d0 157c d1qtil(3) = -d1p(2)/12.d0 am sure this is wrong 158 159 d1qtil(3) = 3.d0*P32/rho53 160 161 d1x(3) = CX2*2.d0*qtil*d1qtil(3) + 162 & CX3*p*d1qtil(3) 163 164 d1mt(3) = rg2*d1x(3) 165 166 167 168C Below we add the MetaGGA correction to the LDA part from above 169 170 if(ldew)func(n) = func(n) + rho43*mt*fac 171 Ex = Ex + rho43*mt*qwght(n)*fac 172 173 Amat(n,1) =Amat(n,1) + (mt*rho13 + rho43*d1mt(1))*fac 174 175 Cmat(n,1) = Cmat(n,1) + 2.d0*(rho43*d1mt(2)*fac) 176c check on this two or one 177 178 Mmat(n,1) = Mmat(n,1) +0.5d0*(rho43*d1mt(3)*fac) 179 180 181 10 continue 182 183 else 184c 185c ======> SPIN-UNRESTRICTED <====== 186 187c 188c use spin density functional theory ie n-->2n 189c Ex=(1/2)Ex[2*alpha] + (1/2)Ex[2*beta] 190 191 do 20 n = 1, nq 192 if (rho(n,1).lt.tol_rho) goto 20 193c 194c Alpha ALPHA ALPHA 195c 196 if (rho(n,2).lt.tol_rho) goto 25 197 rho43 = Ax*(2d0*rho(n,2))**F43 ! spin scaled 198 rrho = 0.5d0/rho(n,2) ! spin scaled 199 rho13 = F43*rho43*rrho !spin scaled & (1/2)factor 200 201 rho53 = (2.d0*rho(n,2))**F53 ! spin scaled 202 rho83 = (2.d0*rho(n,2))**F83 ! spin scaled 203c note all the "rho" quantities ARE spin scaled, for later use 204 205 if (lfac) then 206 Ex = Ex + rho43*qwght(n)*fac*0.5d0 207 Amat(n,1) = Amat(n,1) + rho13*fac 208 if(ldew)func(n) = func(n) + rho43*fac*0.5d0 209 endif 210c 211 gamma = delrho(n,1,1)*delrho(n,1,1) + 212 & delrho(n,2,1)*delrho(n,2,1) + 213 & delrho(n,3,1)*delrho(n,3,1) 214c NOTE gamma above is not spin scaled. that is why 215c -there is 4.d0*gamma in p 216c -there is 2.0 in the gam12 term 217 218 219 220c gam12 = 2d0*dsqrt(gamma) 221c if (.not.(nlfac.and.gam12.gt.tol_rho)) goto 25 222 223c below note factor of two for spin scaling 224 tauN = 2.d0* tau(n,1) 225 226c 227c Evaluate the GC part of Fx, i.e. mt(x) = Fx - 1 228c 229 230 p=0.5d0*P32*(4.d0*gamma)/rho83 231 qtil=(3.d0*tauN*P32/rho53) - G920 - (p/12.d0) 232 233 234 x= CX1*p + CX2*qtil*qtil + CX3*qtil*p+ CX4*p*p 235 if (.not.(nlfac.and.x.gt.tol_rho)) goto 25 236 237 rg2=1.d0/( (1.d0 + x*ruk)*(1.d0 + x*ruk) ) 238 239c ccccccc deriv wrt n, the density 240 241 d1p(1) = p*(-F83)*(2.d0*rrho) ! spin scaled 242 243 T1=3.d0*P32*tauN/rho53 ! spin scaled 244 d1qtil(1) = -F53*T1*2.d0*rrho - d1p(1)/12.d0 245 246 d1x(1) = CX1*d1p(1) + CX2*2.d0*qtil*d1qtil(1) + 247 & CX3*(qtil*d1p(1) + p*d1qtil(1)) + CX4*2d0*p*d1p(1) 248 249 d1mt(1) = rg2*d1x(1) 250 251 252c deriv wrt gamma, the gradient 253 d1p(2) = 0.5d0*P32*4.d0/rho83 ! spin scaled 254 d1qtil(2) = -d1p(2)/12.d0 !spin scaled 255 256 d1x(2) = CX1*d1p(2) + CX2*2d0*qtil*d1qtil(2) + 257 & CX3*(qtil*d1p(2) + p*d1qtil(2)) + CX4*2d0*p*d1p(2) 258 259 d1mt(2) = rg2*d1x(2) 260 261 262c deriv wrt tau, the Kinetic Energy Density 263 264c d1p(3) = 0.d0 term shown for completeness 265 d1qtil(3) = 3.d0*P32/rho53 266 267 d1x(3) = CX2*2.d0*qtil*d1qtil(3) + 268 & CX3*p*d1qtil(3) 269 270 d1mt(3) = rg2*d1x(3) 271 272 mt = uk - uk/(1.d0 + x*ruk) 273 274 Ex = Ex + rho43*mt*qwght(n)*fac*0.5d0 275 if(ldew)func(n) = func(n) + rho43*mt*fac*0.5d0 276 277 Amat(n,1)=Amat(n,1)+(mt*rho13 + 0.5d0*rho43*d1mt(1))*fac 278c note that the (.5) is built into the rho13 term already 279c hence we only need to put it onto the second term in Amat 280 281 282 Cmat(n,1) = Cmat(n,1) + (0.5d0*rho43*d1mt(2)*fac) 283 Mmat(n,1) = Mmat(n,1) +1.0d0*(0.5d0*rho43*d1mt(3)*fac) 284 285 286c 287c Beta BETA BETA 288c 289 25 continue 290 291 if (rho(n,3).lt.tol_rho) goto 20 292 rho43 = Ax*(2d0*rho(n,3))**F43 ! Ax (2 nBeta)^4/3 293 rrho = 0.5d0/rho(n,3) ! 1/(2 nBeta) 294 rho13 = F43*rho43*rrho !spin scaled func deriv of rho43 295 296 rho53 = (2.d0*rho(n,3))**F53 297 rho83 = (2.d0*rho(n,3))**F83 298 299C note all "rho" quantities above are spin scaled for later use 300 301 if (lfac) then 302 Ex = Ex + rho43*qwght(n)*fac*0.5d0 303 Amat(n,2) = Amat(n,2) + rho13*fac 304 if(ldew)func(n) = func(n) + rho43*fac*0.5d0 305 endif 306c 307 308 gamma = delrho(n,1,2)*delrho(n,1,2) + 309 & delrho(n,2,2)*delrho(n,2,2) + 310 & delrho(n,3,2)*delrho(n,3,2) 311c NOTE gamma above is not spin scaled. that is why 312c -there is 4.d0*gamma in p term 313c -there is 2.0 in the gam12 term 314 315 316c gam12 = 2d0*dsqrt(gamma) 317c if (.not.(nlfac.and.gam12.gt.tol_rho)) goto 20 318 319c below note factor of two for spin scaling 320 tauN = 2.d0* tau(n,2) 321 322 323c 324c Evaluate the GC part of F(x), i.e. mt(x) = Fx - 1 325c 326 327 328 p=0.5d0*P32*(4.d0*gamma)/rho83 329 qtil=(3.d0*tauN*P32/rho53) - G920 - (p/12.d0) 330 331 332 x= CX1*p + CX2*qtil*qtil + CX3*qtil*p+ CX4*p*p 333 334 335 if (.not.(nlfac.and.x.gt.tol_rho)) goto 20 336 337 338 rg2=1.d0/( (1.d0 + x*ruk)*(1.d0 + x*ruk) ) 339 340c ccccccc deriv wrt n, the density 341 342 d1p(1) = p*(-F83)*(2.d0*rrho) ! spin scaled 343 344 T1=3.d0*P32*tauN/rho53 ! spin scaled 345 d1qtil(1) = -F53*T1*2.d0*rrho - d1p(1)/12.d0 346 347 d1x(1) = CX1*d1p(1) + CX2*2.d0*qtil*d1qtil(1) + 348 & CX3*(qtil*d1p(1) + p*d1qtil(1)) + CX4*2d0*p*d1p(1) 349 350 d1mt(1) = rg2*d1x(1) 351 352 353c deriv wrt gamma, the gradient 354 d1p(2) = 0.5d0*P32*4.d0/rho83 ! spin scaled 355 d1qtil(2) = -d1p(2)/12.d0 !spin scaled 356 357 d1x(2) = CX1*d1p(2) + CX2*2d0*qtil*d1qtil(2) + 358 & CX3*(qtil*d1p(2) + p*d1qtil(2)) + CX4*2d0*p*d1p(2) 359 360 d1mt(2) = rg2*d1x(2) 361 362 363 364c deriv wrt tau, the Kinetic Energy Density 365 366C d1p(3) = 0.d0 included for completeness 367 d1qtil(3) = 3.d0*P32/rho53 368 369 370 d1x(3) = CX2*2.d0*qtil*d1qtil(3) + 371 & CX3*p*d1qtil(3) 372 373 d1mt(3) = rg2*d1x(3) 374 375 mt = uk - uk/(1.d0 + x*ruk) 376 377 Ex = Ex + rho43*mt*qwght(n)*fac*0.5d0 378 if(ldew)func(n) = func(n) + rho43*mt*fac*0.5d0 379 380 Amat(n,2)=Amat(n,2)+(mt*rho13 + 0.5d0*rho43*d1mt(1))*fac 381c note that the (.5) is built into the rho13 term already 382c hence we only need to put it onto the second term in Amat 383 384 Cmat(n,3) = Cmat(n,3) + 0.5d0*rho43*d1mt(2)*fac 385 Mmat(n,2) = Mmat(n,2) + 0.5d0*rho43*d1mt(3)*fac 386 38720 continue 388 endif 389c 390 return 391 end 392 393 394 395 396 Subroutine xc_xpkzb99_d2() 397 call errquit(' xpkzb99: d2 not coded ',0,0) 398 return 399 end 400 401 402