1* ************************************************ 2* * * 3* * nwpw_GVT4 * 4* * * 5* ************************************************ 6 subroutine nwpw_GVT4(ag,bg,cg,dg,eg,fg,alphag,alphagt, 7 > xg,zg,gammag,G,dGdx,dGdz) 8 implicit none 9 real*8 ag,bg,cg,dg,eg,fg,alphag,alphagt,xg,zg,gammag 10 real*8 G,dGdx,dGdz 11 real*8 xg2,zg2,gammag2 12 13 xg2 = xg*xg 14 zg2 = zg*zg 15 gammag2 = gammag*gammag 16 17 G = (ag + bg*xg + cg*zg 18 & + dg*xg2 + eg*xg*zg + fg*zg2)/gammag 19 20 dGdx = ( -ag*alphag 21 & + bg*(1.0d0 - 2.0d0*alphag*xg) 22 & - 2.0d0*cg*alphag*zg 23 & + dg*(2.0d0*xg - 3.0d0*alphag*xg2) 24 & + eg*(zg - 3.0d0*alphag*zg*xg) 25 & - 3.0d0*fg*alphag*zg2)/gammag2 26 27 dGdz = ( -ag*alphagt 28 & - 2.0d0*bg*alphagt*xg 29 & + cg*(1.0d0 - 2.0d0*alphagt*zg) 30 & - 3.0d0*dg*alphagt*xg2 31 & + eg*(xg - 3.0d0*alphagt*xg*zg) 32 & + fg*(2.0d0*zg - 3.0d0*alphagt*zg2))/gammag2 33 34 return 35 end 36* ************************************************ 37* * * 38* * gen_VS98_restricted * 39* * * 40* ************************************************ 41 42* This function returns the VS98 exchange-correlation 43* energy density, xce, and its derivatives with respect 44* to n, |grad n|, tau. 45 46* 47* Entry - n2ft3d : number of grid points 48* rho_in(*) : density (nup+ndn) 49* agr_in(*): |grad rho_in| 50* tau_in(*): tau 51* x_parameter: scale parameter for exchange 52* c_parameter: scale parameter for correlation 53* 54* Exit - xce(n2ft3d) : VS98 exchange correlation energy density 55* fn(n2ft3d) : d(n*xce)/dn 56* fdn(n2ft3d) : d(n*xce)/d|grad n| 57* fdtau(n2ft3d) : d(n*xce)/dtau 58* 59 subroutine gen_VS98_restricted(n2ft3d,rho_in,agr_in,tau_in, 60 > x_parameter,c_parameter, 61 > xce,fn,fdn,fdtau) 62 implicit none 63* ***** input ***** 64 integer n2ft3d 65 real*8 rho_in(*),agr_in(*),tau_in(*) 66 real*8 x_parameter,c_parameter 67* ***** output ***** 68 real*8 xce(*),fn(*),fdn(*),fdtau(*) 69* ***** local declarations ***** 70 integer i 71 real*8 n,agr,tau 72 real*8 xr,zr,gamma 73 real*8 inv_n,n_13,n_43,n_53,n_83,agr2 74 real*8 x,dx_dn,dx_dagr,z,dz_dn,dz_dtau 75 real*8 GG,dGdx,dGdz,dG_dn,dG_dagr,dG_dtau 76 real*8 ex,fnx,fdnx,fdtaux 77 real*8 ess0c,dess0c_drs,dess0c_dn 78 real*8 eud0c,deud0c_drs 79 real*8 eudc,eud1c,fud1c,dfud1c_dn 80 real*8 dfudc_dn,dfudc_dagr,dfudc_dtau 81 real*8 rs,drs_dn,dummy 82 real*8 D,dD_dx,dD_dz,dD_dn,dD_dagr,dD_dtau 83 real*8 Dt,dD_dDt,dDt_dx,dDt_dz 84 real*8 ec,fnc,fdnc,fdtauc 85* ***** constants ***** 86 real*8 pi,thrd,twthrd,frthrd,fvthrd,etthrd 87 parameter (pi = 3.14159265358979311599d0) 88 parameter (thrd = 1.0d0/3.0d0) 89 parameter (twthrd = 2.0d0/3.0d0) 90 parameter (frthrd = 4.0d0/3.0d0) 91 parameter (fvthrd = 5.0d0/3.0d0) 92 parameter (etthrd = 8.0d0/3.0d0) 93* ***** density cutoff parametersi ***** 94 real*8 tol,ETA,thresd 95 parameter (tol = 1.0d-10) 96 parameter (ETA = 1.0d-20) 97 parameter (thresd = 1.0d-2) 98* ***** VS98 constants ***** 99 real*8 aax,bbx,ccx,ddx,eex,ffx,alphax 100 real*8 aass,bbss,ccss,ddss,eess,ffss,alphass 101 real*8 aaopp,bbopp,ccopp,ddopp,eeopp,ffopp,alphaopp 102 real*8 cf,clda 103 parameter (cf = 9.115599720d0) 104c cf = 3.0d0/5.0d0*(6.0d0*pi*pi)**(2.0d0/3.0d0)) 105 parameter (clda = 0.9305257363491d0) 106 parameter (aax = -9.800683d-1) 107 parameter (bbx = -3.556788d-3) 108 parameter (ccx = 6.250326d-3) 109 parameter (ddx = -2.354518d-5) 110 parameter (eex = -1.282732d-4) 111 parameter (ffx = 3.574822d-4) 112 parameter (alphax = 0.00186726d0) 113 parameter (aass = 3.270912d-1) 114 parameter (bbss = -3.228915d-2) 115 parameter (ccss = -2.942406d-2) 116 parameter (ddss = 2.134222d-3) 117 parameter (eess = -5.451559d-3) 118 parameter (ffss = 1.577575d-2) 119 parameter (alphass = 0.00515088d0) 120 parameter (aaopp = 7.035010d-1) 121 parameter (bbopp = 7.694574d-3) 122 parameter (ccopp = 5.152765d-2) 123 parameter (ddopp = 3.394308d-5) 124 parameter (eeopp = -1.269420d-3) 125 parameter (ffopp = 1.296118d-3) 126 parameter (alphaopp = 0.00304966d0) 127 128 do i=1,n2ft3d 129 n = rho_in(i) + ETA 130 agr = agr_in(i) + ETA 131 tau = 2.0d0*tau_in(i) + ETA 132 133 n = 0.50d0*n 134 agr = 0.50d0*agr 135 tau = 0.50d0*tau 136 137 agr2 = agr*agr 138 inv_n = 1.0d0/n 139 n_13 = n**thrd 140 n_43 = n_13*n 141 n_53 = n_43*n_13 142 n_83 = n_53*n 143 144 x = agr2/n_83 145 dx_dn = -etthrd*x*inv_n 146 dx_dagr = 2.0d0*agr/n_83 147 z = tau/n_53 - cf 148 dz_dn = -fvthrd*tau/n_83 149 dz_dtau = 1.0d0/n_53 150 151* ***** VS98 Exchange ***** 152 gamma = 1.0d0 + alphax*(x + z) 153 xr = x/gamma 154 zr = z/gamma 155 156 call nwpw_GVT4(aax,bbx,ccx,ddx,eex,ffx,alphax,alphax, 157 > xr,zr,gamma,GG,dGdx,dGdz) 158 159 dG_dn = dGdx*dx_dn + dGdz*dz_dn 160 dG_dagr = dGdx*dx_dagr 161 dG_dtau = dGdz*dz_dtau 162 163 ex = clda*n_13*GG 164 fnx = frthrd*ex + clda*n_43*dG_dn 165 fdnx = clda*n_43*dG_dagr 166 fdtaux = clda*n_43*dG_dtau 167 168* ***** VS98 Correlation ***** 169* ***** Same-Spin (alpha-alpha/beta-beta) ***** 170 rs = (0.75d0/(pi*n))**thrd 171 drs_dn = -thrd*rs/n 172 173 call gen_PW91_c_rz(tol,rs,1.0d0,ess0c,dess0c_drs,dummy) 174 dess0c_dn = dess0c_drs*drs_dn 175 176 Dt = 1.0d0 - 0.25d0*x/(z + cf) 177 D = 0.0d0 178 dD_dn = 0.0d0 179 dD_dagr = 0.0d0 180 dD_dtau = 0.0d0 181 182 if (Dt .le. 0.0d0) then 183 D = 0.0d0 184 dD_dn = 0.0d0 185 dD_dagr = 0.0d0 186 dD_dtau = 0.0d0 187 else if (Dt .gt. 0.0d0 .and. Dt .lt. thresd) then 188 D = 2.0d0*Dt*Dt/thresd - Dt*Dt*Dt/(thresd*thresd) 189 dD_dDt = 4.0d0*Dt/thresd - 3.0d0*Dt*Dt/(thresd*thresd) 190 dDt_dx = -0.25d0/(z + cf) 191 dDt_dz = 0.25d0*x/((z + cf)*(z + cf)) 192 dD_dn = dD_dDt*(dDt_dx*dx_dn + dDt_dz*dz_dn) 193 dD_dagr = dD_dDt*dDt_dx*dx_dagr 194 dD_dtau = dD_dDt*dDt_dz*dz_dtau 195 else if(Dt .ge. thresd) then 196 D = Dt 197 dD_dx = -0.25d0/(z + cf) 198 dD_dz = 0.25d0*x/((z + cf)*(z + cf)) 199 dD_dn = dD_dx*dx_dn + dD_dz*dz_dn 200 dD_dagr = dD_dx*dx_dagr 201 dD_dtau = dD_dz*dz_dtau 202 end if 203 204 gamma = 1.0d0 + alphass*(x + z) 205 xr = x/gamma 206 zr = z/gamma 207 208 call nwpw_GVT4(aass,bbss,ccss,ddss,eess,ffss,alphass,alphass, 209 > xr,zr,gamma,GG,dGdx,dGdz) 210 211 dG_dn = dGdx*dx_dn + dGdz*dz_dn 212 dG_dagr = dGdx*dx_dagr 213 dG_dtau = dGdz*dz_dtau 214 215 ec = ess0c*GG*D 216 fnc = ec + n*dess0c_drs*drs_dn*GG*D 217 & + n*ess0c*dG_dn*D 218 & + n*ess0c*GG*dD_dn 219 fdnc = n*ess0c*(dG_dagr*D + GG*dD_dagr) 220 fdtauc = n*ess0c*(dG_dtau*D + GG*dD_dtau) 221 222* ***** Opposite-Spin (alpha-beta) ***** 223 n = 2.0d0*n 224 225 rs = (0.75d0/(pi*n))**thrd 226 drs_dn = -thrd*rs/n 227 228 call gen_PW91_c_rz(tol,rs,0.0d0,eud0c,deud0c_drs,dummy) 229 eud1c = eud0c - ess0c 230 fud1c = n*eud1c 231 dfud1c_dn = eud0c + n*deud0c_drs*drs_dn 232 & - ess0c - 0.50d0*n*dess0c_dn 233 234 x = 2.0d0*x 235 z = 2.0d0*z 236 237 gamma = 1.0d0 + alphaopp*(x + z) 238 xr = x/gamma 239 zr = z/gamma 240 241 call nwpw_GVT4(aaopp,bbopp,ccopp,ddopp,eeopp,ffopp, 242 > alphaopp,alphaopp, 243 > xr,zr,gamma,GG,dGdx,dGdz) 244 245 dG_dn = dGdx*dx_dn + dGdz*dz_dn 246 dG_dagr = dGdx*dx_dagr 247 dG_dtau = dGdz*dz_dtau 248 249 eudc = eud1c*GG 250 dfudc_dn = fud1c*dG_dn + GG*dfud1c_dn 251 dfudc_dagr = fud1c*dG_dagr 252 dfudc_dtau = fud1c*dG_dtau 253 254 ec = ec + eudc 255 fnc = fnc + dfudc_dn 256 fdnc = fdnc + dfudc_dagr 257 fdtauc = fdtauc + dfudc_dtau 258 259 xce(i) = x_parameter*ex + c_parameter*ec 260 fn(i) = x_parameter*fnx + c_parameter*fnc 261 fdn(i) = x_parameter*fdnx + c_parameter*fdnc 262 fdtau(i) = x_parameter*fdtaux + c_parameter*fdtauc 263 end do 264 265 return 266 end 267* ************************************************ 268* * * 269* * gen_VS98_unrestricted * 270* * * 271* ************************************************ 272 273* This function returns the VS98 exchange-correlation 274* energy density, xce, and its derivatives with respect 275* to nup, ndn, |grad nup|, |grad ndn|, tauup, taudn. 276 277* 278* Entry - n2ft3d : number of grid points 279* rho_in(*,2) : density (nup and ndn) 280* agr_in(*,3): |grad rho_in| (nup, ndn and n) 281* tau_in(*,2): tau (nup and ndn) 282* x_parameter: scale parameter for exchange 283* c_parameter: scale parameter for correlation 284* 285* Exit - xce(n2ft3d) : VS98 exchange correlation energy density 286* fn(n2ft3d,2) : d(n*xce)/dnup, d(n*xce)/dndn 287* fdn(n2ft3d,3) : d(n*xce)/d|grad nup|, d(n*xce)/d|grad ndn|, d(n*xce)/d|grad n| 288* fdtau(n2ft3d,2) : d(n*xce)/dtauup, d(n*xce)/dtaudn 289* 290 subroutine gen_VS98_unrestricted(n2ft3d,rho_in,agr_in,tau_in, 291 > x_parameter,c_parameter, 292 > xce,fn,fdn,fdtau) 293 implicit none 294* ***** input ***** 295 integer n2ft3d 296 real*8 rho_in(n2ft3d,2),agr_in(n2ft3d,3),tau_in(n2ft3d,2) 297 real*8 x_parameter,c_parameter 298* ***** output ***** 299 real*8 xce(n2ft3d),fn(n2ft3d,2),fdn(n2ft3d,3),fdtau(n2ft3d,2) 300* ***** local declarations ***** 301 integer i 302 real*8 nup,agrup,tauup 303 real*8 ndn,agrdn,taudn 304 real*8 n,agr,tau 305 real*8 x,dxdn,dxdagr,z,dzdn,dzdtau 306 real*8 xr,zr,gamma 307 real*8 inv_nup,nup_13,nup_43,nup_53,nup_83,agrup2 308 real*8 inv_ndn,ndn_13,ndn_43,ndn_53,ndn_83,agrdn2 309 real*8 xup,dxup_dnup,dxup_dagrup,zup,dzup_dnup,dzup_dtauup 310 real*8 xdn,dxdn_dndn,dxdn_dagrdn,zdn,dzdn_dndn,dzdn_dtaudn 311 real*8 GG,dGdx,dGdz 312 real*8 eupx,fnupx,fdnupx,fdtauupx 313 real*8 ednx,fndnx,fdndnx,fdtaudnx 314 real*8 euu0c,deuu0c_dnup 315 real*8 euuc,dfuuc_dnup,dfuuc_dagrup,dfuuc_dtauup 316 real*8 edd0c,dedd0c_dndn 317 real*8 eddc,dfddc_dndn,dfddc_dagrdn,dfddc_dtaudn 318 real*8 eud0c,deud0c_drs,deud0c_dzeta 319 real*8 eudc,eud1c,fud1c,dfud1c_dnup,dfud1c_dndn 320 real*8 dfudc_dnup,dfudc_dagrup,dfudc_dtauup 321 real*8 dfudc_dndn,dfudc_dagrdn,dfudc_dtaudn 322 real*8 rs,drs_dn,zeta,dzeta_dnup,dzeta_dndn,deuu0c_drs,dedd0c_drs 323 real*8 D,dD_dx,dD_dz,dD_dn,dD_dagr,dD_dtau,dummy 324 real*8 Dt,dD_dDt,dDt_dx,dDt_dz 325 real*8 eupc,fnupc,fdnupc,fdtauupc 326 real*8 ednc,fndnc,fdndnc,fdtaudnc 327 real*8 ex,ec 328* ***** constants ***** 329 real*8 pi,thrd,twthrd,frthrd,fvthrd,etthrd 330 parameter (pi = 3.14159265358979311599d0) 331 parameter (thrd = 1.0d0/3.0d0) 332 parameter (twthrd = 2.0d0/3.0d0) 333 parameter (frthrd = 4.0d0/3.0d0) 334 parameter (fvthrd = 5.0d0/3.0d0) 335 parameter (etthrd = 8.0d0/3.0d0) 336* ***** density cutoff parameters ***** 337 real*8 tol,ETA,thresd 338 parameter (tol = 1.0d-10) 339 parameter (ETA = 1.0d-20) 340 parameter (thresd = 1.0d-2) 341* ***** VS98 constants ***** 342 real*8 aax,bbx,ccx,ddx,eex,ffx,alphax 343 real*8 aass,bbss,ccss,ddss,eess,ffss,alphass 344 real*8 aaopp,bbopp,ccopp,ddopp,eeopp,ffopp,alphaopp 345 real*8 cf,clda 346 parameter (cf = 9.115599720d0) 347c cf = 3.0d0/5.0d0*(6.0d0*pi*pi)**(2.0d0/3.0d0)) 348 parameter (clda = 0.9305257363491d0) 349 parameter (aax = -9.800683d-1) 350 parameter (bbx = -3.556788d-3) 351 parameter (ccx = 6.250326d-3) 352 parameter (ddx = -2.354518d-5) 353 parameter (eex = -1.282732d-4) 354 parameter (ffx = 3.574822d-4) 355 parameter (alphax = 0.00186726d0) 356 parameter (aass = 3.270912d-1) 357 parameter (bbss = -3.228915d-2) 358 parameter (ccss = -2.942406d-2) 359 parameter (ddss = 2.134222d-3) 360 parameter (eess = -5.451559d-3) 361 parameter (ffss = 1.577575d-2) 362 parameter (alphass = 0.00515088d0) 363 parameter (aaopp = 7.035010d-1) 364 parameter (bbopp = 7.694574d-3) 365 parameter (ccopp = 5.152765d-2) 366 parameter (ddopp = 3.394308d-5) 367 parameter (eeopp = -1.269420d-3) 368 parameter (ffopp = 1.296118d-3) 369 parameter (alphaopp = 0.00304966d0) 370 371 do i=1,n2ft3d 372 nup = rho_in(i,1) + ETA 373 agrup = agr_in(i,1) + ETA 374 tauup = tau_in(i,1) + ETA 375 ndn = rho_in(i,2) + ETA 376 agrdn = agr_in(i,2) + ETA 377 taudn = tau_in(i,2) + ETA 378 379 n = nup + ndn 380 381* ***** VS98 Exchange ***** 382* ***** Up ***** 383 agrup2 = agrup*agrup 384 inv_nup = 1.0d0/nup 385 nup_13 = nup**thrd 386 nup_43 = nup_13*nup 387 nup_53 = nup_43*nup_13 388 nup_83 = nup_53*nup 389 390 xup = agrup2/nup_83 391 dxup_dnup = -etthrd*xup*inv_nup 392 dxup_dagrup = 2.0d0*agrup/nup_83 393 zup = tauup/nup_53 - cf 394 dzup_dnup = -fvthrd*tauup/nup_83 395 dzup_dtauup = 1.0d0/nup_53 396 397 gamma = 1.0d0 + alphax*(xup + zup) 398 xr = xup/gamma 399 zr = zup/gamma 400 401 call nwpw_GVT4(aax,bbx,ccx,ddx,eex,ffx,alphax,alphax, 402 > xr,zr,gamma,GG,dGdx,dGdz) 403 404 eupx = clda*nup_13*GG 405 fnupx = frthrd*eupx + clda*nup_43*(dGdx*dxup_dnup 406 & + dGdz*dzup_dnup) 407 fdnupx = clda*nup_43*(dGdx*dxup_dagrup) 408 fdtauupx = clda*nup_43*(dGdz*dzup_dtauup) 409 410* ***** Down ***** 411 agrdn2 = agrdn*agrdn 412 inv_ndn = 1.0d0/ndn 413 ndn_13 = ndn**thrd 414 ndn_43 = ndn_13*ndn 415 ndn_53 = ndn_43*ndn_13 416 ndn_83 = ndn_53*ndn 417 418 xdn = agrdn2/ndn_83 419 dxdn_dndn = -etthrd*xdn*inv_ndn 420 dxdn_dagrdn = 2.0d0*agrdn/ndn_83 421 zdn = taudn/ndn_53 - cf 422 dzdn_dndn = -fvthrd*taudn/ndn_83 423 dzdn_dtaudn = 1.0d0/ndn_53 424 425 gamma = 1.0d0 + alphax*(xdn + zdn) 426 xr = xdn/gamma 427 zr = zdn/gamma 428 429 call nwpw_GVT4(aax,bbx,ccx,ddx,eex,ffx,alphax,alphax, 430 > xr,zr,gamma,GG,dGdx,dGdz) 431 432 ednx = clda*ndn_13*GG 433 fndnx = frthrd*ednx + clda*ndn_43*(dGdx*dxdn_dndn 434 & + dGdz*dzdn_dndn) 435 fdndnx = clda*ndn_43*(dGdx*dxdn_dagrdn) 436 fdtaudnx = clda*ndn_43*(dGdz*dzdn_dtaudn) 437 438 ex = (eupx*nup + ednx*ndn)/n 439 440* ***** VS98 Correlation ***** 441* ***** Same-Spin ***** 442* ***** alpha-alpha ***** 443 rs = (0.75d0/(pi*nup))**thrd 444 drs_dn = -thrd*rs/nup 445 446 call gen_PW91_c_rz(tol,rs,1.0d0,euu0c,deuu0c_drs,dummy) 447 deuu0c_dnup = deuu0c_drs*drs_dn 448 449 Dt = 1.0d0 - 0.25d0*xup/(zup + cf) 450 D = 0.0d0 451 dD_dn = 0.0d0 452 dD_dagr = 0.0d0 453 dD_dtau = 0.0d0 454 455 if (Dt .le. 0.0d0) then 456 D = 0.0d0 457 dD_dn = 0.0d0 458 dD_dagr = 0.0d0 459 dD_dtau = 0.0d0 460 else if (Dt .gt. 0.0d0 .and. Dt .lt. thresd) then 461 D = 2.0d0*Dt*Dt/thresd - Dt*Dt*Dt/(thresd*thresd) 462 dD_dDt = 4.0d0*Dt/thresd - 3.0d0*Dt*Dt/(thresd*thresd) 463 dDt_dx = -0.25d0/(zup + cf) 464 dDt_dz = 0.25d0*xup/((zup + cf)*(zup + cf)) 465 dD_dn = dD_dDt*(dDt_dx*dxup_dnup + dDt_dz*dzup_dnup) 466 dD_dagr = dD_dDt*dDt_dx*dxup_dagrup 467 dD_dtau = dD_dDt*dDt_dz*dzup_dtauup 468 else if(Dt .ge. thresd) then 469 D = Dt 470 dD_dx = -0.25d0/(zup + cf) 471 dD_dz = 0.25d0*xup/((zup + cf)*(zup + cf)) 472 dD_dn = dD_dx*dxup_dnup + dD_dz*dzup_dnup 473 dD_dagr = dD_dx*dxup_dagrup 474 dD_dtau = dD_dz*dzup_dtauup 475 end if 476 477 gamma = 1.0d0 + alphass*(xup + zup) 478 xr = xup/gamma 479 zr = zup/gamma 480 481 call nwpw_GVT4(aass,bbss,ccss,ddss,eess,ffss,alphass,alphass, 482 > xr,zr,gamma,GG,dGdx,dGdz) 483 484 eupc = euu0c*GG*D 485 fnupc = eupc + nup*deuu0c_drs*drs_dn*GG*D 486 & + nup*euu0c*(dGdx*dxup_dnup+dGdz*dzup_dnup)*D 487 & + nup*euu0c*GG*dD_dn 488 fdnupc = nup*euu0c*(dGdx*dxup_dagrup*D + GG*dD_dagr) 489 fdtauupc = nup*euu0c*(dGdz*dzup_dtauup*D + GG*dD_dtau) 490 491* ***** beta-beta ***** 492 rs = (0.75d0/(pi*ndn))**thrd 493 drs_dn = -thrd*rs/ndn 494 495 call gen_PW91_c_rz(tol,rs,1.0d0,edd0c,dedd0c_drs,dummy) 496 dedd0c_dndn = dedd0c_drs*drs_dn 497 498 Dt = 1.0d0 - 0.25d0*xdn/(zdn + cf) 499 D = 0.0d0 500 dD_dn = 0.0d0 501 dD_dagr = 0.0d0 502 dD_dtau = 0.0d0 503 504 if (Dt .le. 0.0d0) then 505 D = 0.0d0 506 dD_dn = 0.0d0 507 dD_dagr = 0.0d0 508 dD_dtau = 0.0d0 509 else if (Dt .gt. 0.0d0 .and. Dt .lt. thresd) then 510 D = 2.0d0*Dt*Dt/thresd - Dt*Dt*Dt/(thresd*thresd) 511 dD_dDt = 4.0d0*Dt/thresd - 3.0d0*Dt*Dt/(thresd*thresd) 512 dDt_dx = -0.25d0/(zdn + cf) 513 dDt_dz = 0.25d0*xdn/((zdn + cf)*(zdn + cf)) 514 dD_dn = dD_dDt*(dDt_dx*dxdn_dndn + dDt_dz*dzdn_dndn) 515 dD_dagr = dD_dDt*dDt_dx*dxdn_dagrdn 516 dD_dtau = dD_dDt*dDt_dz*dzdn_dtaudn 517 else if(Dt .ge. thresd) then 518 D = Dt 519 dD_dx = -0.25d0/(zdn + cf) 520 dD_dz = 0.25d0*xdn/((zdn + cf)*(zdn + cf)) 521 dD_dn = dD_dx*dxdn_dndn + dD_dz*dzdn_dndn 522 dD_dagr = dD_dx*dxdn_dagrdn 523 dD_dtau = dD_dz*dzdn_dtaudn 524 end if 525 526 gamma = 1.0d0 + alphass*(xdn + zdn) 527 xr = xdn/gamma 528 zr = zdn/gamma 529 530 call nwpw_GVT4(aass,bbss,ccss,ddss,eess,ffss,alphass,alphass, 531 > xr,zr,gamma,GG,dGdx,dGdz) 532 533 ednc = edd0c*GG*D 534 fndnc = ednc + ndn*dedd0c_drs*drs_dn*GG*D 535 & + ndn*edd0c*(dGdx*dxdn_dndn+dGdz*dzdn_dndn)*D 536 & + ndn*edd0c*GG*dD_dn 537 fdndnc = ndn*edd0c*(dGdx*dxdn_dagrdn*D + GG*dD_dagr) 538 fdtaudnc = ndn*edd0c*(dGdz*dzdn_dtaudn*D + GG*dD_dtau) 539 540 ec = (eupc*nup + ednc*ndn)/n 541 542* ***** Opposite-Spin (alpha-beta) ***** 543 rs = (0.75d0/(pi*n))**thrd 544 drs_dn = -thrd*rs/n 545 zeta = (nup - ndn)/n 546 dzeta_dnup = ( 1.0d0 - zeta)/n 547 dzeta_dndn = (-1.0d0 - zeta)/n 548 549 call gen_PW91_c_rz(tol,rs,zeta,eud0c,deud0c_drs,deud0c_dzeta) 550 eud1c = eud0c - (euu0c*nup + edd0c*ndn)/n 551 fud1c = n*eud1c 552 dfud1c_dnup = eud0c + n*(deud0c_drs*drs_dn 553 & + deud0c_dzeta*dzeta_dnup) - euu0c 554 & - nup*deuu0c_dnup 555 dfud1c_dndn = eud0c + n*(deud0c_drs*drs_dn 556 & + deud0c_dzeta*dzeta_dndn) - edd0c 557 & - ndn*dedd0c_dndn 558 559 x = xup + xdn 560 z = zup + zdn 561 562 gamma = 1.0d0 + alphaopp*(x + z) 563 xr = x/gamma 564 zr = z/gamma 565 566 call nwpw_GVT4(aaopp,bbopp,ccopp,ddopp,eeopp,ffopp, 567 > alphaopp,alphaopp, 568 > xr,zr,gamma,GG,dGdx,dGdz) 569 570 eudc = eud1c*GG 571 dfudc_dnup = dfud1c_dnup*GG 572 & + fud1c*(dGdx*dxup_dnup + dGdz*dzup_dnup) 573 dfudc_dndn = dfud1c_dndn*GG + 574 & fud1c*(dGdx*dxdn_dndn + dGdz*dzdn_dndn) 575 dfudc_dagrup = fud1c*dGdx*dxup_dagrup 576 dfudc_dagrdn = fud1c*dGdx*dxdn_dagrdn 577 dfudc_dtauup = fud1c*dGdz*dzup_dtauup 578 dfudc_dtaudn = fud1c*dGdz*dzdn_dtaudn 579 580 ec = ec + eudc 581 fnupc = fnupc + dfudc_dnup 582 fndnc = fndnc + dfudc_dndn 583 fdnupc = fdnupc + dfudc_dagrup 584 fdndnc = fdndnc + dfudc_dagrdn 585 fdtauupc = fdtauupc + dfudc_dtauup 586 fdtaudnc = fdtaudnc + dfudc_dtaudn 587 588 xce(i) = x_parameter*ex + c_parameter*ec 589 fn(i,1) = x_parameter*fnupx + c_parameter*fnupc 590 fn(i,2) = x_parameter*fndnx + c_parameter*fndnc 591 fdn(i,1) = x_parameter*fdnupx + c_parameter*fdnupc 592 fdn(i,2) = x_parameter*fdndnx + c_parameter*fdndnc 593 fdn(i,3) = 0.0d0 594 fdtau(i,1) = x_parameter*fdtauupx + c_parameter*fdtauupc 595 fdtau(i,2) = x_parameter*fdtaudnx + c_parameter*fdtaudnc 596 end do 597 598 return 599 end 600 601