1 2* ************************************ 3* * * 4* * gen_BEEF_BW_unrestricted * 5* * * 6* ************************************ 7* 8* This function returns the BEEF exchange-correlation 9* energy density, xce, and its derivatives with respect 10* to nup, ndn, |grad nup|, |grad ndn|, and |grad n|. 11* 12* Entry - n2ft3d : number of grid points 13* dn_in(*,2) : spin densites nup and ndn 14* agr_in(*,3): |grad nup|, |grad ndn|, and |grad n| 15* x_parameter: scale parameter for exchange 16* c_parameter: scale parameter for correlation 17* 18* Exit - xce(*) : PBE96 energy density 19* - fn(*,2) : d(n*xce)/dnup, d(n*xce)/dndn 20* - fdn(*,3): d(n*xce)/d|grad nup|, d(n*xce)/d|grad ndn| 21* d(n*xce)/d|grad n| 22 23 subroutine gen_BEEF_BW_unrestricted(n2ft3d, 24 > dn_in,agr_in, 25 > x_parameter,c_parameter,alphac, 26 > xce,fn,fdn) 27 implicit none 28 29 integer n2ft3d 30 real*8 dn_in(n2ft3d,2) 31 real*8 agr_in(n2ft3d,3) 32 real*8 x_parameter,c_parameter,alphac 33 real*8 xce(n2ft3d) 34 real*8 fn(n2ft3d,2) 35 real*8 fdn(n2ft3d,3) 36 37* **** Density cutoff parameter **** 38 real*8 DNS_CUT,ETA,ETA2,alpha_zeta,alpha_zeta2 39 parameter (DNS_CUT = 1.0d-20) 40 parameter (ETA=1.0d-20) 41 parameter (ETA2=1.0d-14) 42 parameter (alpha_zeta=(1.0d0-ETA2)) 43 parameter (alpha_zeta2=(1.0d0-ETA2)) 44 45c ***** PBE96 GGA exchange constants ****** 46 real*8 MU,KAPPA 47 parameter (MU = 0.2195149727645171d0) 48 parameter (KAPPA = 0.8040000000000000d0) 49 50c ****** PBE96 GGA correlation constants ****** 51 real*8 GAMMA,BETA,BOG 52 parameter (GAMMA = 0.031090690869655d0) 53 parameter (BETA = 0.066724550603149d0) 54 !parameter (BETA = 0.066725d0) 55 parameter (BOG = BETA/GAMMA) 56 57 58c ****** Perdew-Wang92 LDA correlation coefficients ******* 59 real*8 GAM,iGAM,FZZ,iFZZ 60 parameter (GAM = 0.519842099789746329d0) 61 parameter (iGAM = 1.0d0/GAM) 62 parameter (FZZ = (8.0d0/(9.0d0*GAM)) ) 63 parameter (iFZZ = 0.125d0*9.0d0*GAM) 64 65 real*8 A_1,A1_1,B1_1,B2_1,B3_1,B4_1 66 parameter (A_1 = 0.0310907d0) 67 !parameter (A_1 = 0.031091d0) 68 parameter (A1_1 = 0.2137000d0) 69 parameter (B1_1 = 7.5957000d0) 70 parameter (B2_1 = 3.5876000d0) 71 parameter (B3_1 = 1.6382000d0) 72 parameter (B4_1 = 0.4929400d0) 73 74 real*8 A_2,A1_2,B1_2,B2_2,B3_2,B4_2 75 parameter (A_2 = 0.01554535d0) 76 !parameter (A_2 = 0.015545d0) 77 parameter (A1_2 = 0.20548000d0) 78 parameter (B1_2 = 14.11890000d0) 79 parameter (B2_2 = 6.19770000d0) 80 parameter (B3_2 = 3.36620000d0) 81 parameter (B4_2 = 0.62517000d0) 82 83 real*8 A_3,A1_3,B1_3,B2_3,B3_3,B4_3 84 parameter (A_3 = 0.0168869d0) 85 !parameter (A_3 = 0.016887d0) 86 parameter (A1_3 = 0.1112500d0) 87 parameter (B1_3 = 10.3570000d0) 88 parameter (B2_3 = 3.6231000d0) 89 parameter (B3_3 = 0.8802600d0) 90 parameter (B4_3 = 0.4967100d0) 91 92c **** other constants **** 93 real*8 onethird,fourthird,fivethird,onesixthm 94 real*8 twothird,sevensixthm 95 real*8 onethirdm 96 parameter (onethird=1.0d0/3.0d0) 97 parameter (onethirdm=-1.0d0/3.0d0) 98 parameter (twothird=2.0d0/3.0d0) 99 parameter (fourthird=4.0d0/3.0d0) 100 parameter (fivethird=5.0d0/3.0d0) 101 parameter (onesixthm=-1.0d0/6.0d0) 102 parameter (sevensixthm=-7.0d0/6.0d0) 103 104*---- Beef expansion parameters given by Wellendorff et al-----------------* 105 real*8 malphac 106c real*8 alphac,malphac 107c parameter (alphac=0.6001664769d0) 108c parameter (malphac=1.0d0-alphac) 109 real*8 am(0:29) 110 data am(0:29)/1.516501714d0, 4.413532099d-1,-9.182135241d-2, 111 > -2.352754331d-2, 3.418828455d-2, 2.411870076d-3, 112 > -1.416381352d-2, 6.975895581d-4, 9.859205137d-3, 113 > -6.737855051d-3,-1.573330824d-3, 5.036146253d-3, 114 > -2.569472453d-3,-9.874953976d-4, 2.033722895d-3, 115 > -8.018718848d-4,-6.688078723d-4, 1.030936331d-3, 116 > -3.673838660d-4,-4.213635394d-4,5.761607992d-4, 117 > -8.346503735d-5,-4.458447585d-4,4.601290092d-4, 118 > -5.231775398d-6,-4.239570471d-4,3.750190679d-4, 119 > 2.114938125d-5, -1.904911565d-4,7.384362421d-5/ 120 121c **** local variables **** 122 integer i,j 123 real*8 n,agr 124 real*8 nup,agrup 125 real*8 ndn,agrdn 126 real*8 kf,ks,s,n_onethird,pi,rs_scale 127 real*8 rs ! Wigner radius 128 real*8 rss ! rss = sqrt(rs) 129 real*8 rs_n ! rs_n = n*drs/dn 130 real*8 t,t2,t4,t6 131 real*8 t_nup ! t_nup = n*dt/dnup 132 real*8 t_ndn ! t_ndn = n*dt/dndn 133 real*8 t_agr ! t_agr = n*dt/dagr 134 real*8 zet,twoksg 135 real*8 zet_nup ! zet_nup = n*dzet/dnup 136 real*8 zet_ndn ! zet_nup = n*dzet/dnup 137 real*8 zetp_1_3,zetm_1_3 138 real*8 zetpm_1_3,zetmm_1_3 139 real*8 phi,phi3,phi4 140 real*8 phi_zet 141 real*8 A,A2 142 real*8 A_phi,A_ec_lda 143 real*8 Q0,Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8 144 real*8 PON,FZ,z4 145 real*8 tau 146 real*8 F 147 real*8 Fs ! dF/ds 148 real*8 Hpbe 149 real*8 Hpbe_t ! dHpbe/dt 150 real*8 Hpbe_phi ! dHpbe/dphi 151 real*8 Hpbe_ec_lda ! dHpbe/d(ec_lda) 152 real*8 Hpbe_nup,Hpbe_ndn ! n*dHpbe/dnup, n*dHpbe/dndn 153 real*8 Ipbe 154 real*8 Ipbe_t,Ipbe_A ! dIpbe/dt, dIpbe/dA 155 156 real*8 exup,exdn,ex,ex_lda 157 real*8 ecu,ecp,eca,ec,ec_lda 158 real*8 ecu_rs,ecp_rs,eca_rs 159 real*8 ec_lda_rs,ec_lda_zet ! d(ec_lda)/drs, d(ec_lda)/dzet 160 real*8 ec_lda_nup,ec_lda_ndn ! n*d(ec_lda)/dnup, n*d(ec_lda)/dndn 161 real*8 fnxup,fdnxup ! d(n*ex)/dnup, d(n*ex)/dndn 162 real*8 fnxdn,fdnxdn ! d(n*ex)/d|grad nup|, d(n*ex)/d|grad ndn| 163 real*8 fncup,fncdn ! d(n*ec)/dnup, d(n*ec)/dndn 164 real*8 fdnx_const 165 166 real*8 p0,p1,p2,dp,s2,oneovers2,Ft,dp2,sgn 167 168 malphac=1.0d0-alphac 169 170 pi = 4.0d0*datan(1.0d0) 171 rs_scale = (0.75d0/pi)**onethird 172 fdnx_const = -3.0d0/(8.0d0*pi) 173 174 175!$OMP DO 176 do i=1,n2ft3d 177 nup = dn_in(i,1)+ETA 178 agrup = agr_in(i,1) 179 180 ndn = dn_in(i,2)+ETA 181 agrdn = agr_in(i,2) 182 183c **************************************************************** 184c ***** calculate polarized Exchange energies and potentials ***** 185c **************************************************************** 186 187c ************ 188c **** up **** 189c ************ 190 n = 2.0d0*nup 191 agr = 2.0d0*agrup 192 193 n_onethird = (3.0d0*n/pi)**onethird 194 ex_lda = -0.75d0*n_onethird 195 196 kf = (3.0d0*pi*pi*n)**onethird 197 s = agr/(2.0d0*kf*n) 198 199 t = 2.0d0*s*s/(4.0d0+s*s) - 1.0d0 200 if (dabs(t).gt.1.0d0) t = dsign(1.0d0,t) 201 F = 0.0d0 202 Ft = 0.0d0 203 204 s2 = t*t - 1.0d0 205 if (dabs(s2).lt.1.0d-12) then 206 if (t.gt.0.0d0) then 207 do j=0,29 208 F = F + am(j) 209 Ft = Ft + am(j)*0.5d0*dble(j*(j+1)) 210 end do 211 else 212 sgn = 1.0d0 213 do j=0,29 214 F = F + sgn*am(j) 215 Ft = Ft + sgn*am(j)*0.5d0*dble(j*(j+1)) 216 sgn = -sgn 217 end do 218 end if 219 else 220 oneovers2 = 1.0d0/s2 221 p0 = 1.0d0 222 dp = 0 223 F = F + am(0)*p0 224 Ft = Ft + am(0)*dp 225 p1 = t 226 dp = 1.0d0 227 F = F + am(1)*t 228 Ft = Ft + am(1)*dp 229 do j=2,29 230 p2 = (1.0d0/dble(j+1))*((2*j+1)*t*p1 - dble(j)*p0) 231 dp2 = dble(j)*oneovers2*(t*p2-p1) 232 F = F + am(j)*p2 233 Ft = Ft + am(j)*dp2 234 p0 = p1 235 p1 = p2 236 end do 237 end if 238 Fs = (16.0d0*s/(4.0d0+s*s)**2)*Ft 239 240 exup = ex_lda*F 241 fnxup = fourthird*(exup - ex_lda*Fs*s) 242 fdnxup = fdnx_const*Fs 243 244c ************** 245c **** down **** 246c ************** 247 n = 2.0d0*ndn 248 agr = 2.0d0*agrdn 249 250 n_onethird = (3.0d0*n/pi)**onethird 251 ex_lda = -0.75d0*n_onethird 252 253 kf = (3.0d0*pi*pi*n)**onethird 254 s = agr/(2.0d0*kf*n) 255 256 t = 2.0d0*s*s/(4.0d0+s*s) - 1.0d0 257 if (dabs(t).gt.1.0d0) t = dsign(1.0d0,t) 258 F = 0.0d0 259 Ft = 0.0d0 260 261 s2 = t*t - 1.0d0 262 if (dabs(s2).lt.1.0d-12) then 263 if (t.gt.0.0d0) then 264 do j=0,29 265 F = F + am(j) 266 Ft = Ft + am(j)*0.5d0*dble(j*(j+1)) 267 end do 268 else 269 sgn = 1.0d0 270 do j=0,29 271 F = F + sgn*am(j) 272 Ft = Ft + sgn*am(j)*0.5d0*dble(j*(j+1)) 273 sgn = -sgn 274 end do 275 end if 276 else 277 oneovers2 = 1.0d0/s2 278 p0 = 1.0d0 279 dp = 0 280 F = F + am(0)*p0 281 Ft = Ft + am(0)*dp 282 p1 = t 283 dp = 1.0d0 284 F = F + am(1)*t 285 Ft = Ft + am(1)*dp 286 do j=2,29 287 p2 = (1.0d0/dble(j+1))*((2*j+1)*t*p1 - dble(j)*p0) 288 dp2 = dble(j)*oneovers2*(t*p2-p1) 289 F = F + am(j)*p2 290 Ft = Ft + am(j)*dp2 291 p0 = p1 292 p1 = p2 293 end do 294 end if 295 Fs = (16.0d0*s/(4.0d0+s*s)**2)*Ft 296 297 exdn = ex_lda*F 298 fnxdn = fourthird*(exdn - ex_lda*Fs*s) 299 fdnxdn = fdnx_const*Fs 300 301 n = nup+ndn 302 303 ex = (exup*nup+ exdn*ndn)/ n 304 305c ******************************************************************* 306c ***** calculate polarized correlation energies and potentials ***** 307c ******************************************************************* 308 agr = agr_in(i,3) 309 310 zet = (nup-ndn)/n 311c if (zet.gt.0.0d0) zet = zet - ETA2 312c if (zet.lt.0.0d0) zet = zet + ETA2 313c if (dabs(dn_in(i,2)).gt.DNS_CUT) zet_nup = 2*ndn/n**2 314c if (dabs(dn_in(i,1)).gt.DNS_CUT) zet_ndn = -2*nup/n**2 315c if (dabs(dn_in(i,2)).gt.DNS_CUT) zet_nup = 2*ndn/n 316c zet_nup = 2*ndn/n 317c zet_ndn = -2*nup/n 318 zet_nup = -(zet - 1.0d0) 319 zet_ndn = -(zet + 1.0d0) 320 zetpm_1_3 = (1.0d0+zet*alpha_zeta)**onethirdm 321 zetmm_1_3 = (1.0d0-zet*alpha_zeta)**onethirdm 322 zetp_1_3 = (1.0d0+zet*alpha_zeta)*zetpm_1_3**2 323 zetm_1_3 = (1.0d0-zet*alpha_zeta)*zetmm_1_3**2 324 325 326 phi = 0.5d0*( zetp_1_3**2 + zetm_1_3**2) 327 phi_zet = alpha_zeta*( zetpm_1_3 - zetmm_1_3)/3.0d0 328 F =( (1.0d0+zet*alpha_zeta)*zetp_1_3 329 > + (1.0d0-zet*alpha_zeta)*zetm_1_3 330 > - 2.0d0)*iGAM 331 332 FZ = (zetp_1_3 - zetm_1_3)*(alpha_zeta*fourthird*iGAM) 333 334 335 336* **** calculate Wigner radius **** 337 rs = rs_scale/(n**onethird) 338 rss = dsqrt(rs) 339 340* **** calculate n*drs/dn **** 341c rs_n = onethirdm*rs/n 342 rs_n = onethirdm*rs 343 344 345 346c **** calculate t **** 347 kf = (3.0d0*pi*pi*n)**onethird 348 ks = dsqrt(4.0d0*kf/pi) 349 350 twoksg = 2.0d0*ks*phi 351 352 t = agr/(twoksg*n) 353 354* *** calculate n*dt/dnup, n*dt/dndn, n*dt/d|grad n| **** 355 t_nup = sevensixthm*t - (phi_zet)*(zet_nup)*t/phi 356 t_ndn = sevensixthm*t - (phi_zet)*(zet_ndn)*t/phi 357 t_agr = 1.0d0/(twoksg) 358 359 360 361 362c ************************************************** 363c ***** compute LSDA correlation energy density **** 364c ************************************************** 365 call LSDT(A_1,A1_1,B1_1,B2_1,B3_1,B4_1,rss,ecu,ecu_rs) 366 call LSDT(A_2,A1_2,B1_2,B2_2,B3_2,B4_2,rss,ecp,ecp_rs) 367 call LSDT(A_3,A1_3,B1_3,B2_3,B3_3,B4_3,rss,eca,eca_rs) 368 369 z4 = zet**4 370 371 ec_lda = ecu*(1.0d0-F*z4) 372 > + ecp*F*z4 373 > - eca*F*(1.0d0-z4)/FZZ 374 375 ec_lda_rs = ecu_rs*(1.0d0-F*z4) 376 > + ecp_rs*F*z4 377 > - eca_rs*F*(1.0d0-z4)/FZZ 378 379 ec_lda_zet = (4.0d0*(zet**3)*F + FZ*z4)*(ecp-ecu+eca*iFZZ) 380 > - FZ*eca*iFZZ 381 382 383 384 385c ******************************************** 386c **** calculate PBE96 correlation energy **** 387c ******************************************** 388 phi3 = phi**3 389 phi4 = phi3*phi 390 PON = -ec_lda/(phi3*GAMMA) 391 tau = DEXP(PON) 392 393 A = BOG/(tau-1.0d0+ETA) 394 A2 = A*A 395 t2 = t*t 396 t4 = t2*t2 397 t6 = t4*t2 398 Q4 = 1.0d0 + A*t2 399 Q5 = 1.0d0 + 2.0d0*A*t2 400 Q6 = 2.0d0 + A*t2 401 Q7 = 1.0d0+A*t2+A2*t4 402 Q8 = Q7*Q7 403 404 Ipbe = 1.0d0 + BOG*t2*Q4/Q7 405 Hpbe = GAMMA*phi3*DLOG(Ipbe) 406 407 Ipbe_t = BOG*(2.0d0*t)*Q5/Q8 408 Ipbe_A = -BOG*(A*t6) *Q6/Q8 409 410 A_ec_lda = tau/(BETA*phi3)*A2 411 A_phi = -3.0d0*ec_lda*tau/(BETA*phi4)*A2 412 413 414 Hpbe_ec_lda = (GAMMA*phi3/Ipbe)*Ipbe_A*A_ec_lda 415 416 Hpbe_phi = 3.0d0*Hpbe/phi 417 > + (GAMMA*phi3/Ipbe)*Ipbe_A*A_phi 418 419 Hpbe_t = (GAMMA*phi3/Ipbe)*Ipbe_t 420 421 ec_lda_nup = ec_lda_zet 422 > - zet * ec_lda_zet 423 > + rs_n * ec_lda_rs 424 ec_lda_ndn = -ec_lda_zet 425 > - zet * ec_lda_zet 426 > + rs_n * ec_lda_rs 427 428 429 430 Hpbe_nup = ec_lda_nup * Hpbe_ec_lda 431 > + phi_zet*zet_nup * Hpbe_phi 432 > + t_nup * Hpbe_t 433 434 Hpbe_ndn = ec_lda_ndn * Hpbe_ec_lda 435 > + phi_zet*zet_ndn * Hpbe_phi 436 > + t_ndn * Hpbe_t 437 438 439 440 ec = ec_lda + Hpbe 441 442 fncup = ec + (ec_lda_nup + Hpbe_nup) 443 fncdn = ec + (ec_lda_ndn + Hpbe_ndn) 444 445 xce(i) = x_parameter*ex + c_parameter*ec 446 fn(i,1) = x_parameter*fnxup + c_parameter*fncup 447 fn(i,2) = x_parameter*fnxdn + c_parameter*fncdn 448 449 fdn(i,1) = x_parameter*fdnxup 450 fdn(i,2) = x_parameter*fdnxdn 451 fdn(i,3) = c_parameter*t_agr*Hpbe_t 452 453 end do 454!$OMP END DO 455 456 457 458 return 459 end 460 461 462 463 464 465* ************************************ 466* * * 467* * gen_BEEF_BW_restricted * 468* * * 469* ************************************ 470 471* This routine calculates the non-Langreth terms of the BEEF-vdw exchange-correlation 472* potential(xcp) and energy density(xce). 473* 474* 475* Entry - n2ft3d : number of grid points 476* rho_in(*) : density (nup+ndn) 477* agr_in(*): |grad rho_in| 478* x_parameter: scale parameter for exchange 479* c_parameter: scale parameter for correlation 480* 481* Exit - xce(n2ft3d) : PBE96 exchange correlation energy density 482* fn(n2ft3d) : d(n*xce)/dn 483* fdn(n2ft3d) : d(n*xce/d|grad n| 484* 485 subroutine gen_BEEF_BW_restricted(n2ft3d,rho_in,agr_in, 486 > x_parameter,c_parameter,alphac, 487 > xce,fn,fdn) 488* 489 implicit none 490 integer n2ft3d 491 real*8 rho_in(n2ft3d) 492 real*8 agr_in(n2ft3d) 493 real*8 x_parameter,c_parameter,alphac 494 real*8 xce(n2ft3d) 495 real*8 fn(n2ft3d) 496 real*8 fdn(n2ft3d) 497 498* **** Density cutoff parameter **** 499 real*8 DNS_CUT,ETA 500 parameter (DNS_CUT = 1.0d-20) 501 parameter (ETA = 1.0d-20) 502 503c ***** PBE96 GGA exchange constants ****** 504 real*8 MU,KAPPA 505 parameter (MU = 0.2195149727645171d0) 506 parameter (KAPPA = 0.8040000000000000d0) 507 508c ****** PBE96 GGA correlation constants ****** 509 real*8 GAMMA,BETA,BOG 510 parameter (GAMMA = 0.031090690869655d0) 511 parameter (BETA = 0.066724550603149d0) 512 parameter (BOG = BETA/GAMMA) 513 514c ****** Perdew-Wang92 LDA correlation coefficients ******* 515 real*8 GAM,iGAM,FZZ,iFZZ 516 parameter (GAM = 0.519842099789746329d0) 517 parameter (iGAM = 1.0d0/GAM) 518 parameter (FZZ = (8.0d0/(9.0d0*GAM)) ) 519 parameter (iFZZ = 0.125d0*9.0d0*GAM) 520 521 real*8 A_1,A1_1,B1_1,B2_1,B3_1,B4_1 522 parameter (A_1 = 0.0310907d0) 523 !parameter (A_1 = 0.031091d0) 524 parameter (A1_1 = 0.2137000d0) 525 parameter (B1_1 = 7.5957000d0) 526 parameter (B2_1 = 3.5876000d0) 527 parameter (B3_1 = 1.6382000d0) 528 parameter (B4_1 = 0.4929400d0) 529 530 real*8 A_2,A1_2,B1_2,B2_2,B3_2,B4_2 531 parameter (A_2 = 0.01554535d0) 532 !parameter (A_2 = 0.015545d0) 533 parameter (A1_2 = 0.20548000d0) 534 parameter (B1_2 = 14.11890000d0) 535 parameter (B2_2 = 6.19770000d0) 536 parameter (B3_2 = 3.36620000d0) 537 parameter (B4_2 = 0.62517000d0) 538 539 real*8 A_3,A1_3,B1_3,B2_3,B3_3,B4_3 540 parameter (A_3 = 0.0168869d0) 541 !parameter (A_3 = 0.016887d0) 542 parameter (A1_3 = 0.1112500d0) 543 parameter (B1_3 = 10.3570000d0) 544 parameter (B2_3 = 3.6231000d0) 545 parameter (B3_3 = 0.8802600d0) 546 parameter (B4_3 = 0.4967100d0) 547 548c **** other constants **** 549 real*8 onethird,fourthird,fivethird,onesixthm 550 real*8 twothird,sevensixthm,sevensixths 551 real*8 onethirdm 552 parameter (onethird=1.0d0/3.0d0) 553 parameter (onethirdm=-1.0d0/3.0d0) 554 parameter (twothird=2.0d0/3.0d0) 555 parameter (fourthird=4.0d0/3.0d0) 556 parameter (fivethird=5.0d0/3.0d0) 557 parameter (onesixthm=-1.0d0/6.0d0) 558 parameter (sevensixthm=-7.0d0/6.0d0) 559 parameter (sevensixths=7.0d0/6.0d0) 560 561 562*---- Beef expansion parameters given by Wellendorff et al-----------------* 563 real*8 malphac 564c real*8 alphac,malphac 565c parameter (alphac=0.6001664769d0) 566c parameter (malphac=1.0d0-alphac) 567 real*8 am(0:29) 568 data am(0:29)/1.516501714d0, 4.413532099d-1,-9.182135241d-2, 569 > -2.352754331d-2, 3.418828455d-2, 2.411870076d-3, 570 > -1.416381352d-2, 6.975895581d-4, 9.859205137d-3, 571 > -6.737855051d-3,-1.573330824d-3, 5.036146253d-3, 572 > -2.569472453d-3,-9.874953976d-4, 2.033722895d-3, 573 > -8.018718848d-4,-6.688078723d-4, 1.030936331d-3, 574 > -3.673838660d-4,-4.213635394d-4,5.761607992d-4, 575 > -8.346503735d-5,-4.458447585d-4,4.601290092d-4, 576 > -5.231775398d-6,-4.239570471d-4,3.750190679d-4, 577 > 2.114938125d-5, -1.904911565d-4,7.384362421d-5/ 578 579c **** local variables **** 580 integer i,j 581 real*8 n,agr 582 real*8 kf,ks,s,n_onethird,pi,rs_scale 583 real*8 fdnx_const 584 real*8 rs,rss,t,t2,t4,t6 585 real*8 Q0,Q1,Q2,Q3,Q4,Q5,Q8,Q9,B 586 real*8 Ht 587 real*8 B_ec,Hrs,H_B 588 real*8 F,Fs 589 590 real*8 ex_lda,ec_lda 591 real*8 ec_lda_rs 592 real*8 ex,ec,H 593 real*8 fnx,fdnx,fnc,fdnc 594 595 596 real*8 p0,p1,p2,dp,s2,oneovers2,Ft,dp2,sgn 597 598 599 malphac=1.0d0-alphac 600 pi = 4.0d0*datan(1.0d0) 601 rs_scale = (0.75d0/pi)**onethird 602 fdnx_const = -3.0d0/(8.0d0*pi) 603 604!$OMP DO 605 do i=1,n2ft3d 606 n = rho_in(i)+ETA 607 agr = agr_in(i) 608 609c ***** calculate unpolarized Exchange energies and potentials ***** 610 n_onethird = (3.0d0*n/pi)**onethird 611 ex_lda = -0.75d0*n_onethird 612 613 kf = (3.0d0*pi*pi*n)**onethird 614 s = agr/(2.0d0*kf*n) 615 616 t = 2.0d0*s*s/(4.0d0+s*s) - 1.0d0 617 if (dabs(t).gt.1.0d0) t = dsign(1.0d0,t) 618 F = 0.0d0 619 Ft = 0.0d0 620 s2 = t*t - 1.0d0 621 if (dabs(s2).lt.1.0d-12) then 622 if (t.gt.0.0d0) then 623 do j=0,29 624 F = F + am(j) 625 Ft = Ft + am(j)*0.5d0*dble(j*(j+1)) 626 end do 627 else 628 sgn = 1.0d0 629 do j=0,29 630 F = F + sgn*am(j) 631 Ft = Ft + sgn*am(j)*0.5d0*dble(j*(j+1)) 632 sgn = -sgn 633 end do 634 end if 635 else 636 oneovers2 = 1.0d0/s2 637 p0 = 1.0d0 638 dp = 0 639 F = F + am(0)*p0 640 Ft = Ft + am(0)*dp 641 p1 = t 642 dp = 1.0d0 643 F = F + am(1)*t 644 Ft = Ft + am(1)*dp 645 do j=2,29 646 p2 = (1.0d0/dble(j+1))*((2*j+1)*t*p1 - dble(j)*p0) 647 dp2 = dble(j)*oneovers2*(t*p2-p1) 648 F = F + am(j)*p2 649 Ft = Ft + am(j)*dp2 650 p0 = p1 651 p1 = p2 652 end do 653 end if 654 Fs = (16.0d0*s/(4.0d0+s*s)**2)*Ft 655 656 ex = ex_lda*F 657 fnx = fourthird*(ex - ex_lda*Fs*s) 658 fdnx = fdnx_const*Fs 659 660 661* ********************************************************************* 662c ***** calculate unpolarized correlation energies and potentials ***** 663* ********************************************************************* 664 665c **** calculate rs and t **** 666 rs = rs_scale/(n**onethird) 667 rss = dsqrt(rs) 668 669 kf = (3.0d0*pi*pi*n)**onethird 670 ks = dsqrt(4.0d0*kf/pi) 671 t = agr/(2.0d0*ks*n) 672 673 674c **** unpolarized LDA correlation energy **** 675c **** ec_p = correlation energy **** 676c **** ec_p_rs = dec_p/drs **** 677c **** uc_p = dec_p/dn **** 678 call LSDT(A_1,A1_1,B1_1,B2_1,B3_1,B4_1,rss,ec_lda,ec_lda_rs) 679c **** PBE96 correlation energy corrections **** 680 t2 = t*t 681 t4 = t2*t2 682 B = -ec_lda/GAMMA 683 B = BOG/(exp(B)-1.0d0+ETA) 684 Q4 = 1.0d0 + B*t2 685 Q5 = 1.0d0 + B*t2 + B*B*t4 686 H = GAMMA*dlog(1.0d0 + BOG*Q4*t2/Q5) 687 688 689c **** PBE96 correlation fdn and fdnc derivatives **** 690 t6 = t4*t2 691 692 B_ec = (B/BETA)*(BOG+B) 693 694 Q8 = Q5*Q5+BOG*Q4*Q5*t2 695 Q9 = 1.0d0+2*B*t2 696 H_B = -BETA*B*t6*(2.0d0+B*t2)/Q8 697 Hrs = H_B*B_ec*ec_lda_rs 698 699 Ht = 2.0d0*BETA*Q9/Q8*t 700 701 ec = ec_lda + malphac*H 702 fnc = ec - (onethird*rs*ec_lda_rs) 703 > - malphac*(onethird*rs*Hrs) 704 > - malphac*(sevensixths*t*Ht) 705 fdnc = malphac*(0.5d0* Ht/ks) 706 707 xce(i) = x_parameter*ex + c_parameter*ec 708 fn(i) = x_parameter*fnx + c_parameter*fnc 709 fdn(i) = x_parameter*fdnx + c_parameter*fdnc 710 711 end do 712!$OMP END DO 713 714 return 715 end 716 717