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