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