1* ************************************************ 2* * * 3* * nwpw_tpss03_x * 4* * * 5* ************************************************ 6 subroutine nwpw_tpss03_x(pi,thrd,twthrd,frthrd,fvthrd,etthrd, 7 > C920,C1,C2,C3,kappa,mu,b,c,e, 8 > Cx,P23,es, 9 > n,agr,tau, 10 > xe,dfdnx,dfdagrx,dfdtaux) 11 implicit none 12* ***** input ***** 13 real*8 pi,thrd,twthrd,frthrd,fvthrd,etthrd,C920,C1,C2,C3 14 real*8 kappa,mu,b,c,e 15 real*8 Cx,P23,es 16 real*8 n,agr,tau 17* ***** output ***** 18 real*8 xe,dfdnx,dfdagrx,dfdtaux 19* ***** local declarations ***** 20 real*8 n_13,n_53,n_83,inv_n,agr2,tauW,tauU 21 real*8 p,z,fz,alpha,qb,p2,p3,z2,z3,qb2,fa,thresA,dalpha_dfa 22 real*8 x1a,x1,x2,x3a,x3b,x3,x4,x5,x6,x7,xnum,x 23 real*8 dp_dn,dz_dn,dalpha_dn,qba,qbb,dqba_dn,dqbb_dn,dqb_dn 24 real*8 dx1_dn,dx2_dn,dx3_dn,dx4_dn,dx5_dn,dx6_dn,dx7_dn 25 real*8 dx3a_dn,dx3b_dn,dxnum_dn,dx_dn 26 real*8 dp_dagr,dz_dagr,dalpha_dagr,dqba_dagr,dqbb_dagr,dqb_dagr 27 real*8 dx1_dagr,dx2_dagr,dx3_dagr,dx4_dagr,dx5_dagr,dx6_dagr 28 real*8 dx7_dagr,dx3a_dagr,dx3b_dagr,dxnum_dagr,dx_dagr 29 real*8 dz_dtau,dalpha_dtau,dqb_dtau 30 real*8 dx1_dtau,dx2_dtau,dx3_dtau,dx5_dtau,dxnum_dtau,dx_dtau 31 real*8 Fx,ex0,nex0,dFx_dx,dFx_dn,dFx_dagr,dFx_dtau 32 33 34 n_13 = n**thrd 35 n_53 = n_13*n_13*n 36 n_83 = n_53*n 37 inv_n = 1.0d0/n 38 agr2 = agr*agr 39 40 p = agr2/(4.0d0*P23*n_83) 41 dp_dn = -etthrd*p*inv_n 42 dp_dagr = 2.0d0*p/agr 43c dp_dtau = 0.0d0 44 45 p2 = p*p 46 p3 = p2*p 47 48 tauW = 0.125d0*agr2*inv_n 49 tauU = 0.3d0*P23*n_53 50 fz = tauW/tau 51 52 if (fz .gt. 1.0d0) then 53 z = 1.0d0 54 dz_dn = 0.0d0 55 dz_dagr = 0.0d0 56 dz_dtau = 0.0d0 57 else 58 z = fz 59 dz_dn = -z*inv_n 60 dz_dagr = 2.0d0*z/agr 61 dz_dtau = -z/tau 62 end if 63 64 z2 = z*z 65 z3 = z2*z 66 67 alpha = fvthrd*p*(1.0d0/z - 1.0d0) 68c alpha = (tau - tauW)/tauU 69 70 if (alpha .le. 0.0d0) then 71 alpha = 0.0d0 72 dalpha_dn = 0.0d0 73 dalpha_dagr = 0.0d0 74 dalpha_dtau = 0.0d0 75 else 76 dalpha_dn = fvthrd*(-p*dz_dn/z2 + dp_dn*(1.0d0/z - 1.0d0)) 77 dalpha_dagr = (alpha/p)*dp_dagr - fvthrd*(p/z2)*dz_dagr 78 dalpha_dtau = 1.0d0/tauU 79c dalpha_dtau = fvthrd*p*(-1.0d0/z2)*dz_dtau 80 end if 81 82 qb = C920*(alpha - 1.0d0)/dsqrt(1.0d0+b*alpha*(alpha - 1.0d0)) 83 & + twthrd*p 84 85 qb2 = qb*qb 86 87 x1a = (c*z2)/((1.0d0 + z2)**2.0d0) 88 x1 = (C1 + x1a)*p 89 x2 = C2*qb2 90 x3a = C3*qb 91 x3b = dsqrt(0.5d0*(0.36d0*z2 + p2)) 92 x3 = x3a*x3b 93 x4 = C1*C1*p2/kappa 94 x5 = 2.0d0*es*C1*0.36d0*z2 95 x6 = e*mu*p3 96 x7 = 1.0d0/((1.0d0 + es*p)**2.0d0) 97 98 xnum = x1 + x2 + x3 + x4 + x5 + x6 99 x = xnum*x7 100 101* ***** dxdn ***** 102 qba = alpha - 1.0d0 103 qbb = 1.0d0/dsqrt(1.0d0 + b*alpha*(alpha - 1.0d0)) 104 dqba_dn = dalpha_dn 105 dqbb_dn = b*dalpha_dn*(2.0d0*alpha - 1.0d0)* 106 & (-0.5d0)*(qbb**3.0d0) 107 dqb_dn = C920*(qbb*dqba_dn + qba*dqbb_dn) + twthrd*dp_dn 108 109 dx1_dn = (x1/p)*dp_dn 110 & + 2.0d0*c*p*z*dz_dn/((1.0d0 + z2)**3.0d0)*(1.0d0 - z2) 111 dx2_dn = 2.0d0*C2*qb*dqb_dn 112 113 dx3a_dn = C3*dqb_dn 114 dx3b_dn = 0.5d0/x3b*(0.36d0*z*dz_dn + p*dp_dn) 115 dx3_dn = x3b*dx3a_dn+x3a*dx3b_dn 116 117 dx4_dn = (2.0d0*x4/p)*dp_dn 118 dx5_dn = (2.0d0*x5/z)*dz_dn 119 dx6_dn = (3.0d0*x6/p)*dp_dn 120 dx7_dn = -2.0d0*es*dp_dn/(1.0d0 + es*p)**3.0d0 121 122 dxnum_dn = dx1_dn + dx2_dn + dx3_dn + dx4_dn + dx5_dn + dx6_dn 123 dx_dn = x7*dxnum_dn + xnum*dx7_dn 124 125* ***** dxdagr ***** 126 dqba_dagr = dalpha_dagr 127 dqbb_dagr = -0.5d0*dalpha_dagr*b*(2.0d0*alpha - 1.0d0) 128 & *qbb**3.0d0 129 dqb_dagr = C920*(qba*dqbb_dagr + qbb*dqba_dagr) 130 & + twthrd*dp_dagr 131 dx1_dagr = (x1/p)*dp_dagr 132 & + 2.0d0*c*p*z*dz_dagr/((1.0d0 + z2)**3.0d0)*(1.0d0 - z2) 133 134 dx2_dagr = C2*2.0d0*qb*dqb_dagr 135 136 dx3a_dagr = C3*dqb_dagr 137 dx3b_dagr = 0.5d0/x3b*( 0.36d0*z*dz_dagr + p*dp_dagr) 138 dx3_dagr = x3b*dx3a_dagr + x3a*dx3b_dagr 139 140 dx4_dagr = (2.0d0*x4/p)*dp_dagr 141 dx5_dagr = (2.0d0*x5/z)*dz_dagr 142 dx6_dagr = (3.0d0*x6/p)*dp_dagr 143 144 dx7_dagr = -2.0d0*es*dp_dagr/(1.0d0 + es*p)**3.0d0 145 146 dxnum_dagr = dx1_dagr + dx2_dagr + dx3_dagr 147 & + dx4_dagr + dx5_dagr + dx6_dagr 148 dx_dagr = x7*dxnum_dagr + xnum*dx7_dagr 149 150* ***** dxdtau ***** 151 152 dqb_dtau = C920*dalpha_dtau*qbb*(1.0d0 153 & - 0.5d0*b*qba*qbb*qbb*(2.0d0*alpha - 1.0d0)) 154 155 dx1_dtau = c*p*dz_dtau*2.0d0*z*(1.0d0 - z2)/ 156 & ((1.0d0 + z2)**3.0d0) 157 dx2_dtau = 2.0d0*c2*qb*dqb_dtau 158 dx3_dtau = x3*(dqb_dtau/qb 159 & + 0.5d0*0.36d0*z*dz_dtau/(x3b*x3b)) 160 dx5_dtau = 2.0d0*(x5/z)*dz_dtau 161 162 dxnum_dtau = dx1_dtau + dx2_dtau + dx3_dtau + dx5_dtau 163 dx_dtau = x7*dxnum_dtau 164 165 Fx = 1.0d0 + kappa - kappa/(1.0d0 + x/kappa) 166 167 ex0 = Cx*n_13 168 xe = ex0*Fx 169 nex0 = n*ex0 170 171 dFx_dx = 1.0d0/(1.0d0 + x/kappa)**2.0d0 172 dFx_dn = dFx_dx*dx_dn 173 dfdnx = frthrd*xe + nex0*dFx_dn 174 175 dFx_dagr = dFx_dx*dx_dagr 176 dfdagrx = nex0*dFx_dagr 177 178 dFx_dtau = dFx_dx*dx_dtau 179 dfdtaux = nex0*dFx_dtau 180 181 return 182 end 183 184* ************************************************ 185* * * 186* * gen_TPSS03_restricted * 187* * * 188* ************************************************ 189 190* This function returns the TPSS03 exchange-correlation 191* energy density, xce, and its derivatives with respect 192* to n, |grad n|, tau. 193 194* 195* Entry - n2ft3d : number of grid points 196* rho_in(*) : density (nup+ndn) 197* agr_in(*): |grad rho_in| 198* tau_in(*): tau 199* x_parameter: scale parameter for exchange 200* c_parameter: scale parameter for correlation 201* 202* Exit - xce(n2ft3d) : TPSS03 exchange correlation energy density 203* fn(n2ft3d) : d(n*xce)/dn 204* fdn(n2ft3d) : d(n*xce)/d|grad n| 205* fdtau(n2ft3d) : d(n*xce)/dtau 206* 207 208 subroutine gen_TPSS03_restricted(n2ft3d,rho_in,agr_in,tau_in, 209 > x_parameter,c_parameter, 210 > xce,fn,fdn,fdtau) 211 implicit none 212* ***** input ***** 213 integer n2ft3d 214 real*8 rho_in(*),agr_in(*),tau_in(*) 215 real*8 x_parameter,c_parameter 216* ***** output ***** 217 real*8 xce(*),fn(*),fdn(*),fdtau(*) 218* ***** local declarations ***** 219 integer i 220 real*8 n,agr,tau 221 real*8 eup0c,fnup0c,fdnup0c 222 real*8 e0c,fn0c,fdn0c 223 real*8 Cx,P23,es 224 real*8 ex,fnx,fdnx,fdtaux 225 real*8 z,z2,z3,tmp1,tmp2,agr2,c00 226 real*8 tauW,dz2_dn,dz2_dagr,dz2_dtau,dz3_dn,dz3_dagr,dz3_dtau 227 real*8 PKZB1,dPKZB1_dn,dPKZB1_dagr,dPKZB1_dtau 228 real*8 pbec,dpbec_dn,dpbec_dagr 229 real*8 pbeupc,dpbeupc_dn,dpbeupc_dagr 230 real*8 etil,detil_dn,detil_dagr 231 real*8 PKZB2,dPKZB2_dn,dPKZB2_dagr,dPKZB2_dtau 232 real*8 revPKZB,drev_dn,drev_dagr,drev_dtau 233 real*8 revz3,drevz3_dn,drevz3_dagr,drevz3_dtau 234 real*8 ec,fnc,fdnc,fdtauc 235* ***** constants ***** 236 real*8 pi,thrd,twthrd,frthrd,fvthrd,etthrd 237 parameter (pi = 3.14159265358979311599d0) 238 parameter (thrd = 1.0d0/3.0d0) 239 parameter (twthrd = 2.0d0/3.0d0) 240 parameter (frthrd = 4.0d0/3.0d0) 241 parameter (fvthrd = 5.0d0/3.0d0) 242 parameter (etthrd = 8.0d0/3.0d0) 243* ***** density cutoff parameters ***** 244 real*8 tol,ETA 245 parameter (tol = 1.0d-18) 246 parameter (ETA = 1.0d-20) 247* ***** TPSS03 constants ***** 248 real*8 kappa,mu,b,c,d,e,C920,C1,C2,C3 249 parameter (kappa = 0.8040d0) 250 parameter (mu = 0.21951d0) 251 parameter (b = 0.40d0) 252 parameter (c = 1.59096d0) 253 parameter (d = 2.8d0) 254 parameter (e = 1.537d0) 255 parameter (C920 = 9.0d0/20.0d0) 256 parameter (C1 = 10.0d0/81.0d0) 257 parameter (C2 = 146.0d0/2025.0d0) 258 parameter (C3 = -73.0d0/405.0d0) 259 260 Cx = (-0.75d0)*(3.0d0/pi)**thrd 261 P23 = (3.0d0*pi**2.0d0)**twthrd 262 es = dsqrt(e) 263 264 do i=1,n2ft3d 265 n = rho_in(i) + ETA 266 agr = agr_in(i) + ETA 267 tau = 2.0d0*tau_in(i) + ETA 268 269* ***** TPSS03 Exchange ***** 270 call nwpw_tpss03_x(pi,thrd,twthrd,frthrd,fvthrd,etthrd, 271 > C920,C1,C2,C3,kappa,mu,b,c,e, 272 > Cx,P23,es, 273 > n,agr,tau, 274 > ex,fnx,fdnx,fdtaux) 275 276* ***** TPSS03 Correlation ***** 277 call gen_PBE96_c_restricted(rho_in(i),agr_in(i), 278 > e0c,fn0c,fdn0c) 279 280 pbec = e0c 281 dpbec_dn = (fn0c - pbec)/n 282 dpbec_dagr = fdn0c/n 283 284 call gen_PBE96_c_unrestricted(0.50d0*rho_in(i),0.50d0*agr_in(i), 285 > eup0c,fnup0c,fdnup0c) 286 287 pbeupc = eup0c 288 dpbeupc_dn = (fnup0c - pbeupc)/n 289 dpbeupc_dagr = fdnup0c/n 290 291 c00 = 0.53d0 292 293 agr2 = agr*agr 294 295 tauW = 0.125d0*agr2/n 296 z = tauW/tau 297 298 if (z .gt. 1.0d0) then 299 z2 = 1.0d0 300 dz2_dn = 0.0d0 301 dz2_dagr = 0.0d0 302 dz2_dtau = 0.0d0 303 304 z3 = 1.0d0 305 dz3_dn = 0.0d0 306 dz3_dagr = 0.0d0 307 dz3_dtau = 0.0d0 308 else 309 z2 = z*z 310 dz2_dn = -2.0d0*z2/n 311 dz2_dagr = 4.0d0*z2/agr 312 dz2_dtau = -2.0d0*z2/tau 313 314 z3 = z2*z 315 dz3_dn = -3.0d0*z3/n 316 dz3_dagr = 6.0d0*z3/agr 317 dz3_dtau = -3.0d0*z3/tau 318 end if 319 320 tmp1 = 1.0d0 + c00*z2 321 PKZB1 = pbec*tmp1 322 dPKZB1_dn = dpbec_dn*tmp1 + pbec*c00*dz2_dn 323 dPKZB1_dagr = dpbec_dagr*tmp1 + pbec*c00*dz2_dagr 324 dPKZB1_dtau = pbec*c00*dz2_dtau 325 326 if (pbeupc .lt. pbec) then 327 etil = pbec 328 detil_dn = dpbec_dn 329 detil_dagr = dpbec_dagr 330 else 331 etil = pbeupc 332 detil_dn = dpbeupc_dn 333 detil_dagr = dpbeupc_dagr 334 endif 335 336 tmp1 = -(1.0d0 + c00) 337 tmp2 = etil 338 PKZB2 = tmp1*z2*tmp2 339 dPKZB2_dn = tmp1*(dz2_dn*tmp2 + z2*detil_dn) 340 dPKZB2_dagr = tmp1*(dz2_dagr*tmp2 + z2*detil_dagr) 341 dPKZB2_dtau = tmp1*dz2_dtau*tmp2 342 343 revPKZB = PKZB1 + PKZB2 344 drev_dn = dPKZB1_dn + dPKZB2_dn 345 drev_dagr = dPKZB1_dagr + dPKZB2_dagr 346 drev_dtau = dPKZB1_dtau + dPKZB2_dtau 347 348 revz3 = 1.0d0 + d*revPKZB*z3 349 drevz3_dn = d*(drev_dn*z3 + revPKZB*dz3_dn) 350 drevz3_dagr = d*(drev_dagr*z3 + revPKZB*dz3_dagr) 351 drevz3_dtau = d*(drev_dtau*z3 + revPKZB*dz3_dtau) 352 353 ec = revPKZB*revz3 354 fnc = n*(drev_dn*revz3 + revPKZB*drevz3_dn) + ec 355 fdnc = n*(drev_dagr*revz3 + revPKZB*drevz3_dagr) 356 fdtauc = n*(drev_dtau*revz3 + revPKZB*drevz3_dtau) 357 358 xce(i) = x_parameter*ex + c_parameter*ec 359 fn(i) = x_parameter*fnx + c_parameter*fnc 360 fdn(i) = x_parameter*fdnx + c_parameter*fdnc 361 fdtau(i) = x_parameter*fdtaux + c_parameter*fdtauc 362 363 end do 364 return 365 end 366 367* ************************************************ 368* * * 369* * gen_TPSS03_unrestricted * 370* * * 371* ************************************************ 372 373* This function returns the TPSS03 exchange-correlation 374* energy density, xce, and its derivatives with respect 375* to nup, ndn, |grad nup|, |grad ndn|, |grad n|, tauup, taudn. 376 377* 378* Entry - n2ft3d : number of grid points 379* rho_in(*,2) : density (nup and ndn) 380* agr_in(*,3): |grad rho_in| (nup, ndn and n) 381* tau_in(*,2): tau (nup and ndn) 382* x_parameter: scale parameter for exchange 383* c_parameter: scale parameter for correlation 384* 385* Exit - xce(n2ft3d) : TPSS03 exchange correlation energy density 386* fn(n2ft3d,2) : d(n*xce)/dnup, d(n*xce)/dndn 387* fdn(n2ft3d,3) : d(n*xce)/d|grad nup|, d(n*xce)/d|grad ndn|, d(n*xce)/d|grad n| 388* fdtau(n2ft3d,2) : d(n*xce)/dtauup, d(n*xce)/dtaudn 389* 390 391 subroutine gen_TPSS03_unrestricted(n2ft3d,rho_in,agr_in,tau_in, 392 > x_parameter,c_parameter, 393 > xce,fn,fdn,fdtau) 394 implicit none 395* ***** input ***** 396 integer n2ft3d 397 real*8 rho_in(n2ft3d,2),agr_in(n2ft3d,3),tau_in(n2ft3d,2) 398 real*8 x_parameter,c_parameter 399* ***** output ***** 400 real*8 xce(n2ft3d),fn(n2ft3d,2),fdn(n2ft3d,3),fdtau(n2ft3d,2) 401* ***** local declarations ***** 402 integer i 403 real*8 n,agr,tau 404 real*8 nup,agrup,tauup 405 real*8 ndn,agrdn,taudn 406 real*8 Cx,P23,es 407 real*8 z,z2,z3 408 real*8 ex,eupx,fnupx,fdnupx,fdtauupx,ednx,fndnx,fdndnx,fdtaudnx 409 real*8 e0c,fn0c1,fn0c2,fdn0c1,fdn0c2,fdn0c3 410 real*8 eup0c,fnup0c,fdnup0c 411 real*8 edn0c,fndn0c,fdndn0c 412 real*8 zeta,dzeta_dnup,dzeta_dndn 413 real*8 onepzeta,onemzeta,onepzeta2,onemzeta2 414 real*8 zeta2,one2mzeta2 415 real*8 n2,agr2,agraa,agrbb,agrab 416 real*8 tmp1,tmp2,denxi,denxi2,gzeta2 417 real*8 dgzeta2_dnup,dgzeta2_dndn 418 real*8 dgzeta2_dagrup,dgzeta2_dagrdn,dgzeta2_dagr 419 real*8 xi2,dxi2_dnup,dxi2_dndn 420 real*8 dxi2_dagrup,dxi2_dagrdn,dxi2_dagr 421 real*8 cz0,dcz0_dzeta,onezetap43,onezetam73 422 real*8 denczx,ddenczx_dnup,ddenczx_dndn,denczx4,denczx5,denczx8 423 real*8 czx,dczx_dnup,dczx_dndn 424 real*8 dczx_dagrup,dczx_dagrdn,dczx_dagr 425 real*8 tauW,dz2_dn,dz2_dagr,dz2_dtau,dz3_dn,dz3_dagr,dz3_dtau 426 real*8 PKZB1,dPKZB1_dnup,dPKZB1_dndn 427 real*8 dPKZB1_dagrup,dPKZB1_dagrdn,dPKZB1_dagr 428 real*8 dPKZB1_dtauup,dPKZB1_dtaudn 429 real*8 pbec,dpbec_dnup,dpbec_dndn 430 real*8 dpbec_dagrup,dpbec_dagrdn,dpbec_dagr 431 real*8 pbeupc,dpbeupc_dnup,dpbeupc_dndn 432 real*8 dpbeupc_dagrup,dpbeupc_dagrdn,dpbeupc_dagr 433 real*8 pbednc,dpbednc_dnup,dpbednc_dndn 434 real*8 dpbednc_dagrup,dpbednc_dagrdn,dpbednc_dagr 435 real*8 etilup,detilup_dnup,detilup_dndn 436 real*8 detilup_dagrup,detilup_dagrdn,detilup_dagr 437 real*8 etildn,detildn_dnup,detildn_dndn 438 real*8 detildn_dagrup,detildn_dagrdn,detildn_dagr 439 real*8 fa,fb 440 real*8 PKZB2,dPKZB2_dnup,dPKZB2_dndn 441 real*8 dPKZB2_dagrup,dPKZB2_dagrdn,dPKZB2_dagr 442 real*8 dPKZB2_dtauup,dPKZB2_dtaudn 443 real*8 revPKZB,drev_dnup,drev_dndn 444 real*8 drev_dagrup,drev_dagrdn,drev_dagr 445 real*8 drev_dtauup,drev_dtaudn 446 real*8 revz3,drevz3_dnup,drevz3_dndn 447 real*8 drevz3_dagrup,drevz3_dagrdn,drevz3_dagr 448 real*8 drevz3_dtauup,drevz3_dtaudn 449 real*8 ec,fnupc,fndnc,fdnupc,fdndnc,fdnc,fdtauupc,fdtaudnc 450* ***** constants ***** 451 real*8 pi,thrd,twthrd,frthrd,fvthrd,svthrd,etthrd 452 parameter (pi = 3.14159265358979311599d0) 453 parameter (thrd = 1.0d0/3.0d0) 454 parameter (twthrd = 2.0d0/3.0d0) 455 parameter (frthrd = 4.0d0/3.0d0) 456 parameter (fvthrd = 5.0d0/3.0d0) 457 parameter (svthrd = 7.0d0/3.0d0) 458 parameter (etthrd = 8.0d0/3.0d0) 459* ***** density cutoff parameters ***** 460 real*8 tol,ETA 461 parameter (tol = 1.0d-18) 462 parameter (ETA = 1.0d-20) 463* ***** TPSS03 constants ***** 464 real*8 kappa,mu,b,c,d,e,C920,C1,C2,C3 465 parameter (kappa = 0.8040d0) 466 parameter (mu = 0.21951d0) 467 parameter (b = 0.40d0) 468 parameter (c = 1.59096d0) 469 parameter (d = 2.8d0) 470 parameter (e = 1.537d0) 471 parameter (C920 = 9.0d0/20.0d0) 472 parameter (C1 = 10.0d0/81.0d0) 473 parameter (C2 = 146.0d0/2025.0d0) 474 parameter (C3 = -73.0d0/405.0d0) 475 476 Cx = (-0.75d0)*(3.0d0/pi)**thrd 477 P23 = (3.0d0*pi**2.0d0)**twthrd 478 es = dsqrt(e) 479 480 do i=1,n2ft3d 481 nup = rho_in(i,1) + ETA 482 agrup = agr_in(i,1) + ETA 483 tauup = tau_in(i,1) + ETA 484 ndn = rho_in(i,2) + ETA 485 agrdn = agr_in(i,2) + ETA 486 taudn = tau_in(i,2) + ETA 487 488 489* ***** TPSS03 Exchange ***** 490* ***** UP ***** 491 n = 2.0d0*nup 492 agr = 2.0d0*agrup 493 tau = 2.0d0*tauup 494 495 call nwpw_tpss03_x(pi,thrd,twthrd,frthrd,fvthrd,etthrd, 496 > C920,C1,C2,C3,kappa,mu,b,c,e, 497 > Cx,P23,es, 498 > n,agr,tau,eupx,fnupx,fdnupx,fdtauupx) 499 500* ***** DOWN ***** 501 n = 2.0d0*ndn 502 agr = 2.0d0*agrdn 503 tau = 2.0d0*taudn 504 505 call nwpw_tpss03_x(pi,thrd,twthrd,frthrd,fvthrd,etthrd, 506 > C920,C1,C2,C3,kappa,mu,b,c,e, 507 > Cx,P23,es, 508 > n,agr,tau,ednx,fndnx,fdndnx,fdtaudnx) 509 510 n = nup + ndn 511 512 ex = (eupx*nup + ednx*ndn)/n 513 514* ***** TPSS03 Correlation ***** 515 agr = agr_in(i,3) + ETA 516 517 tau = tauup + taudn 518 n2 = n*n 519 agr2 = agr*agr 520 521 call gen_PBE96_c_full_unrestricted(rho_in(i,1),rho_in(i,2), 522 > agr_in(i,1),agr_in(i,2),agr_in(i,3), 523 > e0c,fn0c1,fn0c2,fdn0c1,fdn0c2,fdn0c3) 524 525 pbec = e0c 526 dpbec_dnup = (fn0c1 - pbec)/n 527 dpbec_dndn = (fn0c2 - pbec)/n 528 dpbec_dagrup = fdn0c1/n 529 dpbec_dagrdn = fdn0c2/n 530 dpbec_dagr = fdn0c3/n 531 532 call gen_PBE96_c_unrestricted(rho_in(i,1),agr_in(i,1), 533 > eup0c,fnup0c,fdnup0c) 534 535 pbeupc = eup0c 536 dpbeupc_dnup = (fnup0c - pbeupc)/nup 537 dpbeupc_dndn = 0.0d0 538 dpbeupc_dagrup = fdnup0c/nup 539 dpbeupc_dagrdn = 0.0d0 540 dpbeupc_dagr = 0.0d0 541 542 call gen_PBE96_c_unrestricted(rho_in(i,2),agr_in(i,2), 543 > edn0c,fndn0c,fdndn0c) 544 545 pbednc = edn0c 546 dpbednc_dnup = 0.0d0 547 dpbednc_dndn = (fndn0c - pbednc)/ndn 548 dpbednc_dagrup = 0.0d0 549 dpbednc_dagrdn = fdndn0c/ndn 550 dpbednc_dagr = 0.0d0 551 552c if (dabs(nup - ndn) .lt. 1.0d-7) then 553c czx = 0.53d0 554c dczx_dnup = 0.0d0 555c dczx_dndn = 0.0d0 556c dczx_dagrup = 0.0d0 557c dczx_dagrdn = 0.0d0 558c dczx_dagr = 0.0d0 559c else 560 zeta = (nup - ndn)/n 561 562 onepzeta = 1.0d0 + zeta 563 onemzeta = 1.0d0 - zeta 564 onepzeta2 = onepzeta*onepzeta 565 onemzeta2 = onemzeta*onemzeta 566 567 zeta2 = zeta*zeta 568 one2mzeta2 = 1.0d0 - zeta2 569 570 dzeta_dnup = onemzeta/n 571 dzeta_dndn = -onepzeta/n 572 573 agraa = agrup*agrup 574 agrbb = agrdn*agrdn 575 agrab = agraa + agrbb - agr2 576 577 denxi = 2.0d0*((3.0d0*pi*pi*n)**thrd) 578 denxi2 = denxi*denxi 579 580 gzeta2 = (agraa*onemzeta2 + agrbb*onepzeta2 581 & + agrab*one2mzeta2)/n2 582 583 tmp1 = (-agraa*onemzeta + agrbb*onepzeta 584 & - agrab*zeta)*2.0d0*dzeta_dnup 585 tmp1 = tmp1/n2 586 tmp2 = -2.0d0*gzeta2/n 587 dgzeta2_dnup = tmp1 + tmp2 588 589 tmp1 = (-agraa*onemzeta + agrbb*onepzeta 590 & - agrab*zeta)*2.0d0*dzeta_dndn 591 tmp1 = tmp1/n2 592 tmp2 = -2.0d0*gzeta2/n 593 dgzeta2_dndn = tmp1 + tmp2 594 595 dgzeta2_dagrup = 2.0d0*agrup*(onemzeta2 + one2mzeta2)/n2 596 dgzeta2_dagrdn = 2.0d0*agrdn*(onepzeta2 + one2mzeta2)/n2 597 dgzeta2_dagr = -2.0d0*agr*one2mzeta2/n2 598 599 xi2 = gzeta2/denxi2 600 601 dxi2_dnup = (dgzeta2_dnup - twthrd*gzeta2/n)/denxi2 602 dxi2_dndn = (dgzeta2_dndn - twthrd*gzeta2/n)/denxi2 603 604 dxi2_dagrup = dgzeta2_dagrup/denxi2 605 dxi2_dagrdn = dgzeta2_dagrdn/denxi2 606 dxi2_dagr = dgzeta2_dagr/denxi2 607 608 cz0 = 0.53d0 + 0.87d0*zeta**2.0d0 609 & + 0.50d0*zeta**4.0d0 + 2.26*zeta**6.0d0 610 dcz0_dzeta = 1.74d0*zeta + 2.0d0*zeta**3.0d0 611 & + 13.56d0*zeta**5.0d0 612 onezetap43 = onepzeta**(-frthrd) + onemzeta**(-frthrd) 613 onezetam73 = onemzeta**(-svthrd) - onepzeta**(-svthrd) 614 615 denczx = 1.0d0 + 0.50d0*xi2*onezetap43 616 ddenczx_dnup = 2.0d0*(denczx**3.0d0)*(dxi2_dnup*onezetap43 617 & + xi2*frthrd*onezetam73*dzeta_dnup) 618 ddenczx_dndn = 2.0d0*(denczx**3.0d0)*(dxi2_dndn*onezetap43 619 & + xi2*frthrd*onezetam73*dzeta_dndn) 620 621 denczx4 = denczx**4.0d0 622 denczx5 = denczx4*denczx 623 denczx8 = denczx4*denczx4 624 625 czx = cz0/denczx4 626 dczx_dnup = dcz0_dzeta*dzeta_dnup/denczx4 627 & - cz0*ddenczx_dnup/denczx8 628 dczx_dndn = dcz0_dzeta*dzeta_dndn/denczx4 629 & - cz0*ddenczx_dndn/denczx8 630 dczx_dagrup = -2.0d0*cz0/denczx5*onezetap43*dxi2_dagrup 631 dczx_dagrdn = -2.0d0*cz0/denczx5*onezetap43*dxi2_dagrdn 632 dczx_dagr = -2.0d0*cz0/denczx5*onezetap43*dxi2_dagr 633c end if 634 635 tauW = 0.125d0*agr2/n 636 z = tauW/tau 637 638 if (z .gt. 1.0d0) then 639 z2 = 1.0d0 640 dz2_dn = 0.0d0 641 dz2_dagr = 0.0d0 642 dz2_dtau = 0.0d0 643 644 z3 = 1.0d0 645 dz3_dn = 0.0d0 646 dz3_dagr = 0.0d0 647 dz3_dtau = 0.0d0 648 else 649 z2 = z*z 650 dz2_dn = -2.0d0*z2/n 651c dz2_dnup == dz2_dndn == dz2_dn 652 dz2_dagr = 4.0d0*z2/agr 653 dz2_dtau = -2.0d0*z2/tau 654c dz2_dtauup == dz2_dtaudn == dz2_dtau 655 656 z3 = z2*z 657 dz3_dn = -3.0d0*z3/n 658c dz3_dnup == dz3_dndn == dz3_dn 659 dz3_dagr = 6.0d0*z3/agr 660 dz3_dtau = -3.0d0*z3/tau 661c dz3_dtauup == dz3_dtaudn == dz3_dtau 662 end if 663 664 tmp1 = 1.0d0 + czx*z2 665 PKZB1 = pbec*tmp1 666 dPKZB1_dnup = dpbec_dnup*tmp1 667 & + pbec*(dczx_dnup*z2 + czx*dz2_dn) 668 dPKZB1_dndn = dpbec_dndn*tmp1 669 & + pbec*(dczx_dndn*z2 + czx*dz2_dn) 670 dPKZB1_dagrup = dpbec_dagrup*tmp1 671 & + pbec*dczx_dagrup*z2 672 dPKZB1_dagrdn = dpbec_dagrdn*tmp1 673 & + pbec*dczx_dagrdn*z2 674 dPKZB1_dagr = dpbec_dagr*tmp1 675 & + pbec*(dczx_dagr*z2 + czx*dz2_dagr) 676 dPKZB1_dtauup = pbec*czx*dz2_dtau 677 dPKZB1_dtaudn = dPKZB1_dtauup 678 679 if (pbeupc .lt. pbec) then 680 etilup = pbec 681 detilup_dnup = dpbec_dnup 682 detilup_dndn = dpbec_dndn 683 detilup_dagrup = dpbec_dagrup 684 detilup_dagrdn = dpbec_dagrdn 685 detilup_dagr = dpbec_dagr 686 else 687 etilup = pbeupc 688 detilup_dnup = dpbeupc_dnup 689 detilup_dndn = dpbeupc_dndn 690 detilup_dagrup = dpbeupc_dagrup 691 detilup_dagrdn = dpbeupc_dagrdn 692 detilup_dagr = dpbeupc_dagr 693 endif 694 695 if (pbednc .lt. pbec) then 696 etildn = pbec 697 detildn_dnup = dpbec_dnup 698 detildn_dndn = dpbec_dndn 699 detildn_dagrup = dpbec_dagrup 700 detildn_dagrdn = dpbec_dagrdn 701 detildn_dagr = dpbec_dagr 702 else 703 etildn = pbednc 704 detildn_dnup = dpbednc_dnup 705 detildn_dndn = dpbednc_dndn 706 detildn_dagrup = dpbednc_dagrup 707 detildn_dagrdn = dpbednc_dagrdn 708 detildn_dagr = dpbednc_dagr 709 endif 710 711 fa = nup/n 712 fb = ndn/n 713 714 tmp1 = -(1.0d0 + czx) 715 tmp2 = fa*etilup + fb*etildn 716 PKZB2 = tmp1*z2*tmp2 717 dPKZB2_dnup = -dczx_dnup*z2*tmp2 + tmp1*(dz2_dn*tmp2 718 & + z2*(ndn/n2*(etilup - etildn) 719 & + fa*detilup_dnup + fb*detildn_dnup)) 720 dPKZB2_dndn = -dczx_dndn*z2*tmp2 + tmp1*(dz2_dn*tmp2 721 & + z2*(nup/n2*(etildn - etilup) 722 & + fb*detildn_dndn + fa*detilup_dndn)) 723 dPKZB2_dagrup = -dczx_dagrup*z2*tmp2 724 & + tmp1*z2*(fa*detilup_dagrup + fb*detildn_dagrup) 725 dPKZB2_dagrdn = -dczx_dagrdn*z2*tmp2 726 & + tmp1*z2*(fb*detildn_dagrdn + fa*detilup_dagrdn) 727 dPKZB2_dagr = -dczx_dagr*z2*tmp2 + tmp1*(dz2_dagr*tmp2 728 & + z2*(fa*detilup_dagr + fb*detildn_dagr)) 729 dPKZB2_dtauup = tmp1*dz2_dtau*tmp2 730 dPKZB2_dtaudn = dPKZB2_dtauup 731 732 revPKZB = PKZB1 + PKZB2 733 drev_dnup = dPKZB1_dnup + dPKZB2_dnup 734 drev_dndn = dPKZB1_dndn + dPKZB2_dndn 735 drev_dagrup = dPKZB1_dagrup + dPKZB2_dagrup 736 drev_dagrdn = dPKZB1_dagrdn + dPKZB2_dagrdn 737 drev_dagr = dPKZB1_dagr + dPKZB2_dagr 738 drev_dtauup = dPKZB1_dtauup + dPKZB2_dtauup 739 drev_dtaudn = dPKZB1_dtaudn + dPKZB2_dtaudn 740 741 revz3 = 1.0d0 + d*revPKZB*z3 742 drevz3_dnup = d*(drev_dnup*z3 + revPKZB*dz3_dn) 743 drevz3_dndn = d*(drev_dndn*z3 + revPKZB*dz3_dn) 744 drevz3_dagrup = d*drev_dagrup*z3 745 drevz3_dagrdn = d*drev_dagrdn*z3 746 drevz3_dagr = d*(drev_dagr*z3 + revPKZB*dz3_dagr) 747 drevz3_dtauup = d*(drev_dtauup*z3 + revPKZB*dz3_dtau) 748 drevz3_dtaudn = d*(drev_dtaudn*z3 + revPKZB*dz3_dtau) 749 750 ec = revPKZB*revz3 751 fnupc = n*(drev_dnup*revz3 + revPKZB*drevz3_dnup) + ec 752 fndnc = n*(drev_dndn*revz3 + revPKZB*drevz3_dndn) + ec 753 fdnupc = n*(drev_dagrup*revz3 + revPKZB*drevz3_dagrup) 754 fdndnc = n*(drev_dagrdn*revz3 + revPKZB*drevz3_dagrdn) 755 fdnc = n*(drev_dagr*revz3 + revPKZB*drevz3_dagr) 756 fdtauupc = n*(drev_dtauup*revz3 + revPKZB*drevz3_dtauup) 757 fdtaudnc = n*(drev_dtaudn*revz3 + revPKZB*drevz3_dtaudn) 758 759 760 xce(i) = x_parameter*ex + c_parameter*ec 761 fn(i,1) = x_parameter*fnupx + c_parameter*fnupc 762 fn(i,2) = x_parameter*fndnx + c_parameter*fndnc 763 fdn(i,1) = x_parameter*fdnupx + c_parameter*fdnupc 764 fdn(i,2) = x_parameter*fdndnx + c_parameter*fdndnc 765 fdn(i,3) = c_parameter*fdnc 766 fdtau(i,1) = x_parameter*fdtauupx + c_parameter*fdtauupc 767 fdtau(i,2) = x_parameter*fdtaudnx + c_parameter*fdtaudnc 768 769 end do 770 771 return 772 end 773