1c 2c$Id$ 3c 4#include "dft2drv.fh" 5c 6c Tao,Perdew, Staroverov, Scuseria exchange functional 7c META GGA 8C utilizes ingredients: 9c rho - density 10c delrho - gradient of density 11c tau (tauN)- K.S kinetic energy density 12c tauW - von Weiszacker kinetic energy density 13c tauU - uniform-gas KE density 14 15 16 Subroutine xc_ctpss03(tol_rho, cfac, lcfac, nlcfac, rho, delrho, 17 & nq, ipol, Ec, qwght, ldew, func, 18 & tau, Amat, Cmat, Mmat) 19c References: 20c [a] J. Tao, J.P. Perdew, V.N.Staroverov, G. Scuseria 21c PRL 91, 146401 (2003). 22c [b] J. Tao, J.P. Perdew, V.N.Staroverov, G. Scuseria 23c JCP 120, 6898 (2004). 24 25 Implicit none 26c 27c 28c Input and other parameters 29c 30 integer ipol, nq 31 32 double precision cfac 33 logical ldew,lcfac,nlcfac 34 double precision func(*) ! value of the functional [output] 35 36 double precision fac 37 double precision tol_rho 38c 39c Correlation energy 40c 41 double precision Ec 42c 43c Charge Density 44c 45 double precision rho(nq,ipol*(ipol+1)/2) 46 47c 48c Charge Density Gradient 49c 50 double precision delrho(nq,3,ipol), gammaval, gam12 51 52c 53c Kinetic Energy Density 54c 55 double precision tau(nq,ipol), tauN 56 57c 58c Quadrature Weights 59c 60 double precision qwght(nq) 61c 62c Sampling Matrices for the XC Potential 63c 64 double precision Amat(nq,ipol), Cmat(nq,*) 65 double precision Mmat(nq,*) 66 67 integer n 68 double precision rhoval 69 70c call to the cPBE subroutine 71 72 double precision neGGA, dneGGAdn(2), dneGGAdg(3) 73 double precision rho_t(3), delrho_t(3,2) 74 double precision tauNA,tauNB 75 double precision neFSP, dneFSPdn(2), dneFSPdg(3) 76 double precision delrho_A(3,2), rho_A(3) 77 78c TPSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSs 79 80 double precision THRD, F43, F73 81 double precision zeta, ccc 82 double precision tauw 83 double precision xx2,xx3 84 double precision pbe, en 85 double precision dd 86 double precision rhoa,rhob 87 double precision pbeup, delc, erevc, erevsic, 88 A delca,delcb,pbedown 89 90c derivsssssssssssssssssssssssss 91 92 double precision decggadn, dtwt2dn, dtwt3dn 93 double precision dpbeupdn, drevdn 94 double precision atermn, btermn 95 double precision finaln 96 double precision drevdg, dpbeupdg 97 double precision dtwt2dg, dtwt3dg, decggadg 98 double precision atermg, btermg 99 double precision apartg, cpartg(2), finalg 100 double precision finalgaa,finalgbb,finalgab 101 double precision drevdt, apartt, bpartt(2),finalt 102 double precision drevdta,drevdtb 103 104 double precision dzetadna, dzetadnb 105 double precision dcccdna, dcccdnb 106 double precision dcccdgaa,dcccdgbb,dcccdgab 107 double precision drevdna, drevdnb 108 109 double precision drevdgaa 110 double precision drevdgbb 111 double precision drevdgab 112 double precision dtwt3dt 113 double precision etildea,etildeb, 114 N detiladna,detiladnb,detilbdna,detilbdnb, 115 D detiladgaa,detiladgbb,detilbdgaa,detilbdgbb 116 double precision detiladgab,detilbdgab 117 double precision pi,fabup,fabdown 118 double precision czeta0,czeta1 119 double precision gaa,gbb,gab 120 double precision xi2,delzeta2,dencxi2zeta, 121 D ddez2dna,ddez2dnb,dxi2dna,dxi2dnb, 122 D ddencxi2dna,ddencxi2dnb, 123 D ddencxi2dxi2,ddencxi2dzeta 124 double precision dxi2dgaa 125 double precision dxi2dgbb 126 double precision dxi2dgab 127 double precision denccc,onemzeta,onepzeta 128 double precision rhoval2,rhoval3 129 double precision denxi2,denxi22 130 double precision term1,term2 131 double precision ddelzeta2dna,ddelzeta2dnb 132 double precision ddelzeta2dgaa,ddelzeta2dgbb,ddelzeta2dgab 133 double precision ddencccdna,ddencccdnb 134c 135 parameter (dd = 2.8d0) 136 parameter (THRD = 1d0/3d0) 137 parameter (F43 = 4d0/3d0) 138 parameter (F73 = 7d0/3d0) 139c 140 czeta0(zeta)= 141 & 0.53d0+0.87d0*zeta**2+0.5d0*zeta**4+2.26d0*zeta**6 142 czeta1(zeta)= 143 & 2.d0*0.87d0*zeta+4d0*0.5d0*zeta**3+6d0*2.26d0*zeta**5 144 denccc(zeta,xi2)= 145 & 1.d0+0.5d0*xi2*((1.d0+zeta)**(-F43)+(1.d0-zeta)**(-F43)) 146c 147 pi=acos(-1d0) 148 fac = cfac 149c 150 if (ipol.eq.1 )then 151c ======> SPIN-RESTRICTED <====== 152 153 do 22 n = 1, nq 154 if (rho(n,1).lt.tol_rho) goto 22 155 156 rhoval = rho(n,1) 157 158C set up values to call PBE subroutine 159 rho_t(1) = rho(n,1) 160c do delrho 161 delrho_t(1,1) = delrho(n,1,1) 162 delrho_t(2,1) = delrho(n,2,1) 163 delrho_t(3,1) = delrho(n,3,1) 164 gammaval = delrho(n,1,1)*delrho(n,1,1) + 165 & delrho(n,2,1)*delrho(n,2,1) + 166 & delrho(n,3,1)*delrho(n,3,1) 167 gam12=dsqrt(gammaval) 168c 169c get E_GGA[rho,gamma] 170c 171 neGGA = 0.0d0 !Ec in PBE 172 dneGGAdn(1) = 0.0d0 !Amat in PBE 173 dneGGAdg(1) = 0.0d0 !Cmat in PBE 174 dneGGAdg(2) = 0.0d0 !Cmat in PBE 175 176 call xc_cMpbe96(tol_rho, 177 & rho_t, delrho_t, 178 & dneGGAdn,dneGGAdg, 179 & 1, ipol, neGGA) 180 pbe = neGGA 181 182c 183c epGGA = n*(epsilon_c^GGA) / n =cor. energy per electron 184c epGGA= ec^LDA +H = pbe 185 186 tauN = tau(n,1) 187 188 ccc = 0.53d0 !since zeta=0 189 190 if(sqrt(gammaval).gt.tol_rho) then 191 tauw = 0.125d0*gammaval/rhoval 192 xx2 = (tauw/tauN)**2.d0 193 xx3 = (tauw/tauN)**3.d0 194 else 195 tauw=0d0 196 xx2=0d0 197 xx3=0d0 198 endif 199 en = pbe*(1.d0 + ccc*xx2) 200c 201c set up values to call PBE subroutine as 202c Fully SpinPolarized system 203c 204 205 rho_A(1) = (0.5d0)*rho(n,1) ! total equals (1/2)n_tot 206 rho_A(2) = (0.5d0)*rho(n,1) ! alpha equals (1/2)n_tot 207 rho_A(3) = 0.d0 ! beta equals zero 208 delrho_A(1,1) = (0.5d0)*delrho_t(1,1) ! nabla n_up x 209 delrho_A(2,1) = (0.5d0)*delrho_t(2,1) ! nabla n_up y 210 delrho_A(3,1) = (0.5d0)*delrho_t(3,1) ! nabla n_up z 211 212 delrho_A(1,2) = 0.d0 ! set beta gradient to zero 213 delrho_A(2,2) = 0.d0 ! set beta gradient to zero 214 delrho_A(3,2) = 0.d0 ! set beta gradient to zero 215 216 neFSP = 0.0d0 !Ec in PBE 217 dneFSPdn(1) = 0.0d0 !Amat in PBE 218 dneFSPdn(2) = 0.0d0 !Amat in PBE 219 220 dneFSPdg(1) = 0.0d0 !Cmat in PBE 221 dneFSPdg(2) = 0.0d0 !Cmat in PBE 222 dneFSPdg(3) = 0.0d0 !Cmat in PBE 223 224c 225c get E_GGA[rho_alpha,0,gamma_alpha,0] 226c 227 call xc_cMpbe96(tol_rho, rho_A, delrho_A, 228 & dneFSPdn,dneFSPdg, 229 & 1, 2, neFSP) 230 231 pbeup = neFSP 232 233c functional deriv info below fffffffffffff 234 235 dtwt2dn = -2.d0*xx2/rhoval 236 dtwt3dn = -3.d0*xx3/rhoval 237 decggadn= dneGGAdn(1) 238 if(sqrt(gammaval).gt.tol_rho) then 239 dtwt2dg = 2.d0*0.125d0*tauw/(rhoval*tauN**2) 240 else 241 dtwt2dg = 0d0 242 endif 243 if(abs(taun).gt.tol_rho) then 244 dtwt3dg = 3.d0*xx2*0.125d0/(rhoval*tauN) 245 else 246 dtwt3dg = 0d0 247 endif 248 decggadg= dneGGAdg(1) 249 250 if (pbeup.lt.pbe) then 251 delc= xx2*pbe 252C eps-tilda is eps_c 253c erev = egga*(1-xx2) 254C functional deriv info below fffffffffffffffff 255 256 drevdn= -pbe*dtwt2dn+decggadn*(1.d0 - xx2) 257 drevdg= -pbe*dtwt2dg+decggadg*(1.d0 - xx2) 258 drevdt= 2.d0*pbe*xx2/tauN 259 else 260 delc= xx2*pbeup 261 262C eps-tilda is eps^FSP 263C functional deriv info below fffffffffffffffff 264 265 dpbeupdn = 0.5d0*dneFSPdn(1) 266c above note the .5's. you are taking deriv wrt total density n 267c not deriv wrt n_up 268 dpbeupdg = 0.25d0*dneFSPdg(1) 269c note .25 above is because you want gamma=deln_tot*deln_tot 270c 271 atermn=pbe*ccc*dtwt2dn+(1.d0+ccc*xx2)*decggadn 272 btermn=(1.d0+ccc)*(xx2*dpbeupdn+pbeup*dtwt2dn) 273 drevdn=atermn-btermn 274c 275 atermg=pbe*ccc*dtwt2dg+(1.d0+ccc*xx2)*decggadg 276 btermg=(1.d0+ccc)*(xx2*dpbeupdg+pbeup*dtwt2dg) 277 drevdg=atermg-btermg 278c 279 if(abs(taun).gt.tol_rho) then 280 drevdt=(ccc*pbe-(1.d0+ccc)*pbeup)*xx2*(-2.d0/tauN) 281 else 282 drevdt=0d0 283 endif 284 endif 285 286 delc = -(1.d0 + ccc)*delc 287 erevc = en + delc 288 erevsic = erevc*(1.d0 + dd*erevc*xx3) 289 290 if(ldew) func(n) = func(n) + rhoval*erevsic*fac 291 Ec = Ec + rhoval*erevsic*qwght(n)*fac 292 293c derivs wrt n 294 finaln= rhoval*drevdn + erevc + 295 & dd*(erevc*erevc*xx3 + 296 + rhoval*(xx3*2.d0*erevc*drevdn + 297 + erevc*erevc*dtwt3dn)) 298 Amat(n,1)=Amat(n,1)+(finaln)*fac 299 300c derivs wrt g 301 apartg=rhoval*drevdg 302 cpartg(1)=erevc*erevc*dtwt3dg 303 cpartg(2)=xx3*2.d0*erevc*drevdg 304 finalg=apartg+rhoval*dd*( cpartg(1)+cpartg(2) ) 305 Cmat(n,D1_GAA)=Cmat(n,D1_GAA)+ 2d0*finalg*fac 306 307c derivs wrt t 308 apartt=rhoval*drevdt 309 if(abs(taun).gt.tol_rho) then 310 bpartt(1)=-erevc*erevc*xx3*3.d0/tauN 311 else 312 bpartt(1)=0d0 313 endif 314 bpartt(2)=xx3*2.d0*erevc*drevdt 315 finalt=apartt+dd*rhoval*(bpartt(1)+bpartt(2)) 316 Mmat(n,1)=Mmat(n,1)+0.5d0*finalt*fac 317 31822 continue 319 320 else !ipol=2 321c 322c ======> SPIN-UNRESTRICTED <====== 323c 324 do 20 n = 1, nq 325c 326 if (rho(n,1).lt.tol_rho) goto 20 327 328 rhoval = rho(n,1) 329 330 rho_t(1) = rho(n,1) 331 rho_t(2) = rho(n,2) 332 rho_t(3) = rho(n,3) 333 delrho_t(1,1) = delrho(n,1,1) 334 delrho_t(2,1) = delrho(n,2,1) 335 delrho_t(3,1) = delrho(n,3,1) 336 delrho_t(1,2) = delrho(n,1,2) 337 delrho_t(2,2) = delrho(n,2,2) 338 delrho_t(3,2) = delrho(n,3,2) 339c 340 neGGA = 0.0d0 !Ec in PBE 341 dneGGAdn(1) = 0.0d0 !Amat in PBE (n,1) 342 dneGGAdn(2) = 0.0d0 !Amat in PBE (n,2) 343 dneGGAdg(1) = 0.0d0 !Cmat in PBE--aa 344 dneGGAdg(2) = 0.0d0 !Cmat in PBE--ab 345 dneGGAdg(3) = 0.0d0 !Cmat in PBE--bb 346c 347c get E_GGA[rho,gamma] 348c 349 call xc_cMpbe96(tol_rho, 350 & rho_t, delrho_t, 351 & dneGGAdn,dneGGAdg, 352 & 1, ipol, neGGA) 353 pbe = neGGA 354c 355c epGGA = (epsilon_c^GGA) =cor. energy per electron 356c epGGA= ec^LDA +H = pbe 357c 358 gammaval = delrho(n,1,1)*delrho(n,1,1) + 359 & delrho(n,1,2)*delrho(n,1,2) + 360 & delrho(n,2,1)*delrho(n,2,1) + 361 & delrho(n,2,2)*delrho(n,2,2) + 362 & delrho(n,3,1)*delrho(n,3,1) + 363 & delrho(n,3,2)*delrho(n,3,2) + 364 & 2.d0*(delrho(n,1,1)*delrho(n,1,2) + 365 & delrho(n,2,1)*delrho(n,2,2) + 366 & delrho(n,3,1)*delrho(n,3,2)) 367 gam12=dsqrt(gammaval) 368 tauNa = tau(n,1) 369 tauNb = tau(n,2) 370 taun = tauna+taunb 371 rhoa=rho(n,2) 372 rhob=rho(n,3) 373c 374c Check for small densities (H atom case as well) 375c 376 if ((rhoa.lt.tol_rho).or. 377 & (rhob.lt.tol_rho)) goto 20 378c 379 dcccdna=0.d0 380 dcccdnb=0.d0 381 dcccdgaa=0.d0 382 dcccdgbb=0.d0 383 dcccdgab=0.d0 384 ccc=0d0 385c 386 if(rhoval.gt.tol_rho) then 387c 388 zeta=(rhoa-rhob)/rhoval 389c 390 onepzeta = 1.d0+zeta 391 onemzeta = 1.d0-zeta 392c 393 gaa = delrho(n,1,1)*delrho(n,1,1) + 394 & delrho(n,2,1)*delrho(n,2,1) + 395 & delrho(n,3,1)*delrho(n,3,1) 396 gbb = delrho(n,1,2)*delrho(n,1,2) + 397 & delrho(n,2,2)*delrho(n,2,2) + 398 & delrho(n,3,2)*delrho(n,3,2) 399 gab = delrho(n,1,1)*delrho(n,1,2) + 400 & delrho(n,2,1)*delrho(n,2,2) + 401 & delrho(n,3,1)*delrho(n,3,2) 402c 403 rhoval2 = rhoval*rhoval 404 rhoval3 = rhoval*rhoval*rhoval 405 denxi2 = 2.d0*(3.d0*pi*pi*rhoval)**(1.d0/3.d0) 406 denxi22 = denxi2*denxi2 407 delzeta2 = (gaa*onemzeta*onemzeta + gbb*onepzeta*onepzeta 408 & -2.d0*gab*onemzeta*onepzeta)/rhoval2 409c 410 dzetadna = onemzeta/rhoval 411 dzetadnb = -onepzeta/rhoval 412c 413 term1=-2.d0*gaa*onemzeta+2.d0*gbb*onepzeta+4.d0*gab*zeta 414 term1=(term1/rhoval2)*dzetadna 415 term2=-2.d0*delzeta2/rhoval 416 ddelzeta2dna = term1 + term2 417c 418 term1=-2.d0*gaa*onemzeta+2.d0*gbb*onepzeta+4.d0*gab*zeta 419 term1=(term1/rhoval2)*dzetadnb 420 term2=-2.d0*delzeta2/rhoval 421 ddelzeta2dnb = term1 + term2 422c 423 ddelzeta2dgaa = onemzeta*onemzeta/rhoval2 424 ddelzeta2dgbb = onepzeta*onepzeta/rhoval2 425 ddelzeta2dgab = -2.d0*onepzeta*onemzeta/rhoval2 426c 427 xi2=delzeta2/denxi22 428 dxi2dna=(ddelzeta2dna -(2.d0/3d0)*delzeta2/rhoval)/denxi22 429 dxi2dnb=(ddelzeta2dnb -(2.d0/3d0)*delzeta2/rhoval)/denxi22 430c 431 dxi2dgaa=onemzeta*onemzeta/rhoval2/denxi22 432 dxi2dgbb=onepzeta*onepzeta/rhoval2/denxi22 433 dxi2dgab=-2.d0*onepzeta*onemzeta/rhoval2/denxi22 434c 435 ccc=czeta0(zeta)/(denccc(zeta,xi2)**4) 436 ddencccdna=2.d0*(denccc(zeta,xi2)**3)* 437 & (dxi2dna*(onepzeta**(-F43) + onemzeta**(-F43)) 438 & + xi2*F43*(onemzeta**(-F73) - onepzeta**(-F73))*dzetadna) 439 ddencccdnb=2.d0*(denccc(zeta,xi2)**3)* 440 & (dxi2dnb*(onepzeta**(-F43) + onemzeta**(-F43)) 441 & + xi2*F43*(onemzeta**(-F73) - onepzeta**(-F73))*dzetadnb) 442c 443 dcccdna=(czeta1(zeta)*dzetadna/denccc(zeta,xi2)**4) - 444 & (czeta0(zeta)*ddencccdna/(denccc(zeta,xi2)**8)) 445 dcccdnb=(czeta1(zeta)*dzetadnb/denccc(zeta,xi2)**4) - 446 & (czeta0(zeta)*ddencccdnb/(denccc(zeta,xi2)**8)) 447c 448 dcccdgaa=-(czeta0(zeta)/(denccc(zeta,xi2)**8))* 449 & 2.d0*(denccc(zeta,xi2)**3)* 450 & dxi2dgaa*(onepzeta**(-F43) + onemzeta**(-F43)) 451 dcccdgbb=-(czeta0(zeta)/(denccc(zeta,xi2)**8))* 452 & 2.d0*(denccc(zeta,xi2)**3)* 453 & dxi2dgbb*(onepzeta**(-F43) + onemzeta**(-F43)) 454 dcccdgab=-(czeta0(zeta)/(denccc(zeta,xi2)**8))* 455 & 2.d0*(denccc(zeta,xi2)**3)* 456 & dxi2dgab*(onepzeta**(-F43) + onemzeta**(-F43)) 457 endif 458c 459 tauw = 0.125d0*gammaval/rhoval 460 if(abs(tauw).gt.tol_rho) then 461c 462 xx2 = (tauw/tauN)**2.d0 463 xx3 = (tauw/tauN)**3.d0 464 dtwt2dn = -2.d0*xx2/rhoval 465 dtwt3dn = -3.d0*xx3/rhoval 466 dtwt3dt = -3.d0*xx3/taun 467 dtwt2dg = 2.d0*0.125d0*tauw/(rhoval*tauN**2) 468 dtwt3dg = 3.d0*xx2*0.125d0/(rhoval*tauN) 469 else 470 tauw=0d0 471 xx2=0d0 472 xx3=0d0 473 dtwt2dn = 0d0 474 dtwt3dn = 0d0 475 dtwt3dt = 0d0 476 dtwt2dg = 0d0 477 dtwt3dg = 0d0 478 479 endif 480c 481 en = pbe*(1.d0 + ccc*xx2) 482c 483c Alpha bit 484c set up values to call PBE subroutine as 485c Fully SpinPolarized system for alpha spin 486c to get E_GGA[rho_alpha,0,gamma_alpha,0] 487c 488 rho_A(1) = rhoa 489 rho_A(2) = rhoa 490 rho_A(3) = 0.d0 ! beta equals zero 491 delrho_A(1,1) = delrho_t(1,1) ! nabla n_up x 492 delrho_A(2,1) = delrho_t(2,1) ! nabla n_up y 493 delrho_A(3,1) = delrho_t(3,1) ! nabla n_up z 494 delrho_A(1,2) = 0.d0 ! set beta gradient to zero 495 delrho_A(2,2) = 0.d0 ! set beta gradient to zero 496 delrho_A(3,2) = 0.d0 ! set beta gradient to zero 497 498 neFSP = 0.0d0 !Ec in PBE 499 dneFSPdn(1) = 0.0d0 !Amat in PBE 500 dneFSPdn(2) = 0.0d0 !Amat in PBE 501 502 dneFSPdg(1) = 0.0d0 !Cmat in PBE 503 dneFSPdg(2) = 0.0d0 !Cmat in PBE 504 dneFSPdg(3) = 0.0d0 !Cmat in PBE 505c 506 call xc_cMpbe96(tol_rho, rho_A, delrho_A, 507 & dneFSPdn, dneFSPdg, 1, 2, neFSP) 508 509 pbeup = neFSP 510 511c functional deriv info below fffffffffffff 512 if (pbeup.lt.pbe) then 513 etildea = pbe 514 detiladna = dneggadn(1) 515 detiladnb = dneggadn(2) 516 detiladgaa = dneggadg(D1_GAA) 517 detiladgbb = dneggadg(D1_GBB) 518 detiladgab = dneggadg(D1_GAB) 519 else 520 etildea = pbeup 521 detiladna = dneFSPdn(1) 522 detiladnb = 0.d0 523 detiladgaa = dneFSPdg(D1_GAA) 524 detiladgbb = 0.d0 525 detiladgab = 0.d0 526 endif 527 528c n_sigma/n_total factor 529 fabup=rhoa/rhoval 530 531 delc= xx2*etildea 532 delca = -(1.d0 + ccc)*fabup*delc 533 erevc = en + delca 534c 535c Beta bit 536c set up values to call PBE subroutine as 537c Fully SpinPolarized system for beta spin 538c to get E_GGA[rho_beta,0,gamma_beta,0] 539c 540 rho_A(1) = rhob 541 rho_A(2) = rhob 542 rho_A(3) = 0.d0 ! beta equals zero 543 delrho_A(1,1) = delrho_t(1,2) ! nabla n_up x 544 delrho_A(2,1) = delrho_t(2,2) ! nabla n_up y 545 delrho_A(3,1) = delrho_t(3,2) ! nabla n_up z 546 delrho_A(1,2) = 0.d0 ! set beta gradient to zero 547 delrho_A(2,2) = 0.d0 ! set beta gradient to zero 548 delrho_A(3,2) = 0.d0 ! set beta gradient to zero 549 550 neFSP = 0.0d0 !Ec in PBE 551 dneFSPdn(1) = 0.0d0 !Amat in PBE 552 dneFSPdn(2) = 0.0d0 !Amat in PBE 553 554 dneFSPdg(1) = 0.0d0 !Cmat in PBE 555 dneFSPdg(2) = 0.0d0 !Cmat in PBE 556 dneFSPdg(3) = 0.0d0 !Cmat in PBE 557c 558 call xc_cMpbe96(tol_rho, rho_A, delrho_A, 559 & dneFSPdn,dneFSPdg, 1, 2, neFSP) 560 561 pbedown = neFSP 562 563c functional deriv info below fffffffffffff 564 565 if (pbedown.lt.pbe) then 566 etildeb=pbe 567 detilbdna = dneggadn(1) 568 detilbdnb = dneggadn(2) 569 detilbdgaa = dneggadg(D1_GAA) 570 detilbdgbb = dneggadg(D1_GBB) 571 detilbdgab = dneggadg(D1_GAB) 572 else 573 etildeb = pbedown 574 detilbdna = 0.d0 575 detilbdnb = dneFSPdn(1) 576 detilbdgaa = 0.d0 577 detilbdgbb = dneFSPdg(D1_GAA) 578 detilbdgab = 0.d0 579 endif 580c 581c n_sigma/n_total factor 582 fabdown=rhob/rhoval 583 584 delc= xx2*etildeb 585 delcb = -(1.d0 + ccc)*fabdown*delc 586 erevc = erevc + delcb 587 588 erevsic = erevc*(1.d0 + dd*erevc*xx3) 589 590 if(ldew) func(n) = func(n) + rhoval*erevsic*fac 591 Ec = Ec + rhoval*erevsic*qwght(n)*fac 592c 593c na 594 atermn=pbe*ccc*dtwt2dn+(1.d0+ccc*xx2)*dneggadn(1)+ 595 & pbe*xx2*dcccdna 596 btermn=(1.d0+ccc)*( 597 & dtwt2dn*(fabup*etildea+fabdown*etildeb) + 598 + xx2*( (etildea - etildeb)*rhob/(rhoval*rhoval) + 599 + fabup*detiladna + fabdown*detilbdna) ) + 600 Z xx2*(fabup*etildea+fabdown*etildeb)*dcccdna 601 drevdna=atermn-btermn 602c 603c nb 604 atermn=pbe*ccc*dtwt2dn+(1.d0+ccc*xx2)*dneggadn(2)+ 605 & pbe*xx2*dcccdnb 606 btermn=(1.d0+ccc)*( 607 & dtwt2dn*(fabup*etildea+fabdown*etildeb) + 608 & xx2*( (etildeb-etildea)*rhoa/(rhoval*rhoval) + 609 & fabup*detilbdna+fabdown*detilbdnb) )+ 610 & xx2*(fabup*etildea+fabdown*etildeb)*dcccdnb 611 drevdnb=atermn-btermn 612c 613c gaa 614 atermg=pbe*ccc*dtwt2dg+(1.d0+ccc*xx2)*dneggadg(D1_GAA)+ 615 & pbe*xx2*dcccdgaa 616 btermg=(1.d0+ccc)* 617 & (xx2*(fabup*detiladgaa+fabdown*detilbdgaa)+ 618 & (etildea*fabup+etildeb*fabdown)*dtwt2dg)+ 619 & xx2*(etildea*fabup+etildeb*fabdown)*dcccdgaa 620 drevdgaa=atermg-btermg 621c 622c gbb 623 atermg=pbe*ccc*dtwt2dg+(1.d0+ccc*xx2)*dneggadg(D1_GBB)+ 624 & pbe*xx2*dcccdgbb 625 btermg=(1.d0+ccc)* 626 & (xx2*(fabup*detiladgbb+fabdown*detilbdgbb)+ 627 & (etildea*fabup+etildeb*fabdown)*dtwt2dg)+ 628 & xx2*(etildea*fabup+etildeb*fabdown)*dcccdgbb 629 drevdgbb=atermg-btermg 630c 631c gab 632 atermg=pbe*ccc*2.d0*dtwt2dg+(1.d0+ccc*xx2)*dneggadg(D1_GAB)+ 633 & pbe*xx2*dcccdgab 634 btermg=(1.d0+ccc)* 635 & (xx2*(fabup*detiladgab+fabdown*detilbdgab)+ 636 & (etildea*fabup+etildeb*fabdown)*2.d0*dtwt2dg)+ 637 & xx2*(etildea*fabup+etildeb*fabdown)*dcccdgab 638 drevdgab=atermg-btermg 639c 640c ta,tb 641 if(abs(taun).gt.tol_rho) then 642 drevdta=-2.d0*xx2/tauN* 643 *(ccc*pbe-(1.d0+ccc)*(fabup*etildea+fabdown*etildeb)) 644 drevdtb=-2d0*xx2/tauN* 645 * (ccc*pbe-(1.d0+ccc)*(fabup*etildea+fabdown*etildeb)) 646 else 647 drevdta=0d0 648 drevdtb=0d0 649 endif 650c 651c derivs wrt na 652 finaln= rhoval*drevdna + erevc + 653 & dd*(erevc*erevc*xx3 + rhoval*(xx3*2.d0*erevc*drevdna 654 & + erevc*erevc*dtwt3dn)) 655 Amat(n,1)=Amat(n,1)+finaln*fac 656c 657c derivs wrt nb 658 finaln= rhoval*drevdnb + erevc + 659 & dd*(erevc*erevc*xx3 + rhoval*(xx3*2.d0*erevc*drevdnb 660 & + erevc*erevc*dtwt3dn)) 661 Amat(n,2)=Amat(n,2)+finaln*fac 662c 663c derivs wrt gaa 664 apartg=rhoval*drevdgaa 665 cpartg(1)=erevc*erevc*dtwt3dg 666 cpartg(2)=xx3*2.d0*erevc*drevdgaa 667 finalgaa=apartg+rhoval*dd*(cpartg(1)+cpartg(2)) 668 Cmat(n,D1_GAA)=Cmat(n,D1_GAA)+finalgaa*fac 669c 670c derivs wrt gbb 671 apartg=rhoval*drevdgbb 672 cpartg(1)=erevc*erevc*dtwt3dg 673 cpartg(2)=xx3*2.d0*erevc*drevdgbb 674 finalgbb=apartg+rhoval*dd*(cpartg(1)+cpartg(2)) 675 Cmat(n,D1_GBB)=Cmat(n,D1_GBB)+finalgbb*fac 676c 677c derivs wrt gab 678 apartg=rhoval*drevdgab 679 cpartg(1)=erevc*erevc*2.d0*dtwt3dg 680 cpartg(2)=xx3*2.d0*erevc*drevdgab 681 finalgab=apartg+rhoval*dd*(cpartg(1)+cpartg(2)) 682 Cmat(n,D1_GAB)=Cmat(n,D1_GAB)+finalgab*fac 683c 684c derivs wrt ta 685 apartt=rhoval*drevdta 686 bpartt(1)=erevc*erevc*dtwt3dt 687 bpartt(2)=xx3*2.d0*erevc*drevdta 688 finalt=apartt+dd*rhoval*(bpartt(1)+bpartt(2)) 689 Mmat(n,1)=Mmat(n,1)+0.5d0*finalt*fac 690c 691c derivs wrt tb 692 apartt=rhoval*drevdtb 693 bpartt(1)=erevc*erevc*dtwt3dt 694 bpartt(2)=xx3*2.d0*erevc*drevdtb 695 finalt=apartt+dd*rhoval*(bpartt(1)+bpartt(2)) 696 Mmat(n,2)=Mmat(n,2)+0.5d0*finalt*fac 697c 69820 continue 699 700 endif !end ipol=2 case 701 702 return 703 end 704c 705 Subroutine xc_ctpss03_d2() 706 call errquit(' ctpss03: d2 not coded ',0,0) 707 return 708 end 709