1! 2MODULE exch_gga !<GPU:exch_gga=>exch_gga_gpu> 3! 4CONTAINS 5! 6!----------------------------------------------------------------------- 7SUBROUTINE becke88( rho, grho, sx, v1x, v2x ) !<GPU:DEVICE> 8 !----------------------------------------------------------------------- 9 !! Becke exchange: A.D. Becke, PRA 38, 3098 (1988) 10 !! only gradient-corrected part, no Slater term included 11 ! 12 USE kinds, ONLY: DP 13 ! 14 IMPLICIT NONE 15 ! 16 REAL(DP), INTENT(IN) :: rho, grho 17 REAL(DP), INTENT(OUT) :: sx, v1x, v2x 18 ! 19 ! ... local variables 20 ! 21 REAL(DP) :: rho13, rho43, xs, xs2, sa2b8, shm1, dd, dd2, ee 22 REAL(DP), PARAMETER :: beta=0.0042_DP 23 REAL(DP), PARAMETER :: third=1._DP/3._DP, two13=1.259921049894873_DP 24 ! two13= 2^(1/3) 25 ! 26 rho13 = rho**third 27 rho43 = rho13**4 28 ! 29 xs = two13 * SQRT(grho)/rho43 30 xs2 = xs * xs 31 ! 32 sa2b8 = SQRT(1.0_DP + xs2) 33 shm1 = LOG(xs + sa2b8) 34 ! 35 dd = 1.0_DP + 6.0_DP * beta * xs * shm1 36 dd2 = dd * dd 37 ! 38 ee = 6.0_DP * beta * xs2 / sa2b8 - 1._DP 39 sx = two13 * grho / rho43 * ( - beta / dd) 40 ! 41 v1x = - (4._DP/3._DP) / two13 * xs2 * beta * rho13 * ee / dd2 42 v2x = two13 * beta * (ee-dd) / (rho43 * dd2) 43 ! 44 RETURN 45 ! 46END SUBROUTINE becke88 47! 48! 49!----------------------------------------------------------------------- 50SUBROUTINE ggax( rho, grho, sx, v1x, v2x ) !<GPU:DEVICE> 51 !----------------------------------------------------------------------- 52 !! Perdew-Wang GGA (PW91), exchange part: 53 !! J.P. Perdew et al.,PRB 46, 6671 (1992) 54 ! 55 USE kinds, ONLY: DP 56 ! 57 IMPLICIT NONE 58 ! 59 REAL(DP), INTENT(IN) :: rho, grho 60 REAL(DP), INTENT(OUT) :: sx, v1x, v2x 61 ! 62 ! ... local variables 63 ! 64 REAL(DP) :: rhom43, s, s2, s3, s4, exps, as, sa2b8, shm1, bs, das, & 65 dbs, dls 66 REAL(DP), PARAMETER :: f1=0.19645_DP, f2=7.7956_DP, f3=0.2743_DP, & 67 f4=0.1508_DP, f5=0.004_DP 68 REAL(DP), PARAMETER :: fp1=-0.019292021296426_DP, fp2=0.161620459673995_DP 69 ! fp1= -3/(16 pi)*(3 pi^2)^(-1/3) 70 ! fp2= (1/2)(3 pi^2)**(-1/3) 71 ! 72 rhom43 = rho**(-4.d0/3.d0) 73 s = fp2 * SQRT(grho) * rhom43 74 s2 = s * s 75 s3 = s2 * s 76 s4 = s2 * s2 77 ! 78 exps = f4 * EXP( - 100.d0 * s2) 79 as = f3 - exps - f5 * s2 80 sa2b8 = SQRT(1.0d0 + f2 * f2 * s2) 81 shm1 = LOG(f2 * s + sa2b8) 82 bs = 1.d0 + f1 * s * shm1 + f5 * s4 83 ! 84 das = (200.d0 * exps - 2.d0 * f5) * s 85 dbs = f1 * (shm1 + f2 * s / sa2b8) + 4.d0 * f5 * s3 86 dls = (das / as - dbs / bs) 87 ! 88 sx = fp1 * grho * rhom43 * as / bs 89 v1x = - 4.d0 / 3.d0 * sx / rho * (1.d0 + s * dls) 90 v2x = fp1 * rhom43 * as / bs * (2.d0 + s * dls) 91 ! 92 RETURN 93 ! 94END SUBROUTINE ggax 95! 96! 97!--------------------------------------------------------------- 98SUBROUTINE pbex( rho, grho, iflag, sx, v1x, v2x ) !<GPU:DEVICE> 99 !--------------------------------------------------------------- 100 !! PBE exchange (without Slater exchange): 101 !! iflag=1 J.P.Perdew, K.Burke, M.Ernzerhof, PRL 77, 3865 (1996) 102 !! iflag=2 "revised' PBE: Y. Zhang et al., PRL 80, 890 (1998) 103 !! iflag=3 PBEsol: J.P.Perdew et al., PRL 100, 136406 (2008) 104 !! iflag=4 PBEQ2D: L. Chiodo et al., PRL 108, 126402 (2012) 105 !! iflag=5 optB88: Klimes et al., J. Phys. Cond. Matter, 22, 022201 (2010) 106 !! iflag=6 optB86b: Klimes et al., Phys. Rev. B 83, 195131 (2011) 107 !! iflag=7 ev: Engel and Vosko, PRB 47, 13164 (1991) 108 !! iflag=8 RPBE: B. Hammer, et al., Phys. Rev. B 59, 7413 (1999) 109 !! iflag=9 W31X: D. Chakraborty, K. Berland, and T. Thonhauser, JCTC 16, 5893 (2020) 110 ! 111 USE kinds, ONLY : DP 112 ! 113 IMPLICIT NONE 114 ! 115 INTEGER, INTENT(IN) :: iflag !<GPU:VALUE> 116 REAL(DP), INTENT(IN) :: rho, grho 117 ! input: charge and squared gradient 118 REAL(DP), INTENT(OUT) :: sx, v1x, v2x 119 ! output: energy, potential 120 ! 121 ! ... local variables 122 ! 123 REAL(DP) :: kf, agrho, s1, s2, ds, dsg, exunif, fx, sx_s 124 ! (3*pi2*|rho|)^(1/3) 125 ! |grho| 126 ! |grho|/(2*kf*|rho|) 127 ! s^2 128 ! n*ds/dn 129 ! n*ds/d(gn) 130 ! exchange energy LDA part 131 ! exchange energy gradient part 132 ! auxiliary variable for energy calculation 133 REAL(DP) :: dxunif, dfx, f1, f2, f3, dfx1 134 REAL(DP) :: p, amu, ab, c, dfxdp, dfxds, upbe, uge, s, ak, aa 135 ! numerical coefficients (NB: c2=(3 pi^2)^(1/3) ) 136 REAL(DP), PARAMETER :: pi=3.14159265358979323846d0 137 REAL(DP), PARAMETER :: third=1._DP/3._DP, c1=0.75_DP/pi, & 138 c2=3.093667726280136_DP, c5=4._DP*third, & 139 c6=c2*2.51984210_DP, c7=0.8_DP 140 ! (3pi^2)^(1/3)*2^(4/3) 141 ! parameters of the functional 142 REAL(DP) :: k(9), mu(9), ev(6) 143 ! pbe revpbe pbesol pbeq2d optB88 optB86b 144 ! ev rpbe W31x 145 DATA k / 0.804_DP, 1.2450_DP, 0.804_DP , 0.804_DP, 1.2_DP, 0.0_DP, & 146 0.000_DP, 0.8040_DP, 1.10_DP /, & 147 mu / 0.2195149727645171_DP, 0.2195149727645171_DP, 0.12345679012345679_DP, & 148 0.12345679012345679_DP, 0.22_DP, 0.1234_DP, 0.000_DP, & 149 0.2195149727645171_DP, 0.12345679012345679_DP /, & 150 ev / 1.647127_DP, 0.980118_DP, 0.017399_DP, 1.523671_DP, 0.367229_DP, & 151 0.011282_DP / ! a and b parameters of Engel and Vosko 152 ! 153 SELECT CASE( iflag ) 154 CASE( 4 ) 155 ! 156 agrho = SQRT(grho) 157 kf = c2 * rho**third 158 dsg = 0.5_DP / kf 159 s1 = agrho * dsg / rho 160 p = s1*s1 161 s = s1 162 ak = 0.804_DP 163 amu = 10._DP/81._DP 164 ab = 0.5217_DP 165 c = 2._DP 166 fx = ak - ak / (1.0_DP + amu * p / ak) + p**2 * (1.0_DP + p)/ & 167 (10**c + p**3) * ( -1.0_DP - ak + ak / (1.0_DP + amu * p / ak) & 168 + ab * p ** (-0.1D1/ 0.4D1) ) 169 ! 170 exunif = - c1 * kf 171 sx_s = exunif * fx 172 ! 173 dxunif = exunif * third 174 ! 175 dfxdp = DBLE(1 / (1 + amu * p / ak) ** 2 * amu) + DBLE(2 * p * (1 & 176 + p) / (10 ** c + p ** 3) * (-1 - ak + ak / (1 + amu * p / ak) + ab & 177 * p ** (-0.1d1 / 0.4D1))) + DBLE(p ** 2 / (10 ** c + p ** 3) * ( & 178 -1 - ak + ak / (1 + amu * p / ak) + ab * p ** (-0.1d1 / 0.4D1))) - & 179 DBLE(3 * p ** 4 * (1 + p) / (10 ** c + p ** 3) ** 2 * (-1 - ak + & 180 ak / (1 + amu * p / ak) + ab * p ** (-0.1d1 / 0.4D1))) + DBLE(p ** & 181 2) * DBLE(1 + p) / DBLE(10 ** c + p ** 3) * (-DBLE(1 / (1 + amu * & 182 p / ak) ** 2 * amu) - DBLE(ab * p ** (-0.5d1 / 0.4D1)) / 0.4D1) 183 ! 184 dfxds = dfxdp*2._DP*s 185 dfx = dfxds 186 ds = - c5 * s1 187 ! 188 v1x = sx_s + dxunif * fx + exunif * dfx * ds 189 v2x = exunif * dfx * dsg / agrho 190 sx = sx_s * rho 191 ! 192 CASE( 5, 9 ) 193 ! 194 agrho = SQRT(grho) 195 kf = c2 * rho**third 196 dsg = 0.5_DP / kf 197 s1 = agrho * dsg / rho 198 ab = mu(iflag) / k(iflag) 199 p = s1*c6 200 c = LOG(p + SQRT(p*p+1)) ! asinh(p) 201 dfx1 = 1 + ab*s1*c 202 fx = mu(iflag)*s1*s1/dfx1 203 ! 204 exunif = - c1 * kf 205 sx_s = exunif * fx 206 ! 207 dxunif = exunif * third 208 ! 209 dfx = 2*fx/s1-fx/dfx1*(ab*c+ab*s1/SQRT(p*p+1)*c6) 210 ds = - c5 * s1 211 ! 212 v1x = sx_s + dxunif * fx + exunif * dfx * ds 213 v2x = exunif * dfx * dsg / agrho 214 sx = sx_s * rho 215 ! 216 CASE( 6 ) 217 ! 218 agrho = SQRT(grho) 219 kf = c2 * rho**third 220 dsg = 0.5_DP / kf 221 s1 = agrho * dsg / rho 222 p = mu(iflag)*s1*s1 223 fx = p / ( 1._DP + p )**c7 224 ! 225 exunif = - c1 * kf 226 sx_s = exunif * fx 227 ! 228 dxunif = exunif * third 229 ! 230 dfx = 2*mu(iflag)*s1*fx*(1+(1-c7)*p)/(p*(1+p)) 231 ds = - c5 * s1 232 ! 233 v1x = sx_s + dxunif * fx + exunif * dfx * ds 234 v2x = exunif * dfx * dsg / agrho 235 sx = sx_s * rho 236 ! 237 CASE( 7 ) 238 ! 239 agrho = SQRT(grho) 240 kf = c2 * rho**third 241 dsg = 0.5_DP / kf 242 s1 = agrho * dsg / rho 243 s2 = s1 * s1 244 s = s2*s2 245 f1 = 1._DP + ev(1)*s2 + ev(2)*s + ev(3)*s*s2 246 f2 = 1._DP + ev(4)*s2 + ev(5)*s + ev(6)*s*s2 247 fx = f1 / f2 - 1._DP 248 ! 249 exunif = - c1 * kf 250 sx_s = exunif * fx 251 ! 252 dxunif = exunif * third 253 ds = - c5 * s1 254 ! 255 dfx = ev(1) + 2*ev(2)*s2 + 3*ev(3)*s 256 dfx1 = ev(4) + 2*ev(5)*s2 + 3*ev(6)*s 257 dfx = 2 * s1 * ( dfx - f1*dfx1/f2 ) / f2 258 ! 259 v1x = sx_s + dxunif * fx + exunif * dfx * ds 260 v2x = exunif * dfx * dsg / agrho 261 sx = sx_s * rho 262 ! 263 CASE(8) 264 ! 265 agrho = SQRT(grho) 266 kf = c2 * rho**third 267 dsg = 0.5_DP / kf 268 s1 = agrho * dsg / rho 269 s2 = s1 * s1 270 f1 = exp( - mu(iflag) * s2 / k(iflag) ) 271 f2 = 1._DP - f1 272 fx = k(iflag) * f2 273 ! 274 exunif = - c1 * kf 275 sx_s = exunif * fx 276 ! 277 dxunif = exunif * third 278 ds = - c5 * s1 279 ! 280 dfx = 2._DP * mu(iflag) * s1 * exp( - mu(iflag) * s2 / k(iflag) ) 281 ! 282 v1x = sx_s + dxunif * fx + exunif * dfx * ds 283 v2x = exunif * dfx * dsg / agrho 284 sx = sx_s * rho 285 ! 286 CASE DEFAULT 287 ! 288 agrho = SQRT(grho) 289 kf = c2 * rho**third 290 dsg = 0.5_DP / kf 291 s1 = agrho * dsg / rho 292 s2 = s1 * s1 293 f1 = s2 * mu(iflag) / k(iflag) 294 f2 = 1._DP + f1 295 f3 = k(iflag) / f2 296 fx = k(iflag) - f3 297 ! 298 exunif = - c1 * kf 299 sx_s = exunif * fx 300 ! 301 dxunif = exunif * third 302 ds = - c5 * s1 303 ! 304 dfx1 = f2 * f2 305 dfx = 2._DP * mu(iflag) * s1 / dfx1 306 ! 307 v1x = sx_s + dxunif * fx + exunif * dfx * ds 308 v2x = exunif * dfx * dsg / agrho 309 sx = sx_s * rho 310 ! 311 END SELECT 312 ! 313 ! 314 RETURN 315 ! 316END SUBROUTINE pbex 317! 318! 319!---------------------------------------------------------------------------- 320SUBROUTINE hcth( rho, grho, sx, v1x, v2x ) !<GPU:DEVICE> 321 !-------------------------------------------------------------------------- 322 !! HCTH/120, JCP 109, p. 6264 (1998) 323 !! Parameters set-up after N.L. Doltsisnis & M. Sprik (1999) 324 !! Present release: Mauro Boero, Tsukuba, 11/05/2004 325 ! 326 !! * rhoa = rhob = 0.5 * rho 327 !! * grho is the SQUARE of the gradient of rho! --> gr=sqrt(grho) 328 !! * sx : total exchange correlation energy at point r 329 !! * v1x : d(sx)/drho (eq. dfdra = dfdrb in original) 330 !! * v2x : 1/gr*d(sx)/d(gr) (eq. 0.5 * dfdza = 0.5 * dfdzb in original) 331 ! 332 USE kinds, ONLY: DP 333 ! 334 IMPLICIT NONE 335 ! 336 REAL(DP), INTENT(IN) :: rho, grho 337 REAL(DP), INTENT(OUT) :: sx, v1x, v2x 338 ! 339 ! ... local variables 340 ! 341 REAL(DP), PARAMETER :: pi=3.14159265358979323846d0 342 REAL(DP), PARAMETER :: o3 = 1.0d0/3.0d0, o34 = 4.0d0/3.0d0, fr83 = 8.d0/3.d0 343 REAL(DP) :: cg0(6), cg1(6), caa(6), cab(6), cx(6) 344 REAL(DP) :: r3q2, r3pi, gr, rho_o3, rho_o34, xa, xa2, ra, rab, & 345 dra_drho, drab_drho, g, dg, era1, dera1_dra, erab0, derab0_drab, & 346 ex, dex_drho, uaa, uab, ux, ffaa, ffab, dffaa_drho, dffab_drho, & 347 denaa, denab, denx, f83rho, bygr, gaa, gab, gx, taa, tab, txx, & 348 dgaa_drho, dgab_drho, dgx_drho, dgaa_dgr, dgab_dgr, dgx_dgr 349 ! 350 r3q2 = 2.d0**(-o3) 351 r3pi = (3.d0/pi)**o3 352 ! ... coefficients for pw correlation 353 cg0(1) = 0.031091d0 354 cg0(2) = 0.213700d0 355 cg0(3) = 7.595700d0 356 cg0(4) = 3.587600d0 357 cg0(5) = 1.638200d0 358 cg0(6) = 0.492940d0 359 cg1(1) = 0.015545d0 360 cg1(2) = 0.205480d0 361 cg1(3) =14.118900d0 362 cg1(4) = 6.197700d0 363 cg1(5) = 3.366200d0 364 cg1(6) = 0.625170d0 365 ! ... hcth-19-4 ... 366 caa(1) = 0.489508d+00 367 caa(2) = -0.260699d+00 368 caa(3) = 0.432917d+00 369 caa(4) = -0.199247d+01 370 caa(5) = 0.248531d+01 371 caa(6) = 0.200000d+00 372 cab(1) = 0.514730d+00 373 cab(2) = 0.692982d+01 374 cab(3) = -0.247073d+02 375 cab(4) = 0.231098d+02 376 cab(5) = -0.113234d+02 377 cab(6) = 0.006000d+00 378 cx(1) = 0.109163d+01 379 cx(2) = -0.747215d+00 380 cx(3) = 0.507833d+01 381 cx(4) = -0.410746d+01 382 cx(5) = 0.117173d+01 383 cx(6) = 0.004000d+00 384 ! ... ... ... ... ... 385 ! 386 gr = DSQRT(grho) 387 rho_o3 = rho**(o3) 388 rho_o34 = rho**(o34) 389 xa = 1.25992105d0*gr/rho_o34 390 xa2 = xa*xa 391 ra = 0.781592642d0/rho_o3 392 rab = r3q2*ra 393 dra_drho = -0.260530881d0/rho_o34 394 drab_drho = r3q2*dra_drho 395 CALL pwcorr( ra, cg1, g, dg ) !<GPU:pwcorr=>pwcorr_d> 396 era1 = g 397 dera1_dra = dg 398 CALL pwcorr( rab, cg0, g, dg ) !<GPU:pwcorr=>pwcorr_d> 399 erab0 = g 400 derab0_drab = dg 401 ex = -0.75d0*r3pi*rho_o34 402 dex_drho = -r3pi*rho_o3 403 uaa = caa(6)*xa2 404 uaa = uaa/(1.0d0+uaa) 405 uab = cab(6)*xa2 406 uab = uab/(1.0d0+uab) 407 ux = cx(6)*xa2 408 ux = ux/(1.0d0+ux) 409 ffaa = rho*era1 410 ffab = rho*erab0-ffaa 411 dffaa_drho = era1 + rho*dera1_dra*dra_drho 412 dffab_drho = erab0 + rho*derab0_drab*drab_drho - dffaa_drho 413 ! mb-> i-loop removed 414 denaa = 1.d0 / (1.0d0+caa(6)*xa2) 415 denab = 1.d0 / (1.0d0+cab(6)*xa2) 416 denx = 1.d0 / (1.0d0+cx(6)*xa2) 417 f83rho = fr83 / rho 418 bygr = 2.0d0/gr 419 gaa = caa(1)+uaa*(caa(2)+uaa*(caa(3)+uaa*(caa(4)+uaa*caa(5)))) 420 gab = cab(1)+uab*(cab(2)+uab*(cab(3)+uab*(cab(4)+uab*cab(5)))) 421 gx = cx(1)+ux*(cx(2)+ux*(cx(3)+ux*(cx(4)+ux*cx(5)))) 422 taa = denaa*uaa*(caa(2)+uaa*(2.d0*caa(3)+uaa & 423 *(3.d0*caa(4)+uaa*4.d0*caa(5)))) 424 tab = denab*uab*(cab(2)+uab*(2.d0*cab(3)+uab & 425 *(3.d0*cab(4)+uab*4.d0*cab(5)))) 426 txx = denx*ux*(cx(2)+ux*(2.d0*cx(3)+ux & 427 *(3.d0*cx(4)+ux*4.d0*cx(5)))) 428 dgaa_drho = -f83rho*taa 429 dgab_drho = -f83rho*tab 430 dgx_drho = -f83rho*txx 431 dgaa_dgr = bygr*taa 432 dgab_dgr = bygr*tab 433 dgx_dgr = bygr*txx 434 ! mb 435 sx = ex*gx + ffaa*gaa + ffab*gab 436 v1x = dex_drho*gx + ex*dgx_drho & 437 + dffaa_drho*gaa + ffaa*dgaa_drho & 438 + dffab_drho*gab + ffab*dgab_drho 439 v2x = (ex*dgx_dgr + ffaa*dgaa_dgr + ffab*dgab_dgr) / gr 440 ! 441 RETURN 442 ! 443END SUBROUTINE hcth 444 ! 445 !------------------------------------------------------- 446 SUBROUTINE pwcorr( r, c, g, dg ) !<GPU:DEVICE> 447 !----------------------------------------------------- 448 ! 449 USE kinds, ONLY: DP 450 ! 451 IMPLICIT NONE 452 ! 453 REAL(DP), INTENT(IN) :: r, c(6) 454 REAL(DP), INTENT(OUT) :: g, dg 455 ! 456 ! ... local variables 457 ! 458 REAL(DP) :: r12, r32, r2, rb, drb, sb 459 ! 460 r12 = DSQRT(r) 461 r32 = r*r12 462 r2 = r*r 463 rb = c(3)*r12 + c(4)*r + c(5)*r32 + c(6)*r2 464 sb = 1.0d0 + 1.0d0/(2.0d0*c(1)*rb) 465 g = -2.0d0 * c(1) * (1.0d0+c(2)*r) * DLOG(sb) 466 drb = c(3)/(2.0d0*r12) + c(4) + 1.5d0*c(5)*r12 + 2.0d0*c(6)*r 467 dg = (1.0d0+c(2)*r)*drb/(rb*rb*sb) - 2.0d0*c(1)*c(2)*DLOG(sb) 468 ! 469 RETURN 470 ! 471 END SUBROUTINE pwcorr 472! 473! 474!----------------------------------------------------------------------------- 475SUBROUTINE optx( rho, grho, sx, v1x, v2x ) !<GPU:DEVICE> 476 !--------------------------------------------------------------------------- 477 !! OPTX, Handy et al. JCP 116, p. 5411 (2002) and refs. therein 478 !! Present release: Mauro Boero, Tsukuba, 10/9/2002 479 ! 480 !! rhoa = rhob = 0.5 * rho in LDA implementation 481 !! grho is the SQUARE of the gradient of rho! --> gr=sqrt(grho) 482 !! sx : total exchange correlation energy at point r 483 !! v1x : d(sx)/drho 484 !! v2x : 1/gr*d(sx)/d(gr) 485 ! 486 USE kinds, ONLY: DP 487 ! 488 IMPLICIT NONE 489 ! 490 REAL(DP), INTENT(IN) :: rho, grho 491 REAL(DP), INTENT(OUT) :: sx, v1x, v2x 492 ! 493 ! ... local variables 494 ! 495 REAL(DP), PARAMETER :: small=1.D-30, smal2=1.D-10 496 ! ... coefficients and exponents 497 REAL(DP), PARAMETER :: o43=4.0d0/3.0d0, two13=1.259921049894873D0, & 498 two53=3.174802103936399D0, gam=0.006D0, a1cx=0.9784571170284421D0, & 499 a2=1.43169D0 500 REAL(DP) :: gr, rho43, xa, gamx2, uden, uu 501 ! 502 ! ... OPTX in compact form 503 ! 504 gr = MAX(grho,smal2) 505 rho43 = rho**o43 506 xa = two13*DSQRT(gr)/rho43 507 gamx2 = gam*xa*xa 508 uden = 1.d+00/(1.d+00+gamx2) 509 uu = a2*gamx2*gamx2*uden*uden 510 uden = rho43*uu*uden 511 sx = -rho43*(a1cx+uu)/two13 512 v1x = o43*(sx+two53*uden)/rho 513 v2x = -two53*uden/gr 514 ! 515 RETURN 516 ! 517END SUBROUTINE optx 518! 519! 520!--------------------------------------------------------------- 521SUBROUTINE wcx( rho, grho, sx, v1x, v2x ) !<GPU:DEVICE> 522 !--------------------------------------------------------------- 523 !! Wu-Cohen exchange (without Slater exchange): 524 !! Z. Wu and R. E. Cohen, PRB 73, 235116 (2006) 525 ! 526 USE kinds, ONLY: DP 527 ! 528 IMPLICIT NONE 529 ! 530 REAL(DP), INTENT(IN) :: rho, grho 531 REAL(DP), INTENT(OUT) :: sx, v1x, v2x 532 ! 533 ! ... local variables 534 ! 535 REAL(DP) :: kf, agrho, s1, s2, es2, ds, dsg, exunif, fx 536 ! (3*pi2*|rho|)^(1/3) 537 ! |grho| 538 ! |grho|/(2*kf*|rho|) 539 ! s^2 540 ! n*ds/dn 541 ! n*ds/d(gn) 542 ! exchange energy LDA part 543 ! exchange energy gradient part 544 REAL(DP) :: dxunif, dfx, f1, f2, f3, dfx1, x1, x2, x3, & 545 dxds1, dxds2, dxds3, sx_s 546 ! numerical coefficients (NB: c2=(3 pi^2)^(1/3) ) 547 REAL(DP), PARAMETER :: pi=3.14159265358979323846d0 548 REAL(DP), PARAMETER :: third=1.d0 / 3.d0, c1=0.75d0/pi , & 549 c2=3.093667726280136d0, c5=4.d0*third, & 550 teneightyone = 0.123456790123d0 551 ! parameters of the functional 552 REAL(DP), PARAMETER :: k=0.804d0, mu=0.2195149727645171d0, & 553 cwc=0.00793746933516d0 554 ! 555 agrho = SQRT(grho) 556 kf = c2 * rho**third 557 dsg = 0.5d0 / kf 558 s1 = agrho * dsg / rho 559 s2 = s1 * s1 560 es2 = EXP(-s2) 561 ds = - c5 * s1 562 ! 563 ! Energy 564 ! x = 10/81 s^2 + (mu - 10/81) s^2 e^-s^2 + ln (1 + c s^4) 565 x1 = teneightyone * s2 566 x2 = (mu - teneightyone) * s2 * es2 567 x3 = LOG(1.d0 + cwc * s2 * s2) 568 f1 = (x1 + x2 + x3) / k 569 f2 = 1.d0 + f1 570 f3 = k / f2 571 fx = k - f3 572 exunif = - c1 * kf 573 sx_s = exunif * fx 574 ! 575 ! Potential 576 dxunif = exunif * third 577 dfx1 = f2 * f2 578 dxds1 = teneightyone 579 dxds2 = (mu - teneightyone) * es2 * (1.d0 - s2) 580 dxds3 = 2.d0 * cwc * s2 / (1.d0 + cwc * s2 *s2) 581 dfx = 2.d0 * s1 * (dxds1 + dxds2 + dxds3) / dfx1 582 ! 583 v1x = sx_s + dxunif * fx + exunif * dfx * ds 584 v2x = exunif * dfx * dsg / agrho 585 sx = sx_s * rho 586 ! 587 RETURN 588 ! 589END SUBROUTINE wcx 590! 591! 592!----------------------------------------------------------------------- 593SUBROUTINE pbexsr( rho, grho, sxsr, v1xsr, v2xsr, omega ) !<GPU:DEVICE> 594 !--------------------------------------------------------------------- 595 ! INCLUDE 'cnst.inc' 596 USE kinds, ONLY: DP 597 ! 598 IMPLICIT NONE 599 ! 600 REAL(DP), INTENT(IN) :: omega !<GPU:VALUE> 601 REAL(DP), INTENT(IN) :: rho, grho 602 REAL(DP), INTENT(OUT) :: sxsr, v1xsr, v2xsr 603 ! 604 ! ... local variables 605 ! 606 REAL(DP) :: rs, vx, aa, rr, ex, s2, s, d1x, d2x, fx, dsdn, dsdg 607 REAL(DP), PARAMETER :: small=1.D-20, smal2=1.D-08 608 REAL(DP), PARAMETER :: us=0.161620459673995492D0, ax=-0.738558766382022406D0, & 609 um=0.2195149727645171D0, uk=0.8040D0, ul=um/uk 610 REAL(DP), PARAMETER :: f1=-1.10783814957303361_DP, alpha=2.0_DP/3.0_DP 611 ! 612 ! CALL XC(RHO,EX,EC,VX,VC) 613 ! 614 rs = rho**(1.0_DP/3.0_DP) 615 vx = (4.0_DP/3.0_DP)*f1*alpha*rs 616 ! 617 ! aa = dmax1(grho,smal2) 618 aa = grho 619 ! rr = rho**(-4.0_DP/3.0_DP) 620 rr = 1.0_DP/(rho*rs) 621 ex = ax/rr 622 s2 = aa*rr*rr*us*us 623 ! 624 s = SQRT(s2) 625 IF (s > 8.3D0) THEN 626 s = 8.572844D0 - 18.796223D0/s2 627 ENDIF 628 ! 629 CALL wpbe_analy_erfc_approx_grad( rho, s, omega, fx, d1x, d2x ) !<GPU:wpbe_analy_erfc_approx_grad=>wpbe_analy_erfc_approx_grad_d> 630 ! 631 sxsr = ex*fx ! - ex 632 dsdn = -4.D0/3.D0*s/rho 633 v1xsr = vx*fx + (dsdn*d2x+d1x)*ex ! - VX 634 dsdg = us*rr 635 v2xsr = ex*1.D0/SQRT(aa)*dsdg*d2x 636 ! 637 ! NOTE, here sx is the total energy density, 638 ! not just the gradient correction energy density as e.g. in pbex() 639 ! And the same goes for the potentials V1X, V2X 640 ! 641 RETURN 642 ! 643END SUBROUTINE pbexsr 644! 645! 646!----------------------------------------------------------------------- 647SUBROUTINE rPW86( rho, grho, sx, v1x, v2x ) !<GPU:DEVICE> 648 !--------------------------------------------------------------------- 649 !! PRB 33, 8800 (1986) and J. Chem. Theory comp. 5, 2754 (2009). 650 ! 651 USE kinds, ONLY: DP 652 ! 653 IMPLICIT NONE 654 ! 655 REAL(DP), INTENT(IN) :: rho, grho 656 REAL(DP), INTENT(OUT) :: sx, v1x, v2x 657 ! 658 ! ... local variables 659 ! 660 REAL(DP) :: s, s_2, s_3, s_4, s_5, s_6, fs, grad_rho, df_ds 661 REAL(DP), PARAMETER :: a=1.851_DP, b=17.33_DP, c=0.163_DP, & 662 s_prefactor=6.18733545256027_DP, & 663 Ax=-0.738558766382022_DP, four_thirds=4._DP/3._DP 664 ! 665 grad_rho = SQRT(grho) 666 ! 667 s = grad_rho/(s_prefactor*rho**(four_thirds)) 668 ! 669 s_2 = s**2 670 s_3 = s_2 * s 671 s_4 = s_2**2 672 s_5 = s_3 * s_2 673 s_6 = s_2 * s_4 674 ! 675 ! Calculation of energy 676 fs = (1 + a*s_2 + b*s_4 + c*s_6)**(1._DP/15._DP) 677 sx = Ax * rho**(four_thirds) * (fs -1.0_DP) 678 ! 679 ! Calculation of the potential 680 df_ds = (1._DP/(15._DP*fs**(14.0_DP)))*(2*a*s + 4*b*s_3 + 6*c*s_5) 681 ! 682 v1x = Ax*(four_thirds)*(rho**(1._DP/3._DP)*(fs -1.0_DP) & 683 -grad_rho/(s_prefactor * rho)*df_ds) 684 ! 685 v2x = Ax * df_ds/(s_prefactor*grad_rho) 686 ! 687END SUBROUTINE rPW86 688! 689! 690!----------------------------------------------------------------- 691SUBROUTINE c09x( rho, grho, sx, v1x, v2x ) !<GPU:DEVICE> 692 !--------------------------------------------------------------- 693 !! Cooper '09 exchange for vdW-DF (without Slater exchange): 694 !! V. R. Cooper, Phys. Rev. B 81, 161104(R) (2010) 695 ! 696 !! Developed thanks to the contribution of 697 !! Ikutaro Hamada - ikutaro@wpi-aimr.tohoku.ac.jp 698 !! WPI-Advanced Institute of Materials Research, Tohoku University 699 ! 700 USE kinds, ONLY: DP 701 ! 702 IMPLICIT NONE 703 ! 704 REAL(DP), INTENT(IN) :: rho, grho 705 REAL(DP), INTENT(OUT) :: sx, v1x, v2x 706 ! 707 ! ... local variables 708 ! 709 REAL(DP) :: kf, agrho, s1, s2, sx_s, ds, dsg, exunif, fx 710 ! (3*pi2*|rho|)^(1/3) 711 ! |grho| 712 ! |grho|/(2*kf*|rho|) 713 ! s^2 714 ! n*ds/dn 715 ! n*ds/d(gn) 716 ! exchange energy LDA part 717 ! exchange energy gradient part 718 REAL(DP) :: dxunif, dfx, f1, f2, f3, dfx1, dfx2 719 ! numerical coefficients (NB: c2=(3 pi^2)^(1/3) ) 720 REAL(DP), PARAMETER :: pi=3.14159265358979323846d0 721 REAL(DP), PARAMETER :: third=1._DP/3._DP, c1=0.75_DP/pi, & 722 c2=3.093667726280136_DP, c5=4._DP*third 723 ! parameters of the functional 724 REAL(DP) :: kappa, mu, alpha 725 DATA kappa / 1.245_DP /, & 726 mu / 0.0617_DP /, & 727 alpha / 0.0483_DP / 728 ! 729 agrho = SQRT(grho) 730 kf = c2 * rho**third 731 dsg = 0.5_DP / kf 732 s1 = agrho * dsg / rho 733 s2 = s1 * s1 734 ds = - c5 * s1 735 ! 736 ! ... Energy 737 ! 738 f1 = EXP( - alpha * s2 ) 739 f2 = EXP( - alpha * s2 / 2.0_DP ) 740 f3 = mu * s2 * f1 741 fx = f3 + kappa * ( 1.0_DP - f2 ) 742 exunif = - c1 * kf 743 sx_s = exunif * fx 744 ! 745 ! ... Potential 746 ! 747 dxunif = exunif * third 748 dfx1 = 2.0_DP * mu * s1 * ( 1.0_DP - alpha * s2 ) * f1 749 dfx2 = kappa * alpha * s1 * f2 750 dfx = dfx1 + dfx2 751 v1x = sx_s + dxunif * fx + exunif * dfx * ds 752 v2x = exunif * dfx * dsg / agrho 753 ! 754 sx = sx_s * rho 755 ! 756 RETURN 757 ! 758END SUBROUTINE c09x 759! 760! 761!--------------------------------------------------------------- 762SUBROUTINE sogga( rho, grho, sx, v1x, v2x ) !<GPU:DEVICE> 763 !------------------------------------------------------------- 764 !! SOGGA exchange 765 ! 766 USE kinds, ONLY: DP 767 ! 768 IMPLICIT NONE 769 ! 770 REAL(DP), INTENT(IN) :: rho, grho 771 REAL(DP), INTENT(OUT) :: sx, v1x, v2x 772 ! input: charge and abs gradient 773 ! output: energy and potential 774 ! 775 ! ... local variables 776 ! 777 REAL(DP) :: rho43, xs, xs2, dxs2_drho, dxs2_dgrho2 778 REAL(DP) :: CX, denom, C1, C2, Fso, Fpbe, ex, Fx, dFx_dxs2, dex_drho 779 ! 780 REAL(DP), PARAMETER :: one = 1.0_DP, two = 2.0_DP, three = 3.0_DP, & 781 & four = 4.0_DP, eight = 8.0_DP, & 782 & f13 = one/three, f23 = two/three, f43 = four/three, & 783 & f34 = three/four,f83 = eight/three, f12 = one/two 784 ! 785 REAL(DP), PARAMETER :: mu=0.12346_DP, kapa=0.552_DP 786 REAL(DP), PARAMETER :: pi=3.14159265358979323846d0 787 ! 788 ! Cx LDA 789 CX = f34 * (three/pi)**f13 790 denom = four * (three*pi**two)**f23 791 C1 = mu / denom 792 C2 = mu / (kapa * denom) 793 ! 794 rho43 = rho**f43 795 xs = grho / rho43 796 xs2 = xs * xs 797 ! 798 dxs2_drho = -f83 * xs2 / rho 799 dxs2_dgrho2 = one /rho**f83 800 ! 801 ex = - CX * rho43 802 dex_drho = - f43 * CX * rho**f13 803 ! 804 Fso = kapa * (one - EXP(-C2*xs2)) 805 Fpbe = C1 * xs2 / (one + C2*xs2) 806 ! 807 Fx = f12 * (Fpbe + Fso) 808 dFx_dxs2 = f12 * (C1 / ((one + C2*xs2)**2) + C1*EXP(-C2*xs2)) 809 ! 810 ! Energy 811 sx = Fx * ex 812 ! 813 ! Potential 814 v1x = dex_drho * Fx + ex * dFx_dxs2 * dxs2_drho 815 v2x = two * ex * dFx_dxs2 * dxs2_dgrho2 816 ! 817END SUBROUTINE sogga 818! 819! 820!------------------------------------------------------------------------- 821SUBROUTINE pbexgau( rho, grho, sxsr, v1xsr, v2xsr, alpha_gau ) !<GPU:DEVICE> 822 !----------------------------------------------------------------------- 823 ! 824 USE kinds, ONLY: DP 825 ! 826 IMPLICIT NONE 827 ! 828 REAL(DP), INTENT(IN) :: alpha_gau !<GPU:VALUE> 829 REAL(DP), INTENT(IN) :: rho, grho 830 REAL(DP), INTENT(OUT) :: sxsr, v1xsr, v2xsr 831 ! 832 ! ... local variables 833 ! 834 REAL(DP) :: rs, vx, aa, rr, ex, s2, s, d1x, d2x, fx, dsdn, dsdg 835 ! 836 REAL(DP), PARAMETER :: small=1.D-20, smal2=1.D-08 837 REAL(DP), PARAMETER :: us=0.161620459673995492D0, ax=-0.738558766382022406D0, & 838 um=0.2195149727645171D0, uk=0.8040D0, ul=um/uk 839 REAL(DP), PARAMETER :: f1 = -1.10783814957303361_DP, alpha = 2.0_DP/3.0_DP 840 ! 841 rs = rho**(1.0_DP/3.0_DP) 842 vx = (4.0_DP/3.0_DP)*f1*alpha*rs 843 aa = grho 844 rr = 1.0_DP/(rho*rs) 845 ex = ax/rr 846 ! AX is 3/4/PI*(3*PI*PI)**(1/3). This is the same as -c1*c2 in pbex(). 847 s2 = aa*rr*rr*us*us 848 s = SQRT(s2) 849 IF (s > 10.D0) THEN 850 s = 10.D0 851 ENDIF 852 CALL pbe_gauscheme( rho, s, alpha_gau, fx, d1x, d2x ) !<GPU:pbe_gauscheme=>pbe_gauscheme_d> 853 sxsr = ex*fx ! - EX 854 dsdn = -4.D0/3.D0*s/rho 855 v1xsr = vx*fx + (dsdn*d2x+d1x)*ex ! - VX 856 dsdg = us*rr 857 v2xsr = ex*1.D0/SQRT(aa)*dsdg*d2x 858 ! 859 ! NOTE, here sx is the total energy density, 860 ! not just the gradient correction energy density as e.g. in pbex() 861 ! And the same goes for the potentials V1X, V2X 862 ! 863 RETURN 864 ! 865END SUBROUTINE pbexgau 866 ! 867 !----------------------------------------------------------------------- 868SUBROUTINE pbe_gauscheme( rho, s, alpha_gau, Fx, dFxdr, dFxds ) !<GPU:DEVICE> 869 !-------------------------------------------------------------------- 870 ! 871 IMPLICIT NONE 872 ! 873 REAL*8 rho,s,alpha_gau,Fx,dFxdr,dFxds 874 ! input: charge and squared gradient and alpha_gau 875 ! output: GGA enhancement factor of gau-PBE 876 ! output: d(Fx)/d(s), d(Fx)/d(rho) 877 ! 878 REAL*8 Kx, Nx 879 ! PBE96 GGA enhancement factor 880 ! GGA enhancement factor of Gaussian Function 881 ! 882 REAL*8 bx, cx, PI, sqrtpial, Prefac, term_PBE, Third, KsF 883 REAL*8 d1sdr, d1Kxds, d1Kxdr, d1bxdr, d1bxds, d1bxdKx, & 884 d1Nxdbx,d1Nxdr, d1Nxds 885 ! 886 REAL*8 Zero,One,Two,Three,Four,Five,Six,Seven,Eight,Nine,Ten 887 ! 888 SAVE Zero,One,Two,Three,Four,Five,Six,Seven,Eight,Nine,Ten 889 DATA Zero,One,Two,Three,Four,Five,Six,Seven,Eight,Nine,Ten & 890 / 0D0,1D0,2D0,3D0,4D0,5D0,6D0,7D0,8D0,9D0,10D0 / 891 ! 892 REAL*8 k , mu 893 DATA k / 0.804d0 / , mu / 0.21951d0 / 894 ! parameters of PBE functional 895 ! 896 Third = One/Three 897 PI = ACOS(-One) 898 KsF = (Three*PI*PI*rho)**Third 899 sqrtpial = SQRT(PI/alpha_gau) 900 Prefac = Two * SQRT(PI/alpha_gau) / Three 901 ! 902 ! PBE96 GGA enhancement factor part 903 term_PBE = One / (One + s*s*mu/k) 904 Kx = One + k - k * term_PBE 905 ! 906 ! GGA enhancement factor of Gaussian Function part 907 bx = SQRT(Kx*alpha_gau) / KsF 908 ! 909 ! cx = exp(-One/Four/bx/bx) - One 910 IF (ABS(One/bx/bx) < 1.0D-4) THEN 911 cx = TayExp(-One/bx/bx) !<GPU:TayExp=>TayExp_d> 912 ELSE 913 cx = EXP(-One/bx/bx) - One 914 ENDIF 915 ! 916 Nx = bx * Prefac * ( SQRT(PI) * qe_erf(One/bx) + & !<GPU:qe_erf=>qe_erf_d> 917 (bx - Two*bx*bx*bx)*cx - Two*bx ) 918 ! 919 ! for convergency 920 IF (ABS(Nx) < 1.0D-15) THEN 921 Nx = Zero 922 ELSEIF ((One - ABS(Nx)) < 1.0D-15) THEN 923 Nx = One 924 ELSE 925 Nx = Nx 926 ENDIF 927 ! for convergency end 928 ! 929 Fx = Kx * Nx 930 ! 931 ! 1st derivatives 932 d1sdr = - Four / Three * s / rho 933 ! 934 d1Kxds = Two * s * mu * term_PBE * term_PBE 935 d1Kxdr = d1Kxds * d1sdr 936 d1bxdKx = bx / (Two* Kx) 937 ! 938 d1bxdr = - bx /(Three*rho) + d1Kxdr * d1bxdKx 939 ! 940 d1bxds = d1bxdKx * d1Kxds 941 ! 942 d1Nxdbx = Nx/bx - Prefac * bx * Three * & 943 ( cx*(One + Two*bx*bx) + Two ) 944 ! 945 d1Nxdr = d1Nxdbx * d1bxdr 946 d1Nxds = d1Nxdbx * d1bxds 947 ! 948 dFxdr = d1Kxdr * Nx + Kx * d1Nxdr 949 dFxds = d1Kxds * Nx + Kx * d1Nxds 950 ! 951 RETURN 952 ! 953END SUBROUTINE pbe_gauscheme 954! 955! 956!------------------------------------------------- 957FUNCTION TayExp(X) !<GPU:DEVICE> 958 !------------------------------------------- 959 USE kinds, ONLY: DP 960 IMPLICIT NONE 961 REAL(DP), INTENT(IN) :: X 962 REAL(DP) :: TAYEXP !<GPU:TAYEXP=>TAYEXP_d> 963 INTEGER :: NTERM,I 964 REAL(DP) :: SUMVAL,IVAL,COEF 965 PARAMETER (NTERM=16) 966 ! 967 SUMVAL = X 968 IVAL = X 969 COEF = 1.0D0 970 DO 10 I = 2, NTERM 971 COEF = COEF * I 972 IVAL = IVAL * (X / COEF) 973 SUMVAL = SUMVAL + IVAL 97410 CONTINUE 975 TAYEXP = SUMVAL !<GPU:TAYEXP=>TAYEXP_d> 976 ! 977 RETURN 978 ! 979END FUNCTION TayExp 980! 981! 982! 983!------------------------------------------------------------------------- 984SUBROUTINE PW86( rho, grho, sx, v1x, v2x ) !<GPU:DEVICE> 985 !----------------------------------------------------------------------- 986 !! Perdew-Wang 1986 exchange gradient correction: PRB 33, 8800 (1986) 987 ! 988 USE kinds, ONLY: DP 989 ! 990 IMPLICIT NONE 991 ! 992 REAL(DP), INTENT(IN) :: rho, grho 993 REAL(DP), INTENT(OUT) :: sx, v1x, v2x 994 ! 995 ! ... local variables 996 ! 997 REAL(DP) :: s, s_2, s_3, s_4, s_5, s_6, fs, grad_rho, df_ds 998 REAL(DP), PARAMETER :: a=1.296_DP, b=14._DP, c=0.2_DP, & 999 s_prefactor=6.18733545256027_DP, & 1000 Ax=-0.738558766382022_DP, four_thirds=4._DP/3._DP 1001 ! 1002 grad_rho = SQRT(grho) 1003 ! 1004 s = grad_rho / ( s_prefactor*rho**(four_thirds) ) 1005 ! 1006 s_2 = s**2 1007 s_3 = s_2 * s 1008 s_4 = s_2**2 1009 s_5 = s_3 * s_2 1010 s_6 = s_2 * s_4 1011 ! 1012 ! Calculation of energy 1013 fs = (1 + a*s_2 + b*s_4 + c*s_6)**(1._DP/15._DP) 1014 sx = Ax * rho**(four_thirds) * (fs-1._DP) 1015 ! 1016 ! Calculation of the potential 1017 df_ds = (1._DP/(15._DP*fs**(14._DP)))*(2*a*s + 4*b*s_3 + 6*c*s_5) 1018 ! 1019 v1x = Ax*(four_thirds)*( rho**(1._DP/3._DP)*(fs-1._DP) & 1020 -grad_rho/(s_prefactor * rho)*df_ds ) 1021 ! 1022 v2x = Ax * df_ds/(s_prefactor*grad_rho) 1023 ! 1024END SUBROUTINE PW86 1025! 1026! 1027!----------------------------------------------------------------------- 1028SUBROUTINE becke86b( rho, grho, sx, v1x, v2x ) !<GPU:DEVICE> 1029 !----------------------------------------------------------------------- 1030 !! Becke 1986 gradient correction to exchange 1031 !! A.D. Becke, J. Chem. Phys. 85 (1986) 7184 1032 ! 1033 USE kinds, ONLY: DP 1034 ! 1035 IMPLICIT NONE 1036 ! 1037 REAL(DP), INTENT(IN) :: rho, grho 1038 REAL(DP), INTENT(OUT) :: sx, v1x, v2x 1039 ! 1040 ! ... local variables 1041 ! 1042 REAL(DP) :: arho, agrho 1043 REAL(DP) :: sgp1, sgp1_45, sgp1_95 1044 REAL(DP) :: rdg2_43, rdg2_73, rdg2_83, rdg2_4, rdg4_5 1045 REAL(DP), PARAMETER :: beta=0.00375_DP, gamma=0.007_DP 1046 ! 1047 arho = 0.5_DP * rho 1048 agrho = 0.25_DP * grho 1049 ! 1050 rdg2_43 = agrho / arho**(4d0/3d0) 1051 rdg2_73 = rdg2_43 / arho 1052 rdg2_83 = rdg2_43 * rdg2_43 / agrho 1053 rdg2_4 = rdg2_43 * rdg2_83 / agrho 1054 rdg4_5 = rdg2_73 * rdg2_83 1055 ! 1056 sgp1 = 1d0 + gamma * rdg2_83 1057 sgp1_45 = sgp1**(-4d0/5d0) 1058 sgp1_95 = sgp1_45 / sgp1 1059 ! 1060 sx = -2d0 * beta * agrho / arho**(4d0/3d0) * sgp1_45 1061 v1x = -beta * (-4d0/3d0*rdg2_73*sgp1_45 + 32d0/15d0*gamma*rdg4_5*sgp1_95) 1062 v2x = -beta * (sgp1_45*rdg2_43/agrho - 4d0/5d0 *gamma*rdg2_4*sgp1_95) 1063 ! 1064END SUBROUTINE becke86b 1065! 1066! 1067!--------------------------------------------------------------- 1068SUBROUTINE b86b( rho, grho, iflag, sx, v1x, v2x ) !<GPU:DEVICE> 1069 !------------------------------------------------------------- 1070 !! Becke exchange (without Slater exchange): 1071 !! iflag=1: A. D. Becke, J. Chem. Phys. 85, 7184 (1986) (B86b) 1072 !! iflag=2: J. Klimes, Phys. Rev. B 83, 195131 (2011). (OptB86b) 1073 !! iflag=3: I. Hamada, Phys. Rev. B 89, 121103(R) (B86R) 1074 !! iflag=4: D. Chakraborty, K. Berland, and T. Thonhauser, JCTC 16, 5893 (2020) 1075 ! 1076 !! Ikutaro Hamada - HAMADA.Ikutaro@nims.go.jp 1077 !! National Institute for Materials Science 1078 ! 1079 USE kinds, ONLY : DP 1080 IMPLICIT NONE 1081 ! 1082 INTEGER, INTENT(IN) :: iflag !<GPU:VALUE> 1083 REAL(DP), INTENT(IN) :: rho, grho 1084 REAL(DP), INTENT(OUT) :: sx, v1x, v2x 1085 ! 1086 ! ... local variables 1087 ! 1088 REAL(DP) :: kf, agrho, s1, s2, sx_s, ds, dsg, exunif, fx 1089 ! (3*pi2*|rho|)^(1/3) 1090 ! |grho| 1091 ! |grho|/(2*kf*|rho|) 1092 ! s^2 1093 ! n*ds/dn 1094 ! n*ds/d(gn) 1095 ! exchange energy LDA part 1096 ! exchange energy gradient part 1097 REAL(DP) :: dxunif, dfx, f1, f2, f3, dfx1 1098 ! numerical coefficients (NB: c2=(3 pi^2)^(1/3) ) 1099 REAL(DP), PARAMETER :: pi=3.14159265358979323846d0 1100 REAL(DP), PARAMETER :: third=1._DP/3._DP, c1=0.75_DP/pi, & 1101 c2=3.093667726280136_DP, c5=4._DP*third 1102 ! parameters of the functional 1103 REAL(DP) :: k(4), mu(4) 1104 DATA k / 0.5757_DP, 1.0000_DP, 0.711357_DP, 0.58_DP /, & 1105 mu/ 0.2449_DP, 0.1234_DP, 0.1234_DP, 0.12345679012345679_DP / 1106 ! 1107 agrho = SQRT(grho) 1108 kf = c2 * rho**third 1109 dsg = 0.5_DP / kf 1110 s1 = agrho * dsg / rho 1111 s2 = s1 * s1 1112 ds = - c5 * s1 1113 ! 1114 ! ... Energy 1115 ! 1116 f1 = mu(iflag)*s2 1117 f2 = 1._DP + mu(iflag)*s2/k(iflag) 1118 f3 = f2**(4._DP/5._DP) 1119 fx = f1/f3 1120 exunif = - c1 * kf 1121 sx_s = exunif * fx 1122 ! 1123 ! ... Potential 1124 ! 1125 dxunif = exunif * third 1126 dfx1 = 1._DP + (1._DP/5._DP)*mu(iflag)*s2 / k(iflag) 1127 dfx = 2._DP * mu(iflag) * s1 * dfx1 / (f2 * f3) 1128 v1x = sx_s + dxunif * fx + exunif * dfx * ds 1129 v2x = exunif * dfx * dsg / agrho 1130 sx = sx_s * rho 1131 ! 1132 RETURN 1133 ! 1134END SUBROUTINE b86b 1135! 1136! 1137!----------------------------------------------------------------------- 1138SUBROUTINE cx13( rho, grho, sx, v1x, v2x ) !<GPU:DEVICE> 1139 !----------------------------------------------------------------------- 1140 !! The new exchange partner for a vdW-DF1-cx suggested 1141 !! by K. Berland and P. Hyldgaard, see PRB 89, 035412 (2014), 1142 !! to test the plasmon nature of the vdW-DF1 inner functional. 1143 ! 1144 USE kinds, ONLY : DP 1145 ! 1146 IMPLICIT NONE 1147 ! 1148 REAL(DP), INTENT(IN) :: rho, grho 1149 REAL(DP), INTENT(OUT) :: sx, v1x, v2x 1150 ! 1151 ! ... local variables 1152 ! 1153 REAL(DP) :: s, s_2, s_3, s_4, s_5, s_6, fs, fs_rPW86, df_rPW86_ds, grad_rho, df_ds 1154 REAL(DP), PARAMETER :: alp=0.021789_DP, beta=1.15_DP, a=1.851_DP, b=17.33_DP, & 1155 c=0.163_DP, mu_LM=0.09434_DP, & 1156 s_prefactor=6.18733545256027_DP, & 1157 Ax = -0.738558766382022_DP, four_thirds = 4._DP/3._DP 1158 ! 1159 grad_rho = SQRT(grho) 1160 ! 1161 s = grad_rho/(s_prefactor*rho**(four_thirds)) 1162 ! 1163 s_2 = s*s 1164 s_3 = s_2 * s 1165 s_4 = s_2 * s_2 1166 s_5 = s_3 * s_2 1167 s_6 = s_2 * s_2 *s_2 1168 ! 1169 ! ... Energy 1170 fs_rPW86 = (1._DP + a*s_2 + b*s_4 + c*s_6)**(1._DP/15._DP) 1171 fs = 1._DP/(1._DP + alp*s_6) * (1._DP + mu_LM *s_2) & 1172 + alp*s_6/(beta+alp*s_6)*fs_rPW86 1173 ! 1174 sx = Ax * rho**(four_thirds) * (fs-1._DP) 1175 ! 1176 ! ... Potential 1177 df_rPW86_ds = (1._DP/(15._DP*fs_rPW86**(14._DP)))*(2*a*s + 4*b*s_3 + 6*c*s_5) 1178 ! 1179 df_ds = 1._DP/(1._DP+alp*s_6)**2*( 2._DP*mu_LM*s*(1._DP+alp*s_6) & 1180 - 6._DP*alp*s_5*( 1._DP+mu_LM*s_2) ) & 1181 + alp*s_6/(beta+alp*s_6)*df_rPW86_ds & 1182 + 6._DP*alp*s_5*fs_rPW86/(beta+alp*s_6)*(1._DP-alp*s_6/(beta + alp*s_6)) 1183 ! 1184 v1x = Ax*(four_thirds)*(rho**(1._DP/3._DP)*(fs-1._DP) & 1185 -grad_rho/(s_prefactor * rho)*df_ds) 1186 v2x = Ax * df_ds/(s_prefactor*grad_rho) 1187 ! 1188END SUBROUTINE cx13 1189! 1190! 1191! 1192! ===========> SPIN <=========== 1193! 1194!----------------------------------------------------------------------- 1195SUBROUTINE becke88_spin( rho_up, rho_dw, grho_up, grho_dw, sx_up, sx_dw, v1x_up, v1x_dw, v2x_up, v2x_dw ) !<GPU:DEVICE> 1196 !----------------------------------------------------------------------- 1197 !! Becke exchange: A.D. Becke, PRA 38, 3098 (1988) - Spin polarized case 1198 ! 1199 USE kinds, ONLY: DP 1200 ! 1201 IMPLICIT NONE 1202 ! 1203 REAL(DP), INTENT(IN) :: rho_up, rho_dw 1204 !! charge 1205 REAL(DP), INTENT(IN) :: grho_up, grho_dw 1206 !! gradient 1207 REAL(DP), INTENT(OUT) :: sx_up, sx_dw 1208 !! the up and down energies 1209 REAL(DP), INTENT(OUT) :: v1x_up, v1x_dw 1210 !! first part of the potential 1211 REAL(DP), INTENT(OUT) :: v2x_up, v2x_dw 1212 !! second part of the potential 1213 ! 1214 ! ... local variables 1215 ! 1216 INTEGER :: is 1217 REAL(DP), PARAMETER :: beta = 0.0042_DP, third = 1._DP/3._DP 1218 REAL(DP) :: rho13, rho43, xs, xs2, sa2b8, shm1, dd, dd2, ee 1219 ! 1220 ! 1221 !DO is = 1, 2 1222 rho13 = rho_up**third 1223 rho43 = rho13**4 1224 xs = SQRT(grho_up) / rho43 1225 xs2 = xs * xs 1226 sa2b8 = SQRT(1.0d0 + xs2) 1227 shm1 = LOG(xs + sa2b8) 1228 dd = 1.0d0 + 6.0d0 * beta * xs * shm1 1229 dd2 = dd * dd 1230 ee = 6.0d0 * beta * xs2 / sa2b8 - 1.d0 1231 sx_up = grho_up / rho43 * (-beta/dd) 1232 v1x_up = -(4.d0/3.d0) * xs2 * beta * rho13 * ee / dd2 1233 v2x_up = beta * (ee-dd) / (rho43*dd2) 1234 1235 rho13 = rho_dw**third 1236 rho43 = rho13**4 1237 xs = SQRT(grho_dw) / rho43 1238 xs2 = xs * xs 1239 sa2b8 = SQRT(1.0d0 + xs2) 1240 shm1 = LOG(xs + sa2b8) 1241 dd = 1.0d0 + 6.0d0 * beta * xs * shm1 1242 dd2 = dd * dd 1243 ee = 6.0d0 * beta * xs2 / sa2b8 - 1.d0 1244 sx_dw = grho_dw / rho43 * (-beta/dd) 1245 v1x_dw = -(4.d0/3.d0) * xs2 * beta * rho13 * ee / dd2 1246 v2x_dw = beta * (ee-dd) / (rho43*dd2) 1247 !ENDDO 1248 ! 1249 RETURN 1250 ! 1251END SUBROUTINE becke88_spin 1252! 1253! 1254!----------------------------------------------------------------------------- 1255SUBROUTINE wpbe_analy_erfc_approx_grad( rho, s, omega, Fx_wpbe, d1rfx, d1sfx ) !<GPU:DEVICE> 1256 !----------------------------------------------------------------------- 1257 ! 1258 ! wPBE Enhancement Factor (erfc approx.,analytical, gradients) 1259 ! 1260 !-------------------------------------------------------------------- 1261 ! 1262 USE kinds, ONLY: DP 1263 IMPLICIT NONE 1264 ! 1265 REAL(DP) rho,s,omega,Fx_wpbe,d1sfx,d1rfx 1266 ! 1267 REAL(DP) f12,f13,f14,f18,f23,f43,f32,f72,f34,f94,f1516,f98 1268 REAL(DP) pi,pi2,pi_23,srpi 1269 REAL(DP) Three_13 1270 ! 1271 REAL(DP) ea1,ea2,ea3,ea4,ea5,ea6,ea7,ea8 1272 REAL(DP) eb1 1273 REAL(DP) A,B,C,D,E 1274 REAL(DP) Ha1,Ha2,Ha3,Ha4,Ha5 1275 REAL(DP) Fc1,Fc2 1276 REAL(DP) EGa1,EGa2,EGa3 1277 REAL(DP) EGscut,wcutoff,expfcutoff 1278 ! 1279 REAL(DP) xkf, xkfrho 1280 REAL(DP) w,w2,w3,w4,w5,w6,w7,w8 1281 REAL(DP) d1rw 1282 REAL(DP) A2,A3,A4,A12,A32,A52,A72 1283 REAL(DP) X 1284 REAL(DP) s2,s3,s4,s5,s6 1285 ! 1286 REAL(DP) H,F 1287 REAL(DP) Hnum,Hden,d1sHnum,d1sHden 1288 REAL(DP) d1sH,d1sF 1289 REAL(DP) G_a,G_b,EG 1290 REAL(DP) d1sG_a,d1sG_b,d1sEG 1291 ! 1292 REAL(DP) Hsbw,Hsbw2,Hsbw3,Hsbw4,Hsbw12,Hsbw32,Hsbw52,Hsbw72 1293 REAL(DP) DHsbw,DHsbw2,DHsbw3,DHsbw4,DHsbw5 1294 REAL(DP) DHsbw12,DHsbw32,DHsbw52,DHsbw72,DHsbw92 1295 REAL(DP) d1sHsbw,d1rHsbw 1296 REAL(DP) d1sDHsbw,d1rDHsbw 1297 REAL(DP) HsbwA94,HsbwA9412 1298 REAL(DP) HsbwA942,HsbwA943,HsbwA945 1299 REAL(DP) piexperf,expei 1300 REAL(DP) piexperfd1,expeid1 1301 REAL(DP) d1spiexperf,d1sexpei 1302 REAL(DP) d1rpiexperf,d1rexpei 1303 REAL(DP) expei1,expei2,expei3,expei4 1304 ! 1305 REAL(DP) DHs,DHs2,DHs3,DHs4,DHs72,DHs92,DHsw,DHsw2,DHsw52,DHsw72 1306 REAL(DP) d1sDHs,d1rDHsw 1307 ! 1308 REAL(DP) np1,np2 1309 REAL(DP) d1rnp1,d1rnp2 1310 REAL(DP) t1,t2t9,t10,t10d1 1311 REAL(DP) f2,f3,f4,f5,f6,f7,f8,f9 1312 REAL(DP) f2d1,f3d1,f4d1,f5d1,f6d1,f8d1,f9d1 1313 REAL(DP) d1sf2,d1sf3,d1sf4,d1sf5,d1sf6,d1sf7,d1sf8,d1sf9 1314 REAL(DP) d1rf2,d1rf3,d1rf4,d1rf5,d1rf6,d1rf7,d1rf8,d1rf9 1315 REAL(DP) d1st1,d1rt1 1316 REAL(DP) d1st2t9,d1rt2t9 1317 REAL(DP) d1st10,d1rt10 1318 REAL(DP) d1sterm1,d1rterm1,term1d1 1319 REAL(DP) d1sterm2 1320 REAL(DP) d1sterm3,d1rterm3 1321 REAL(DP) d1sterm4,d1rterm4 1322 REAL(DP) d1sterm5,d1rterm5 1323 ! 1324 REAL(DP) term1,term2,term3,term4,term5 1325 ! 1326 REAL(DP) ax,um,uk,ul 1327 REAL(DP) gc1,gc2 1328 ! 1329 ! REAL(DP) ei 1330 ! 1331 REAL(DP) Zero,One,Two,Three,Four,Five,Six,Seven,Eight,Nine,Ten 1332 REAL(DP) Fifteen,Sixteen 1333 REAL(DP) r12,r64,r36,r81,r256,r384,r864,r1944,r4374 1334 REAL(DP) r20,r25,r27,r48,r120,r128,r144,r288,r324,r512,r729 1335 REAL(DP) r30,r32,r75,r243,r2187,r6561,r40,r105,r54,r135 1336 REAL(DP) r1215,r15309 1337 ! 1338 SAVE Zero,One,Two,Three,Four,Five,Six,Seven,Eight,Nine,Ten 1339 DATA Zero,One,Two,Three,Four,Five,Six,Seven,Eight,Nine,Ten & 1340 / 0D0,1D0,2D0,3D0,4D0,5D0,6D0,7D0,8D0,9D0,10D0 / 1341 SAVE Fifteen,Sixteen 1342 DATA Fifteen,Sixteen / 1.5D1, 1.6D1 / 1343 SAVE r36,r64,r81,r256,r384,r864,r1944,r4374 1344 DATA r36,r64,r81,r256,r384,r864,r1944,r4374 & 1345 / 3.6D1,6.4D1,8.1D1,2.56D2,3.84D2,8.64D2,1.944D3,4.374D3 / 1346 SAVE r27,r48,r120,r128,r144,r288,r324,r512,r729 1347 DATA r27,r48,r120,r128,r144,r288,r324,r512,r729 & 1348 / 2.7D1,4.8D1,1.2D2,1.28D2,1.44D2,2.88D2,3.24D2,5.12D2,7.29D2 / 1349 SAVE r20,r32,r243,r2187,r6561,r40 1350 DATA r20,r32,r243,r2187,r6561,r40 & 1351 / 2.0d1,3.2D1,2.43D2,2.187D3,6.561D3,4.0d1 / 1352 SAVE r12,r25,r30,r54,r75,r105,r135,r1215,r15309 1353 DATA r12,r25,r30,r54,r75,r105,r135,r1215,r15309 & 1354 / 1.2D1,2.5d1,3.0d1,5.4D1,7.5d1,1.05D2,1.35D2,1.215D3,1.5309D4 / 1355 ! 1356 ! ... General constants 1357 ! 1358 f12 = 0.5d0 1359 f13 = One/Three 1360 f14 = 0.25d0 1361 f18 = 0.125d0 1362 ! 1363 f23 = Two * f13 1364 f43 = Two * f23 1365 ! 1366 f32 = 1.5d0 1367 f72 = 3.5d0 1368 f34 = 0.75d0 1369 f94 = 2.25d0 1370 f98 = 1.125d0 1371 f1516 = Fifteen / Sixteen 1372 ! 1373 pi = ACOS(-One) 1374 pi2 = pi*pi 1375 pi_23 = pi2**f13 1376 srpi = SQRT(pi) 1377 ! 1378 Three_13 = Three**f13 1379 ! 1380 ! Constants from fit 1381 ! 1382 ea1 = -1.128223946706117d0 1383 ea2 = 1.452736265762971d0 1384 ea3 = -1.243162299390327d0 1385 ea4 = 0.971824836115601d0 1386 ea5 = -0.568861079687373d0 1387 ea6 = 0.246880514820192d0 1388 ea7 = -0.065032363850763d0 1389 ea8 = 0.008401793031216d0 1390 ! 1391 eb1 = 1.455915450052607d0 1392 ! 1393 ! Constants for PBE hole 1394 ! 1395 A = 1.0161144d0 1396 B = -3.7170836d-1 1397 C = -7.7215461d-2 1398 D = 5.7786348d-1 1399 E = -5.1955731d-2 1400 X = - Eight/Nine 1401 ! 1402 ! Constants for fit of H(s) (PBE) 1403 ! 1404 Ha1 = 9.79681d-3 1405 Ha2 = 4.10834d-2 1406 Ha3 = 1.87440d-1 1407 Ha4 = 1.20824d-3 1408 Ha5 = 3.47188d-2 1409 ! 1410 ! Constants for F(H) (PBE) 1411 ! 1412 Fc1 = 6.4753871d0 1413 Fc2 = 4.7965830d-1 1414 ! 1415 ! Constants for polynomial expansion for EG for small s 1416 ! 1417 EGa1 = -2.628417880d-2 1418 EGa2 = -7.117647788d-2 1419 EGa3 = 8.534541323d-2 1420 ! 1421 ! Constants for large x expansion of exp(x)*ei(-x) 1422 ! 1423 expei1 = 4.03640D0 1424 expei2 = 1.15198D0 1425 expei3 = 5.03627D0 1426 expei4 = 4.19160D0 1427 ! 1428 ! Cutoff criterion below which to use polynomial expansion 1429 ! 1430 EGscut = 8.0d-2 1431 wcutoff = 1.4D1 1432 expfcutoff = 7.0D2 1433 ! 1434 ! Calculate prelim variables 1435 ! 1436 xkf = (Three*pi2*rho) ** f13 1437 xkfrho = xkf * rho 1438 ! 1439 A2 = A*A 1440 A3 = A2*A 1441 A4 = A3*A 1442 A12 = SQRT(A) 1443 A32 = A12*A 1444 A52 = A32*A 1445 A72 = A52*A 1446 ! 1447 w = omega / xkf 1448 w2 = w * w 1449 w3 = w2 * w 1450 w4 = w2 * w2 1451 w5 = w3 * w2 1452 w6 = w5 * w 1453 w7 = w6 * w 1454 w8 = w7 * w 1455 ! 1456 d1rw = -(One/(Three*rho))*w 1457 ! 1458 X = - Eight/Nine 1459 ! 1460 s2 = s*s 1461 s3 = s2*s 1462 s4 = s2*s2 1463 s5 = s4*s 1464 s6 = s5*s 1465 ! 1466 ! Calculate wPBE enhancement factor 1467 ! 1468 Hnum = Ha1*s2 + Ha2*s4 1469 Hden = One + Ha3*s4 + Ha4*s5 + Ha5*s6 1470 ! 1471 H = Hnum/Hden 1472 ! 1473 d1sHnum = Two*Ha1*s + Four*Ha2*s3 1474 d1sHden = Four*Ha3*s3 + Five*Ha4*s4 + Six*Ha5*s5 1475 ! 1476 d1sH = (Hden*d1sHnum - Hnum*d1sHden) / (Hden*Hden) 1477 ! 1478 F = Fc1*H + Fc2 1479 d1sF = Fc1*d1sH 1480 ! 1481 ! Change exponent of Gaussian if we're using the simple approx. 1482 ! 1483 IF (w > wcutoff) eb1 = 2.0d0 1484 ! 1485 ! Calculate helper variables (should be moved later on...) 1486 ! 1487 Hsbw = s2*H + eb1*w2 1488 Hsbw2 = Hsbw*Hsbw 1489 Hsbw3 = Hsbw2*Hsbw 1490 Hsbw4 = Hsbw3*Hsbw 1491 Hsbw12 = SQRT(Hsbw) 1492 Hsbw32 = Hsbw12*Hsbw 1493 Hsbw52 = Hsbw32*Hsbw 1494 Hsbw72 = Hsbw52*Hsbw 1495 ! 1496 d1sHsbw = d1sH*s2 + Two*s*H 1497 d1rHsbw = Two*eb1*d1rw*w 1498 ! 1499 DHsbw = D + s2*H + eb1*w2 1500 DHsbw2 = DHsbw*DHsbw 1501 DHsbw3 = DHsbw2*DHsbw 1502 DHsbw4 = DHsbw3*DHsbw 1503 DHsbw5 = DHsbw4*DHsbw 1504 DHsbw12 = SQRT(DHsbw) 1505 DHsbw32 = DHsbw12*DHsbw 1506 DHsbw52 = DHsbw32*DHsbw 1507 DHsbw72 = DHsbw52*DHsbw 1508 DHsbw92 = DHsbw72*DHsbw 1509 ! 1510 HsbwA94 = f94 * Hsbw / A 1511 HsbwA942 = HsbwA94*HsbwA94 1512 HsbwA943 = HsbwA942*HsbwA94 1513 HsbwA945 = HsbwA943*HsbwA942 1514 HsbwA9412 = SQRT(HsbwA94) 1515 ! 1516 DHs = D + s2*H 1517 DHs2 = DHs*DHs 1518 DHs3 = DHs2*DHs 1519 DHs4 = DHs3*DHs 1520 DHs72 = DHs3*SQRT(DHs) 1521 DHs92 = DHs72*DHs 1522 ! 1523 d1sDHs = Two*s*H + s2*d1sH 1524 ! 1525 DHsw = DHs + w2 1526 DHsw2 = DHsw*DHsw 1527 DHsw52 = SQRT(DHsw)*DHsw2 1528 DHsw72 = DHsw52*DHsw 1529 ! 1530 d1rDHsw = Two*d1rw*w 1531 ! 1532 IF (s > EGscut) THEN 1533 ! 1534 G_a = srpi * (Fifteen*E + Six*C*(One+F*s2)*DHs + & 1535 Four*B*(DHs2) + Eight*A*(DHs3)) & 1536 * (One / (Sixteen * DHs72)) & 1537 - f34*pi*SQRT(A) * EXP(f94*H*s2/A) * & 1538 (One - qe_erf(f32*s*SQRT(H/A))) !<GPU:qe_erf=>qe_erf_d> 1539 ! 1540 d1sG_a = (One/r32)*srpi * & 1541 ((r36*(Two*H + d1sH*s) / (A12*SQRT(H/A))) & 1542 + (One/DHs92) * & 1543 (-Eight*A*d1sDHs*DHs3 - r105*d1sDHs*E & 1544 -r30*C*d1sDHs*DHs*(One+s2*F) & 1545 +r12*DHs2*(-B*d1sDHs + C*s*(d1sF*s + Two*F))) & 1546 - ((r54*EXP(f94*H*s2/A)*srpi*s*(Two*H+d1sH*s)* & 1547 qe_erfc(f32*SQRT(H/A)*s)) & !<GPU:qe_erfc=>qe_erfc_d> 1548 / A12)) 1549 ! 1550 G_b = (f1516 * srpi * s2) / DHs72 1551 ! 1552 d1sG_b = (Fifteen*srpi*s*(Four*DHs - Seven*d1sDHs*s)) & 1553 / (r32*DHs92) 1554 ! 1555 EG = - (f34*pi + G_a) / G_b 1556 ! 1557 d1sEG = (-Four*d1sG_a*G_b + d1sG_b*(Four*G_a + Three*pi)) & 1558 / (Four*G_b*G_b) 1559 ! 1560 ELSE 1561 ! 1562 EG = EGa1 + EGa2*s2 + EGa3*s4 1563 d1sEG = Two*EGa2*s + Four*EGa3*s3 1564 ! 1565 ENDIF 1566 ! 1567 ! Calculate the terms needed in any case 1568 ! 1569 term2 = (DHs2*B + DHs*C + Two*E + DHs*s2*C*F + Two*s2*EG) / & 1570 (Two*DHs3) 1571 ! 1572 d1sterm2 = (-Six*d1sDHs*(EG*s2 + E) & 1573 + DHs2 * (-d1sDHs*B + s*C*(d1sF*s + Two*F)) & 1574 + Two*DHs * (Two*EG*s - d1sDHs*C & 1575 + s2 * (d1sEG - d1sDHs*C*F))) & 1576 / (Two*DHs4) 1577 1578 term3 = - w * (Four*DHsw2*B + Six*DHsw*C + Fifteen*E & 1579 + Six*DHsw*s2*C*F + Fifteen*s2*EG) / & 1580 (Eight*DHs*DHsw52) 1581 ! 1582 d1sterm3 = w * (Two*d1sDHs*DHsw * (Four*DHsw2*B & 1583 + Six*DHsw*C + Fifteen*E & 1584 + Three*s2*(Five*EG + Two*DHsw*C*F)) & 1585 + DHs * (r75*d1sDHs*(EG*s2 + E) & 1586 + Four*DHsw2*(d1sDHs*B & 1587 - Three*s*C*(d1sF*s + Two*F)) & 1588 - Six*DHsw*(-Three*d1sDHs*C & 1589 + s*(Ten*EG + Five*d1sEG*s & 1590 - Three*d1sDHs*s*C*F)))) & 1591 / (Sixteen*DHs2*DHsw72) 1592 ! 1593 d1rterm3 = (-Two*d1rw*DHsw * (Four*DHsw2*B & 1594 + Six*DHsw*C + Fifteen*E & 1595 + Three*s2*(Five*EG + Two*DHsw*C*F)) & 1596 + w * d1rDHsw * (r75*(EG*s2 + E) & 1597 + Two*DHsw*(Two*DHsw*B + Nine*C & 1598 + Nine*s2*C*F))) & 1599 / (Sixteen*DHs*DHsw72) 1600 1601 term4 = - w3 * (DHsw*C + Five*E + DHsw*s2*C*F + Five*s2*EG) / & 1602 (Two*DHs2*DHsw52) 1603 ! 1604 d1sterm4 = (w3 * (Four*d1sDHs*DHsw * (DHsw*C + Five*E & 1605 + s2 * (Five*EG + DHsw*C*F)) & 1606 + DHs * (r25*d1sDHs*(EG*s2 + E) & 1607 - Two*DHsw2*s*C*(d1sF*s + Two*F) & 1608 + DHsw * (Three*d1sDHs*C + s*(-r20*EG & 1609 - Ten*d1sEG*s & 1610 + Three*d1sDHs*s*C*F))))) & 1611 / (Four*DHs3*DHsw72) 1612 ! 1613 d1rterm4 = (w2 * (-Six*d1rw*DHsw * (DHsw*C + Five*E & 1614 + s2 * (Five*EG + DHsw*C*F)) & 1615 + w * d1rDHsw * (r25*(EG*s2 + E) + & 1616 Three*DHsw*C*(One + s2*F)))) & 1617 / (Four*DHs2*DHsw72) 1618 ! 1619 term5 = - w5 * (E + s2*EG) / & 1620 (DHs3*DHsw52) 1621 ! 1622 d1sterm5 = (w5 * (Six*d1sDHs*DHsw*(EG*s2 + E) & 1623 + DHs * (-Two*DHsw*s * (Two*EG + d1sEG*s) & 1624 + Five*d1sDHs * (EG*s2 + E)))) & 1625 / (Two*DHs4*DHsw72) 1626 ! 1627 d1rterm5 = (w4 * Five*(EG*s2 + E) * (-Two*d1rw*DHsw & 1628 + d1rDHsw * w)) & 1629 / (Two*DHs3*DHsw72) 1630 ! 1631 ! 1632 IF ((s > 0.0d0).OR.(w > 0.0d0)) THEN 1633 ! 1634 t10 = (f12)*A*LOG(Hsbw / DHsbw) 1635 t10d1 = f12*A*(One/Hsbw - One/DHsbw) 1636 d1st10 = d1sHsbw*t10d1 1637 d1rt10 = d1rHsbw*t10d1 1638 ! 1639 ENDIF 1640 ! 1641 ! Calculate exp(x)*f(x) depending on size of x 1642 ! 1643 IF (HsbwA94 < expfcutoff) THEN 1644 ! 1645 piexperf = pi*EXP(HsbwA94)*qe_erfc(HsbwA9412) !<GPU:qe_erfc=>qe_erfc_d> 1646 ! expei = Exp(HsbwA94)*Ei(-HsbwA94) 1647 expei = EXP(HsbwA94)*(-expint(1,HsbwA94)) !<GPU:expint=>expint_d> 1648 1649 ELSE 1650 ! 1651 ! print *,rho,s," LARGE HsbwA94" 1652 ! 1653 piexperf = pi*(One/(srpi*HsbwA9412) & 1654 - One/(Two*SQRT(pi*HsbwA943)) & 1655 + Three/(Four*SQRT(pi*HsbwA945))) 1656 ! 1657 expei = - (One/HsbwA94) * & 1658 (HsbwA942 + expei1*HsbwA94 + expei2) / & 1659 (HsbwA942 + expei3*HsbwA94 + expei4) 1660 1661 ENDIF 1662 ! 1663 ! Calculate the derivatives (based on the orig. expression) 1664 ! --> Is this ok? ==> seems to be ok... 1665 ! 1666 piexperfd1 = - (Three*srpi*SQRT(Hsbw/A))/(Two*Hsbw) & 1667 + (Nine*piexperf)/(Four*A) 1668 d1spiexperf = d1sHsbw*piexperfd1 1669 d1rpiexperf = d1rHsbw*piexperfd1 1670 1671 expeid1 = f14*(Four/Hsbw + (Nine*expei)/A) 1672 d1sexpei = d1sHsbw*expeid1 1673 d1rexpei = d1rHsbw*expeid1 1674 ! 1675 IF (w == Zero) THEN 1676 ! 1677 ! Fall back to original expression for the PBE hole 1678 ! 1679 t1 = -f12*A*expei 1680 d1st1 = -f12*A*d1sexpei 1681 d1rt1 = -f12*A*d1rexpei 1682 ! 1683 ! write(*,*) s, t1, t10, d1st1,d1rt1,d1rt10 1684 ! 1685 IF (s > 0.0D0) THEN 1686 ! 1687 term1 = t1 + t10 1688 d1sterm1 = d1st1 + d1st10 1689 d1rterm1 = d1rt1 + d1rt10 1690 ! 1691 Fx_wpbe = X * (term1 + term2) 1692 ! 1693 d1sfx = X * (d1sterm1 + d1sterm2) 1694 d1rfx = X * d1rterm1 1695 ! 1696 ELSE 1697 ! 1698 Fx_wpbe = 1.0d0 1699 ! 1700 ! TODO This is checked to be true for term1 1701 ! How about the other terms??? 1702 ! 1703 d1sfx = 0.0d0 1704 d1rfx = 0.0d0 1705 ! 1706 ENDIF 1707 ! 1708 ! 1709 ELSEIF (w > wcutoff) THEN 1710 ! 1711 ! Use simple Gaussian approximation for large w 1712 ! 1713 ! print *,rho,s," LARGE w" 1714 ! 1715 term1 = -f12*A*(expei+LOG(DHsbw)-LOG(Hsbw)) 1716 1717 term1d1 = - A/(Two*DHsbw) - f98*expei 1718 d1sterm1 = d1sHsbw*term1d1 1719 d1rterm1 = d1rHsbw*term1d1 1720 1721 Fx_wpbe = X * (term1 + term2 + term3 + term4 + term5) 1722 1723 d1sfx = X * (d1sterm1 + d1sterm2 + d1sterm3 & 1724 + d1sterm4 + d1sterm5) 1725 1726 d1rfx = X * (d1rterm1 + d1rterm3 + d1rterm4 + d1rterm5) 1727 ! 1728 ELSE 1729 ! 1730 ! For everything else, use the full blown expression 1731 ! 1732 ! First, we calculate the polynomials for the first term 1733 ! 1734 np1 = -f32*ea1*A12*w + r27*ea3*w3/(Eight*A12) & 1735 - r243*ea5*w5/(r32*A32) + r2187*ea7*w7/(r128*A52) 1736 ! 1737 d1rnp1 = - f32*ea1*d1rw*A12 + (r81*ea3*d1rw*w2)/(Eight*A12) & 1738 - (r1215*ea5*d1rw*w4)/(r32*A32) & 1739 + (r15309*ea7*d1rw*w6)/(r128*A52) 1740 ! 1741 np2 = -A + f94*ea2*w2 - r81*ea4*w4/(Sixteen*A) & 1742 + r729*ea6*w6/(r64*A2) - r6561*ea8*w8/(r256*A3) 1743 ! 1744 ! 1745 d1rnp2 = f12*(Nine*ea2*d1rw*w) & 1746 - (r81*ea4*d1rw*w3)/(Four*A) & 1747 + (r2187*ea6*d1rw*w5)/(r32*A2) & 1748 - (r6561*ea8*d1rw*w7)/(r32*A3) 1749 ! 1750 ! The first term is 1751 ! 1752 t1 = f12*(np1*piexperf + np2*expei) 1753 d1st1 = f12*(d1spiexperf*np1 + d1sexpei*np2) 1754 d1rt1 = f12*(d1rnp2*expei + d1rpiexperf*np1 + & 1755 d1rexpei*np2 + d1rnp1*piexperf) 1756 ! 1757 ! The factors for the main polynomoal in w and their derivatives 1758 ! 1759 f2 = (f12)*ea1*srpi*A / DHsbw12 1760 f2d1 = - ea1*srpi*A / (Four*DHsbw32) 1761 d1sf2 = d1sHsbw*f2d1 1762 d1rf2 = d1rHsbw*f2d1 1763 ! 1764 f3 = (f12)*ea2*A / DHsbw 1765 f3d1 = - ea2*A / (Two*DHsbw2) 1766 d1sf3 = d1sHsbw*f3d1 1767 d1rf3 = d1rHsbw*f3d1 1768 ! 1769 f4 = ea3*srpi*(-f98 / Hsbw12 & 1770 + f14*A / DHsbw32) 1771 f4d1 = ea3*srpi*((Nine/(Sixteen*Hsbw32))- & 1772 (Three*A/(Eight*DHsbw52))) 1773 d1sf4 = d1sHsbw*f4d1 1774 d1rf4 = d1rHsbw*f4d1 1775 ! 1776 f5 = ea4*(One/r128) * (-r144*(One/Hsbw) & 1777 + r64*(One/DHsbw2)*A) 1778 f5d1 = ea4*((f98/Hsbw2)-(A/DHsbw3)) 1779 d1sf5 = d1sHsbw*f5d1 1780 d1rf5 = d1rHsbw*f5d1 1781 ! 1782 f6 = ea5*(Three*srpi*(Three*DHsbw52*(Nine*Hsbw-Two*A) & 1783 + Four*Hsbw32*A2)) & 1784 / (r32*DHsbw52*Hsbw32*A) 1785 f6d1 = ea5*srpi*((r27/(r32*Hsbw52))- & 1786 (r81/(r64*Hsbw32*A))- & 1787 ((Fifteen*A)/(Sixteen*DHsbw72))) 1788 d1sf6 = d1sHsbw*f6d1 1789 d1rf6 = d1rHsbw*f6d1 1790 ! 1791 f7 = ea6*(((r32*A)/DHsbw3 & 1792 + (-r36 + (r81*s2*H)/A)/Hsbw2)) / r32 1793 d1sf7 = ea6*(Three*(r27*d1sH*DHsbw4*Hsbw*s2 + & 1794 Eight*d1sHsbw*A*(Three*DHsbw4 - Four*Hsbw3*A) + & 1795 r54*DHsbw4*s*(Hsbw - d1sHsbw*s)*H))/ & 1796 (r32*DHsbw4*Hsbw3*A) 1797 d1rf7 = ea6*d1rHsbw*((f94/Hsbw3)-((Three*A)/DHsbw4) & 1798 -((r81*s2*H)/(Sixteen*Hsbw3*A))) 1799 ! 1800 f8 = ea7*(-Three*srpi*(-r40*Hsbw52*A3 & 1801 +Nine*DHsbw72*(r27*Hsbw2-Six*Hsbw*A+Four*A2))) & 1802 / (r128 * DHsbw72*Hsbw52*A2) 1803 f8d1 = ea7*srpi*((r135/(r64*Hsbw72)) + (r729/(r256*Hsbw32*A2)) & 1804 -(r243/(r128*Hsbw52*A)) & 1805 -((r105*A)/(r32*DHsbw92))) 1806 d1sf8 = d1sHsbw*f8d1 1807 d1rf8 = d1rHsbw*f8d1 1808 ! 1809 f9 = (r324*ea6*eb1*DHsbw4*Hsbw*A & 1810 + ea8*(r384*Hsbw3*A3 + DHsbw4*(-r729*Hsbw2 & 1811 + r324*Hsbw*A - r288*A2))) / (r128*DHsbw4*Hsbw3*A2) 1812 f9d1 = -((r81*ea6*eb1)/(Sixteen*Hsbw3*A)) & 1813 + ea8*((r27/(Four*Hsbw4))+(r729/(r128*Hsbw2*A2)) & 1814 -(r81/(Sixteen*Hsbw3*A)) & 1815 -((r12*A/DHsbw5))) 1816 d1sf9 = d1sHsbw*f9d1 1817 d1rf9 = d1rHsbw*f9d1 1818 ! 1819 t2t9 = f2*w + f3*w2 + f4*w3 + f5*w4 + f6*w5 & 1820 + f7*w6 + f8*w7 + f9*w8 1821 d1st2t9 = d1sf2*w + d1sf3*w2 + d1sf4*w3 + d1sf5*w4 & 1822 + d1sf6*w5 + d1sf7*w6 + d1sf8*w7 & 1823 + d1sf9*w8 1824 d1rt2t9 = d1rw*f2 + d1rf2*w + Two*d1rw*f3*w & 1825 + d1rf3*w2 + Three*d1rw*f4*w2 & 1826 + d1rf4*w3 + Four*d1rw*f5*w3 & 1827 + d1rf5*w4 + Five*d1rw*f6*w4 & 1828 + d1rf6*w5 + Six*d1rw*f7*w5 & 1829 + d1rf7*w6 + Seven*d1rw*f8*w6 & 1830 + d1rf8*w7 + Eight*d1rw*f9*w7 + d1rf9*w8 1831 ! 1832 ! The final value of term1 for 0 < omega < wcutoff is: 1833 ! 1834 term1 = t1 + t2t9 + t10 1835 ! 1836 d1sterm1 = d1st1 + d1st2t9 + d1st10 1837 d1rterm1 = d1rt1 + d1rt2t9 + d1rt10 1838 ! 1839 ! The final value for the enhancement factor and its 1840 ! derivatives is: 1841 ! 1842 Fx_wpbe = X * (term1 + term2 + term3 + term4 + term5) 1843 ! 1844 d1sfx = X * (d1sterm1 + d1sterm2 + d1sterm3 & 1845 + d1sterm4 + d1sterm5) 1846 ! 1847 d1rfx = X * (d1rterm1 + d1rterm3 + d1rterm4 + d1rterm5) 1848 ! 1849 ENDIF 1850 1851END SUBROUTINE wpbe_analy_erfc_approx_grad 1852! 1853!--------------------------------------------------------------------- 1854function qe_erf(x) !<GPU:DEVICE> 1855 !--------------------------------------------------------------------- 1856 ! Error function - computed from the rational approximations of 1857 ! W. J. Cody, Math. Comp. 22 (1969), pages 631-637. 1858 ! 1859 ! for abs(x) le 0.47 erf is calculated directly 1860 ! for abs(x) gt 0.47 erf is calculated via erf(x)=1-erfc(x) 1861 USE kinds, ONLY: DP 1862 implicit none 1863 REAL(DP), intent(in) :: x 1864 REAL(DP) :: x2, p1 (4), q1 (4) 1865 REAL(DP) :: qe_erf !<GPU:qe_erf=>qe_erf_d> 1866 data p1 / 2.426679552305318E2, 2.197926161829415E1, & 1867 6.996383488619136d0, -3.560984370181538E-2 / 1868 data q1 / 2.150588758698612E2, 9.116490540451490E1, & 1869 1.508279763040779E1, 1.000000000000000d0 / 1870 ! 1871 if (abs (x) > 6.0d0) then 1872 ! 1873 ! erf(6)=1-10^(-17) cannot be distinguished from 1 1874 ! 1875 qe_erf = sign (1.0d0, x) !<GPU:qe_erf=>qe_erf_d> 1876 else 1877 if (abs (x) <= 0.47d0) then 1878 x2 = x**2 1879 qe_erf=x *(p1 (1) + x2 * (p1 (2) + x2 * (p1 (3) + x2 * p1 (4) ) ) ) & !<GPU:qe_erf=>qe_erf_d> 1880 / (q1 (1) + x2 * (q1 (2) + x2 * (q1 (3) + x2 * q1 (4) ) ) ) 1881 else 1882 qe_erf = 1.0d0 - qe_erfc(x) !<GPU:qe_erf=>qe_erf_d,qe_erfc=>qe_erfc_d> 1883 endif 1884 endif 1885 ! 1886 return 1887end function qe_erf 1888! 1889!--------------------------------------------------------------------- 1890function qe_erfc(x) !<GPU:DEVICE> 1891 !--------------------------------------------------------------------- 1892 ! 1893 ! erfc(x) = 1-erf(x) - See comments in erf 1894 ! 1895 USE kinds, ONLY: DP 1896 implicit none 1897 ! 1898 REAL(DP),intent(in) :: x 1899 REAL(DP) :: qe_erfc !<GPU:qe_erfc=>qe_erfc_d> 1900 REAL(DP) :: ax, x2, xm2, p2 (8), q2 (8), p3 (5), q3 (5), pim1 1901 ! 1902 data p2 / 3.004592610201616E2, 4.519189537118719E2, & 1903 3.393208167343437E2, 1.529892850469404E2, & 1904 4.316222722205674E1, 7.211758250883094d0, & 1905 5.641955174789740E-1,-1.368648573827167E-7 / 1906 data q2 / 3.004592609569833E2, 7.909509253278980E2, & 1907 9.313540948506096E2, 6.389802644656312E2, & 1908 2.775854447439876E2, 7.700015293522947E1, & 1909 1.278272731962942E1, 1.000000000000000d0 / 1910 data p3 /-2.996107077035422E-3,-4.947309106232507E-2, & 1911 -2.269565935396869E-1,-2.786613086096478E-1, & 1912 -2.231924597341847E-2 / 1913 data q3 / 1.062092305284679E-2, 1.913089261078298E-1, & 1914 1.051675107067932d0, 1.987332018171353d0, & 1915 1.000000000000000d0 / 1916 1917 data pim1 / 0.56418958354775629d0 / 1918 ! ( pim1= sqrt(1/pi) ) 1919 ax = abs (x) 1920 if (ax > 26.0d0) then 1921 ! 1922 ! erfc(26.0)=10^(-296); erfc( 9.0)=10^(-37); 1923 ! 1924 qe_erfc = 0.0d0 !<GPU:qe_erfc=>qe_erfc_d> 1925 elseif (ax > 4.0d0) then 1926 x2 = x**2 1927 xm2 = (1.0d0 / ax) **2 1928 qe_erfc = (1.0d0 / ax) * exp ( - x2) * (pim1 + xm2 * (p3 (1) & !<GPU:qe_erfc=>qe_erfc_d> 1929 + xm2 * (p3 (2) + xm2 * (p3 (3) + xm2 * (p3 (4) + xm2 * p3 (5) & 1930 ) ) ) ) / (q3 (1) + xm2 * (q3 (2) + xm2 * (q3 (3) + xm2 * & 1931 (q3 (4) + xm2 * q3 (5) ) ) ) ) ) 1932 elseif (ax > 0.47d0) then 1933 x2 = x**2 1934 qe_erfc = exp ( - x2) * (p2 (1) + ax * (p2 (2) + ax * (p2 (3) & !<GPU:qe_erfc=>qe_erfc_d> 1935 + ax * (p2 (4) + ax * (p2 (5) + ax * (p2 (6) + ax * (p2 (7) & 1936 + ax * p2 (8) ) ) ) ) ) ) ) / (q2 (1) + ax * (q2 (2) + ax * & 1937 (q2 (3) + ax * (q2 (4) + ax * (q2 (5) + ax * (q2 (6) + ax * & 1938 (q2 (7) + ax * q2 (8) ) ) ) ) ) ) ) 1939 else 1940 qe_erfc = 1.0d0 - qe_erf(ax) !<GPU:qe_erfc=>qe_erfc_d, qe_erf=>qe_erf_d> 1941 endif 1942 ! 1943 ! erf(-x)=-erf(x) => erfc(-x) = 2-erfc(x) 1944 ! 1945 if (x < 0.0d0) qe_erfc = 2.0d0 - qe_erfc !<GPU:qe_erfc=>qe_erfc_d> 1946 ! 1947 return 1948end function qe_erfc 1949! 1950 1951!------------------------------------------------------------------ 1952FUNCTION EXPINT(n, x) !<GPU:DEVICE> 1953!----------------------------------------------------------------------- 1954! Evaluates the exponential integral E_n(x) 1955! Parameters: maxit is the maximum allowed number of iterations, 1956! eps is the desired relative error, not smaller than the machine precision, 1957! big is a number near the largest representable floating-point number, 1958! Inspired from Numerical Recipes 1959! 1960 USE kinds, ONLY: DP 1961 IMPLICIT NONE 1962 INTEGER, INTENT(IN) :: n 1963 REAL(DP), INTENT(IN) :: x 1964 REAL(DP) :: expint !<GPU:expint=>expint_d> 1965 INTEGER, parameter :: maxit=200 1966 REAL(DP), parameter :: eps=1E-12, big=huge(x)*eps 1967 REAL(DP), parameter :: euler = 0.577215664901532860606512d0 1968! EPS=1E-9, FPMIN=1E-30 1969 1970 INTEGER :: i, nm1, k 1971 REAL(DP) :: a,b,c,d,del,fact,h,iarsum 1972 1973 IF (.NOT. ((n >= 0).AND.(x >= 0.0).AND.((x > 0.0).OR.(n > 1)))) THEN 1974 !CALL errore('expint','bad arguments', 1) 1975 STOP 1976 END IF 1977 1978 IF (n == 0) THEN 1979 expint = exp(-x)/x !<GPU:expint=>expint_d> 1980 RETURN 1981 END IF 1982 nm1 = n-1 1983 IF (x == 0.0d0) THEN 1984 expint = 1.0d0/nm1 !<GPU:expint=>expint_d> 1985 ELSE IF (x > 1.0d0) THEN 1986 b = x+n 1987 c = big 1988 d = 1.0d0/b 1989 h = d 1990 DO i=1,maxit 1991 a = -i*(nm1+i) 1992 b = b+2.0d0 1993 d = 1.0d0/(a*d+b) 1994 c = b+a/c 1995 del = c*d 1996 h = h*del 1997 IF (ABS(del-1.0d0) <= EPS) EXIT 1998 END DO 1999 IF (i > maxit) STOP !CALL errore('expint','continued fraction failed',1) 2000 expint = h*EXP(-x) !<GPU:expint=>expint_d> 2001 ELSE 2002 IF (nm1 /= 0) THEN 2003 expint = 1.0d0/nm1 !<GPU:expint=>expint_d> 2004 ELSE 2005 expint = -LOG(x)-euler !<GPU:expint=>expint_d> 2006 END IF 2007 fact = 1.0d0 2008 do i=1,maxit 2009 fact = -fact*x/i 2010 IF (i /= nm1) THEN 2011 del = -fact/(i-nm1) 2012 ELSE 2013 2014 iarsum = 0.0d0 2015 do k=1,nm1 2016 iarsum = iarsum + 1.0d0/k 2017 end do 2018 2019 del = fact*(-LOG(x)-euler+iarsum) 2020! del = fact*(-LOG(x)-euler+sum(1.0d0/arth(1,1,nm1))) 2021 END IF 2022 expint = expint + del !<GPU:expint=>expint_d> 2023 IF (ABS(del) < ABS(expint)*eps) EXIT !<GPU:expint=>expint_d> 2024 END DO 2025 IF (i > maxit) STOP !CALL errore('expint','series failed',1) 2026 END IF 2027END FUNCTION EXPINT 2028! 2029END MODULE 2030 2031