1* 2* $Id$ 3* 4 5* ************************************ 6* * * 7* * gen_revPBE_unrestricted * 8* * * 9* ************************************ 10* 11* This function returns the revPBE exchange-correlation 12* energy density, xce, and its derivatives with respect 13* to nup, ndn, |grad nup|, |grad ndn|, and |grad n|. 14* 15* Entry - n2ft3d : number of grid points 16* dn_in(*,2) : spin densites nup and ndn 17* agr_in(*,3): |grad nup|, |grad ndn|, and |grad n| 18* x_parameter: scale parameter for exchange 19* c_parameter: scale parameter for correlation 20* 21* Exit - xce(*) : revPBE energy density 22* - fn(*,2) : d(n*xce)/dnup, d(n*xce)/dndn 23* - fdn(*,3): d(n*xce)/d|grad nup|, d(n*xce)/d|grad ndn| 24* d(n*xce)/d|grad n| 25 26 subroutine gen_revPBE_BW_unrestricted(n2ft3d, 27 > dn_in,agr_in, 28 > x_parameter,c_parameter, 29 > xce,fn,fdn) 30 implicit none 31 32 integer n2ft3d 33 real*8 dn_in(n2ft3d,2) 34 real*8 agr_in(n2ft3d,3) 35 real*8 x_parameter,c_parameter 36 real*8 xce(n2ft3d) 37 real*8 fn(n2ft3d,2) 38 real*8 fdn(n2ft3d,3) 39 40* **** Density cutoff parameter **** 41 real*8 DNS_CUT,ETA,ETA2,alpha_zeta,alpha_zeta2 42 parameter (DNS_CUT = 1.0d-20) 43 parameter (ETA=1.0d-20) 44 parameter (ETA2=1.0d-14) 45 parameter (alpha_zeta=(1.0d0-ETA2)) 46 parameter (alpha_zeta2=(1.0d0-ETA2)) 47 48c ***** revPBE GGA exchange constants ****** 49 real*8 MU,KAPPA 50 parameter (MU = 0.2195149727645171d0) 51 !parameter (KAPPA = 0.8040000000000000d0) 52 parameter (KAPPA = 1.2450000000000000d0) 53 54c ****** PBE96 GGA correlation constants ****** 55 real*8 GAMMA,BETA,BOG 56 parameter (GAMMA = 0.031090690869655d0) 57 parameter (BETA = 0.066724550603149d0) 58 !parameter (BETA = 0.066725d0) 59 parameter (BOG = BETA/GAMMA) 60 61 62c ****** Perdew-Wang92 LDA correlation coefficients ******* 63 real*8 GAM,iGAM,FZZ,iFZZ 64 parameter (GAM = 0.519842099789746329d0) 65 parameter (iGAM = 1.0d0/GAM) 66 parameter (FZZ = (8.0d0/(9.0d0*GAM)) ) 67 parameter (iFZZ = 0.125d0*9.0d0*GAM) 68 69 real*8 A_1,A1_1,B1_1,B2_1,B3_1,B4_1 70 parameter (A_1 = 0.0310907d0) 71 !parameter (A_1 = 0.031091d0) 72 parameter (A1_1 = 0.2137000d0) 73 parameter (B1_1 = 7.5957000d0) 74 parameter (B2_1 = 3.5876000d0) 75 parameter (B3_1 = 1.6382000d0) 76 parameter (B4_1 = 0.4929400d0) 77 78 real*8 A_2,A1_2,B1_2,B2_2,B3_2,B4_2 79 parameter (A_2 = 0.01554535d0) 80 !parameter (A_2 = 0.015545d0) 81 parameter (A1_2 = 0.20548000d0) 82 parameter (B1_2 = 14.11890000d0) 83 parameter (B2_2 = 6.19770000d0) 84 parameter (B3_2 = 3.36620000d0) 85 parameter (B4_2 = 0.62517000d0) 86 87 real*8 A_3,A1_3,B1_3,B2_3,B3_3,B4_3 88 parameter (A_3 = 0.0168869d0) 89 !parameter (A_3 = 0.016887d0) 90 parameter (A1_3 = 0.1112500d0) 91 parameter (B1_3 = 10.3570000d0) 92 parameter (B2_3 = 3.6231000d0) 93 parameter (B3_3 = 0.8802600d0) 94 parameter (B4_3 = 0.4967100d0) 95 96c **** other constants **** 97 real*8 onethird,fourthird,fivethird,onesixthm 98 real*8 twothird,sevensixthm 99 real*8 onethirdm 100 parameter (onethird=1.0d0/3.0d0) 101 parameter (onethirdm=-1.0d0/3.0d0) 102 parameter (twothird=2.0d0/3.0d0) 103 parameter (fourthird=4.0d0/3.0d0) 104 parameter (fivethird=5.0d0/3.0d0) 105 parameter (onesixthm=-1.0d0/6.0d0) 106 parameter (sevensixthm=-7.0d0/6.0d0) 107 108c **** local variables **** 109 integer i 110 real*8 n,agr 111 real*8 nup,agrup 112 real*8 ndn,agrdn 113 real*8 kf,ks,s,P0,n_onethird,pi,rs_scale 114 real*8 rs ! Wigner radius 115 real*8 rss ! rss = sqrt(rs) 116 real*8 rs_n ! rs_n = n*drs/dn 117 real*8 t,t2,t4,t6 118 real*8 t_nup ! t_nup = n*dt/dnup 119 real*8 t_ndn ! t_ndn = n*dt/dndn 120 real*8 t_agr ! t_agr = n*dt/dagr 121 real*8 zet,twoksg 122 real*8 zet_nup ! zet_nup = n*dzet/dnup 123 real*8 zet_ndn ! zet_nup = n*dzet/dnup 124 real*8 zetp_1_3,zetm_1_3 125 real*8 zetpm_1_3,zetmm_1_3 126 real*8 phi,phi3,phi4 127 real*8 phi_zet 128 real*8 A,A2 129 real*8 A_phi,A_ec_lda 130 real*8 Q0,Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8 131 real*8 PON,FZ,z4 132 real*8 tau 133 real*8 F 134 real*8 Fs ! dF/ds 135 real*8 Hpbe 136 real*8 Hpbe_t ! dHpbe/dt 137 real*8 Hpbe_phi ! dHpbe/dphi 138 real*8 Hpbe_ec_lda ! dHpbe/d(ec_lda) 139 real*8 Hpbe_nup,Hpbe_ndn ! n*dHpbe/dnup, n*dHpbe/dndn 140 real*8 Ipbe 141 real*8 Ipbe_t,Ipbe_A ! dIpbe/dt, dIpbe/dA 142 143 real*8 exup,exdn,ex,ex_lda 144 real*8 ecu,ecp,eca,ec,ec_lda 145 real*8 ecu_rs,ecp_rs,eca_rs 146 real*8 ec_lda_rs,ec_lda_zet ! d(ec_lda)/drs, d(ec_lda)/dzet 147 real*8 ec_lda_nup,ec_lda_ndn ! n*d(ec_lda)/dnup, n*d(ec_lda)/dndn 148 real*8 fnxup,fdnxup ! d(n*ex)/dnup, d(n*ex)/dndn 149 real*8 fnxdn,fdnxdn ! d(n*ex)/d|grad nup|, d(n*ex)/d|grad ndn| 150 real*8 fncup,fncdn ! d(n*ec)/dnup, d(n*ec)/dndn 151 real*8 fdnx_const 152 153 pi = 4.0d0*datan(1.0d0) 154 rs_scale = (0.75d0/pi)**onethird 155 fdnx_const = -3.0d0/(8.0d0*pi) 156 157 158!$OMP DO 159 do i=1,n2ft3d 160 nup = dn_in(i,1)+ETA 161 agrup = agr_in(i,1) 162 163 ndn = dn_in(i,2)+ETA 164 agrdn = agr_in(i,2) 165 166c **************************************************************** 167c ***** calculate polarized Exchange energies and potentials ***** 168c **************************************************************** 169 170c ************ 171c **** up **** 172c ************ 173 n = 2.0d0*nup 174 agr = 2.0d0*agrup 175 176 n_onethird = (3.0d0*n/pi)**onethird 177 ex_lda = -0.75d0*n_onethird 178 179 kf = (3.0d0*pi*pi*n)**onethird 180 s = agr/(2.0d0*kf*n) 181 P0 = 1.0d0 + (MU/KAPPA)*s*s 182 183 F = (1.0d0 + KAPPA - KAPPA/P0) 184 Fs = 2.0d0*MU/(P0*P0)*s 185 186 exup = ex_lda*F 187 fnxup = fourthird*(exup - ex_lda*Fs*s) 188 fdnxup = fdnx_const*Fs 189 190c ************** 191c **** down **** 192c ************** 193 n = 2.0d0*ndn 194 agr = 2.0d0*agrdn 195 196 n_onethird = (3.0d0*n/pi)**onethird 197 ex_lda = -0.75d0*n_onethird 198 199 kf = (3.0d0*pi*pi*n)**onethird 200 s = agr/(2.0d0*kf*n) 201 P0 = 1.0d0 + (MU/KAPPA)*s*s 202 203 F = (1.0d0 + KAPPA - KAPPA/P0) 204 Fs = 2.0d0*MU/(P0*P0)*s 205 206 exdn = ex_lda*F 207 fnxdn = fourthird*(exdn - ex_lda*Fs*s) 208 fdnxdn = fdnx_const*Fs 209 210 n = nup+ndn 211 212 ex = (exup*nup+ exdn*ndn)/ n 213 214c ******************************************************************* 215c ***** calculate polarized correlation energies and potentials ***** 216c ******************************************************************* 217 agr = agr_in(i,3) 218 219 zet = (nup-ndn)/n 220c if (zet.gt.0.0d0) zet = zet - ETA2 221c if (zet.lt.0.0d0) zet = zet + ETA2 222c if (dabs(dn_in(i,2)).gt.DNS_CUT) zet_nup = 2*ndn/n**2 223c if (dabs(dn_in(i,1)).gt.DNS_CUT) zet_ndn = -2*nup/n**2 224c if (dabs(dn_in(i,2)).gt.DNS_CUT) zet_nup = 2*ndn/n 225c zet_nup = 2*ndn/n 226c zet_ndn = -2*nup/n 227 zet_nup = -(zet - 1.0d0) 228 zet_ndn = -(zet + 1.0d0) 229 zetpm_1_3 = (1.0d0+zet*alpha_zeta)**onethirdm 230 zetmm_1_3 = (1.0d0-zet*alpha_zeta)**onethirdm 231 zetp_1_3 = (1.0d0+zet*alpha_zeta)*zetpm_1_3**2 232 zetm_1_3 = (1.0d0-zet*alpha_zeta)*zetmm_1_3**2 233 234 235 phi = 0.5d0*( zetp_1_3**2 + zetm_1_3**2) 236 phi_zet = alpha_zeta*( zetpm_1_3 - zetmm_1_3)/3.0d0 237 F =( (1.0d0+zet*alpha_zeta)*zetp_1_3 238 > + (1.0d0-zet*alpha_zeta)*zetm_1_3 239 > - 2.0d0)*iGAM 240 241 FZ = (zetp_1_3 - zetm_1_3)*(alpha_zeta*fourthird*iGAM) 242 243 244 245* **** calculate Wigner radius **** 246 rs = rs_scale/(n**onethird) 247 rss = dsqrt(rs) 248 249* **** calculate n*drs/dn **** 250c rs_n = onethirdm*rs/n 251 rs_n = onethirdm*rs 252 253 254 255c **** calculate t **** 256 kf = (3.0d0*pi*pi*n)**onethird 257 ks = dsqrt(4.0d0*kf/pi) 258 259 twoksg = 2.0d0*ks*phi 260 261 t = agr/(twoksg*n) 262 263* *** calculate n*dt/dnup, n*dt/dndn, n*dt/d|grad n| **** 264 t_nup = sevensixthm*t - (phi_zet)*(zet_nup)*t/phi 265 t_ndn = sevensixthm*t - (phi_zet)*(zet_ndn)*t/phi 266 t_agr = 1.0d0/(twoksg) 267 268 269 270 271c ************************************************** 272c ***** compute LSDA correlation energy density **** 273c ************************************************** 274 call LSDT2(A_1,A1_1,B1_1,B2_1,B3_1,B4_1,rss,ecu,ecu_rs) 275 call LSDT2(A_2,A1_2,B1_2,B2_2,B3_2,B4_2,rss,ecp,ecp_rs) 276 call LSDT2(A_3,A1_3,B1_3,B2_3,B3_3,B4_3,rss,eca,eca_rs) 277 278 z4 = zet**4 279 280 ec_lda = ecu*(1.0d0-F*z4) 281 > + ecp*F*z4 282 > - eca*F*(1.0d0-z4)/FZZ 283 284 ec_lda_rs = ecu_rs*(1.0d0-F*z4) 285 > + ecp_rs*F*z4 286 > - eca_rs*F*(1.0d0-z4)/FZZ 287 288 ec_lda_zet = (4.0d0*(zet**3)*F + FZ*z4)*(ecp-ecu+eca*iFZZ) 289 > - FZ*eca*iFZZ 290 291 292 293 294c ******************************************** 295c **** calculate PBE96 correlation energy **** 296c ******************************************** 297 phi3 = phi**3 298 phi4 = phi3*phi 299 PON = -ec_lda/(phi3*GAMMA) 300 tau = DEXP(PON) 301 302 A = BOG/(tau-1.0d0+ETA) 303 A2 = A*A 304 t2 = t*t 305 t4 = t2*t2 306 t6 = t4*t2 307 Q4 = 1.0d0 + A*t2 308 Q5 = 1.0d0 + 2.0d0*A*t2 309 Q6 = 2.0d0 + A*t2 310 Q7 = 1.0d0+A*t2+A2*t4 311 Q8 = Q7*Q7 312 313 Ipbe = 1.0d0 + BOG*t2*Q4/Q7 314 Hpbe = GAMMA*phi3*DLOG(Ipbe) 315 316 Ipbe_t = BOG*(2.0d0*t)*Q5/Q8 317 Ipbe_A = -BOG*(A*t6) *Q6/Q8 318 319 A_ec_lda = tau/(BETA*phi3)*A2 320 A_phi = -3.0d0*ec_lda*tau/(BETA*phi4)*A2 321 322 323 Hpbe_ec_lda = (GAMMA*phi3/Ipbe)*Ipbe_A*A_ec_lda 324 325 Hpbe_phi = 3.0d0*Hpbe/phi 326 > + (GAMMA*phi3/Ipbe)*Ipbe_A*A_phi 327 328 Hpbe_t = (GAMMA*phi3/Ipbe)*Ipbe_t 329 330 ec_lda_nup = ec_lda_zet 331 > - zet * ec_lda_zet 332 > + rs_n * ec_lda_rs 333 ec_lda_ndn = -ec_lda_zet 334 > - zet * ec_lda_zet 335 > + rs_n * ec_lda_rs 336 337 338 339 Hpbe_nup = ec_lda_nup * Hpbe_ec_lda 340 > + phi_zet*zet_nup * Hpbe_phi 341 > + t_nup * Hpbe_t 342 343 Hpbe_ndn = ec_lda_ndn * Hpbe_ec_lda 344 > + phi_zet*zet_ndn * Hpbe_phi 345 > + t_ndn * Hpbe_t 346 347 348 349 ec = ec_lda + Hpbe 350 351 fncup = ec + (ec_lda_nup + Hpbe_nup) 352 fncdn = ec + (ec_lda_ndn + Hpbe_ndn) 353 354 xce(i) = x_parameter*ex + c_parameter*ec 355 fn(i,1) = x_parameter*fnxup + c_parameter*fncup 356 fn(i,2) = x_parameter*fnxdn + c_parameter*fncdn 357 358 fdn(i,1) = x_parameter*fdnxup 359 fdn(i,2) = x_parameter*fdnxdn 360 fdn(i,3) = c_parameter*t_agr*Hpbe_t 361 362 end do 363!$OMP END DO 364 365 366 367 return 368 end 369 370 371 372 373 374* ************************************ 375* * * 376* * gen_revPBE_BW_restricted * 377* * * 378* ************************************ 379* 380* This routine calculates the revPBE exchange-correlation 381* potential(xcp) and energy density(xce). 382* 383* 384* Entry - n2ft3d : number of grid points 385* rho_in(*) : density (nup+ndn) 386* agr_in(*): |grad rho_in| 387* x_parameter: scale parameter for exchange 388* c_parameter: scale parameter for correlation 389* 390* Exit - xce(n2ft3d) : revPBE exchange correlation energy density 391* fn(n2ft3d) : d(n*xce)/dn 392* fdn(n2ft3d) : d(n*xce/d|grad n| 393* 394 subroutine gen_revPBE_BW_restricted(n2ft3d,rho_in,agr_in, 395 > x_parameter,c_parameter, 396 > xce,fn,fdn) 397 implicit none 398 399 integer n2ft3d 400 real*8 rho_in(n2ft3d) 401 real*8 agr_in(n2ft3d) 402 real*8 x_parameter,c_parameter 403 real*8 xce(n2ft3d) 404 real*8 fn(n2ft3d) 405 real*8 fdn(n2ft3d) 406 407 408* **** Density cutoff parameter **** 409 real*8 DNS_CUT,ETA 410 parameter (DNS_CUT = 1.0d-20) 411 parameter (ETA = 1.0d-20) 412 413c ***** revPBE GGA exchange constants ****** 414 real*8 MU,KAPPA 415 parameter (MU = 0.2195149727645171d0) 416 !parameter (KAPPA = 0.8040000000000000d0) 417 parameter (KAPPA = 1.2450000000000000d0) 418 419c ****** revPBE GGA correlation constants ****** 420 real*8 GAMMA,BETA,BOG 421 parameter (GAMMA = 0.031090690869655d0) 422 parameter (BETA = 0.066724550603149d0) 423 parameter (BOG = BETA/GAMMA) 424 425 426c ****** Perdew-Wang92 LDA correlation coefficients ******* 427 real*8 A_1,A1_1,B1_1,B2_1,B3_1,B4_1 428 parameter (A_1 = 0.0310907d0) 429 parameter (A1_1 = 0.2137000d0) 430 parameter (B1_1 = 7.5957000d0) 431 parameter (B2_1 = 3.5876000d0) 432 parameter (B3_1 = 1.6382000d0) 433 parameter (B4_1 = 0.4929400d0) 434 435 real*8 A_2,A1_2,B1_2,B2_2,B3_2,B4_2 436 parameter (A_2 = 0.01554535d0) 437 parameter (A1_2 = 0.20548000d0) 438 parameter (B1_2 = 14.11890000d0) 439 parameter (B2_2 = 6.19770000d0) 440 parameter (B3_2 = 3.36620000d0) 441 parameter (B4_2 = 0.62517000d0) 442 443 real*8 A_3,A1_3,B1_3,B2_3,B3_3,B4_3 444 parameter (A_3 = 0.0168869d0) 445 parameter (A1_3 = 0.1112500d0) 446 parameter (B1_3 = 10.3570000d0) 447 parameter (B2_3 = 3.6231000d0) 448 parameter (B3_3 = 0.8802600d0) 449 parameter (B4_3 = 0.4967100d0) 450 451c **** other constants **** 452 real*8 onethird,fourthird,sevensixths 453 parameter (onethird=1.0d0/3.0d0) 454 parameter (fourthird=4.0d0/3.0d0) 455 parameter (sevensixths=7.0d0/6.0d0) 456 457c **** local variables **** 458 integer i 459 real*8 n,agr 460 real*8 kf,ks,s,P0,n_onethird,pi,rs_scale 461 real*8 fdnx_const 462 real*8 rs,rss,t,t2,t4,t6 463 real*8 Q0,Q1,Q2,Q3,Q4,Q5,Q8,Q9,B 464 real*8 Ht 465 real*8 B_ec,Hrs,H_B 466 real*8 F,Fs 467 468 real*8 ex_lda,ec_lda 469 real*8 ec_lda_rs 470 real*8 ex,ec,H 471 real*8 fnx,fdnx,fnc,fdnc 472 473 474 pi = 4.0d0*datan(1.0d0) 475 rs_scale = (0.75d0/pi)**onethird 476 fdnx_const = -3.0d0/(8.0d0*pi) 477 478!$OMP DO 479 do i=1,n2ft3d 480 n = rho_in(i)+ETA 481 agr = agr_in(i) 482 483c ***** calculate unpolarized Exchange energies and potentials ***** 484 n_onethird = (3.0d0*n/pi)**onethird 485 ex_lda = -0.75d0*n_onethird 486 487 kf = (3.0d0*pi*pi*n)**onethird 488 s = agr/(2.0d0*kf*n) 489 P0 = 1.0d0 + (MU/KAPPA)*s*s 490 491c if (n.gt.DNS_CUT) then 492c F = (1.0d0 + KAPPA - KAPPA/P0) 493c Fs = 2.0d0*MU/(P0*P0)*s 494c else 495c F = 1.0d0 496c Fs = 0.0d0 497c end if 498 F = (1.0d0 + KAPPA - KAPPA/P0) 499 Fs = 2.0d0*MU/(P0*P0)*s 500 501 ex = ex_lda*F 502 fnx = fourthird*(ex - ex_lda*Fs*s) 503 fdnx = fdnx_const*Fs 504 505 506* ********************************************************************* 507c ***** calculate unpolarized correlation energies and potentials ***** 508* ********************************************************************* 509 510c **** calculate rs and t **** 511 rs = rs_scale/(n**onethird) 512 rss = dsqrt(rs) 513 514 kf = (3.0d0*pi*pi*n)**onethird 515 ks = dsqrt(4.0d0*kf/pi) 516 t = agr/(2.0d0*ks*n) 517 518 519c **** unpolarized LDA correlation energy **** 520c **** ec_p = correlation energy **** 521c **** ec_p_rs = dec_p/drs **** 522c **** uc_p = dec_p/dn **** 523 call LSDT2(A_1,A1_1,B1_1,B2_1,B3_1,B4_1,rss,ec_lda,ec_lda_rs) 524c **** PBE96 correlation energy corrections **** 525 t2 = t*t 526 t4 = t2*t2 527 B = -ec_lda/GAMMA 528 B = BOG/(exp(B)-1.0d0+ETA) 529 Q4 = 1.0d0 + B*t2 530 Q5 = 1.0d0 + B*t2 + B*B*t4 531 H = GAMMA*dlog(1.0d0 + BOG*Q4*t2/Q5) 532 533 534c **** PBE96 correlation fdn and fdnc derivatives **** 535 t6 = t4*t2 536 537 B_ec = (B/BETA)*(BOG+B) 538 539 Q8 = Q5*Q5+BOG*Q4*Q5*t2 540 Q9 = 1.0d0+2*B*t2 541 H_B = -BETA*B*t6*(2.0d0+B*t2)/Q8 542 Hrs = H_B*B_ec*ec_lda_rs 543 544 Ht = 2.0d0*BETA*Q9/Q8*t 545 546 ec = ec_lda + H 547 fnc = ec - (onethird*rs*ec_lda_rs) 548 > - (onethird*rs*Hrs) 549 > - (sevensixths*t*Ht) 550 fdnc = 0.5d0* Ht/ks 551 552 xce(i) = x_parameter*ex + c_parameter*ec 553 fn(i) = x_parameter*fnx + c_parameter*fnc 554 fdn(i) = x_parameter*fdnx + c_parameter*fdnc 555 556 557c write(*,*) "pbe96:",i,ec,fnc,fdnc 558 559 560 end do 561!$OMP END DO 562 563 return 564 end 565ccccccccccccccccccccccccccccccccccccccccccccccccccccccc 566 subroutine LSDT2(a,a1,b1,b2,b3,b4,srs,ec,ec_rs) 567 real*8 a,a1,b1,b2,b3,b4,srs,ec,ec_rs 568 real*8 q0,q1,q1p,qd,ql 569 q0= -2.0d0*a*(1.0d0+a1*srs*srs) 570 q1= 2.0d0*a*srs*(b1+srs*(b2+srs*(b3+srs*b4))) 571 q1p= a*((b1/srs)+2.0d0*b2+srs*(3.0d0*b3+srs*4.0d0*b4)) 572 qd=1.0d0/(q1*q1+q1) 573 ql= -dlog(qd*q1*q1) 574 ec= q0*ql 575 ec_rs= -2.0d0*a*a1*ql-q0*q1p*qd 576 return 577 end 578 579