1! 2MODULE corr_gga !<GPU:corr_gga=>corr_gga_gpu> 3 ! 4 USE corr_lda, ONLY : pw, pw_spin !<GPU:pw_spin=>pw_spin_d,pw=>pw_d,corr_lda=>corr_lda_gpu> 5 ! 6CONTAINS 7! 8!----------------------------------------------------------------------- 9SUBROUTINE perdew86( rho, grho, sc, v1c, v2c ) !<GPU:DEVICE> 10 !----------------------------------------------------------------------- 11 !! Perdew gradient correction on correlation: PRB 33, 8822 (1986). 12 ! 13 USE kinds, ONLY : DP 14 ! 15 IMPLICIT NONE 16 ! 17 REAL(DP), INTENT(IN) :: rho, grho 18 REAL(DP), INTENT(OUT) :: sc, v1c, v2c 19 ! 20 ! ... local variables 21 ! 22 REAL(DP), PARAMETER :: p1=0.023266_DP, p2=7.389d-6, p3=8.723_DP, & 23 p4=0.472_DP 24 REAL(DP), PARAMETER :: pc1=0.001667_DP, pc2=0.002568_DP, pci=pc1 + pc2 25 REAL(DP), PARAMETER :: third=1._DP/3._DP, pi34=0.6203504908994_DP 26 ! pi34=(3/4pi)^(1/3) 27 REAL(DP) :: rho13, rho43, rs, rs2, rs3, cna, cnb, cn, drs 28 REAL(DP) :: dcna, dcnb, dcn, phi, ephi 29 ! 30 rho13 = rho**third 31 rho43 = rho13**4 32 rs = pi34 / rho13 33 rs2 = rs * rs 34 rs3 = rs * rs2 35 ! 36 cna = pc2 + p1 * rs + p2 * rs2 37 cnb = 1._DP + p3 * rs + p4 * rs2 + 1.d4 * p2 * rs3 38 cn = pc1 + cna / cnb 39 ! 40 drs = - third * pi34 / rho43 41 dcna = (p1 + 2._DP * p2 * rs) * drs 42 dcnb = (p3 + 2._DP * p4 * rs + 3.d4 * p2 * rs2) * drs 43 dcn = dcna / cnb - cna / (cnb * cnb) * dcnb 44 ! 45 phi = 0.192_DP * pci / cn * SQRT(grho) * rho**(-7._DP/6._DP) 46 ! SdG: in the original paper 1.745*0.11=0.19195 is used 47 ephi = EXP( -phi ) 48 ! 49 sc = grho / rho43 * cn * ephi 50 v1c = sc * ( (1._DP+phi) * dcn / cn - ((4._DP/3._DP)-(7._DP/ & 51 6._DP)*phi) / rho ) 52 v2c = cn * ephi / rho43 * (2._DP - phi) 53 ! 54 RETURN 55 ! 56END SUBROUTINE perdew86 57! 58! 59!----------------------------------------------------------------------- 60SUBROUTINE ggac( rho, grho, sc, v1c, v2c ) !<GPU:DEVICE> 61 !----------------------------------------------------------------------- 62 !! Perdew-Wang GGA (PW91) correlation part 63 ! 64 USE kinds, ONLY: DP 65 ! 66 IMPLICIT NONE 67 ! 68 REAL(DP), INTENT(IN) :: rho, grho 69 REAL(DP), INTENT(OUT) :: sc, v1c, v2c 70 ! 71 ! ... local variables 72 ! 73 REAL(DP) :: rs, ec, vc 74 ! 75 REAL(DP), PARAMETER :: al=0.09_DP, pa=0.023266_DP, pb=7.389d-6, & 76 pc=8.723_DP, pd=0.472_DP, & 77 cx=-0.001667_DP, cxc0=0.002568_DP, cc0=-cx+cxc0 78 ! 79 REAL(DP), PARAMETER :: third=1._DP/3._DP, pi34=0.6203504908994_DP, & 80 nu=15.755920349483144_DP, be=nu*cc0, & 81 xkf=1.919158292677513_DP, xks=1.128379167095513_DP 82 ! pi34=(3/4pi)^(1/3), nu=(16/pi)*(3 pi^2)^(1/3) 83 ! xkf=(9 pi/4)^(1/3), xks= sqrt(4/pi) 84 REAL(DP) :: kf, ks, rs1, rs2, rs3, t, expe, af, bf, y, xy, qy, s1 85 REAL(DP) :: h0, dh0, ddh0, ee, cn, dcn, cna, dcna, cnb, dcnb, h1, & 86 dh1, ddh1 87 ! 88 rs = pi34 / rho**third 89 ! 90 CALL pw( rs, 1, ec, vc ) !<GPU:pw=>pw_d> 91 ! 92 rs1 = rs 93 rs2 = rs1 * rs1 94 rs3 = rs1 * rs2 95 ! 96 kf = xkf / rs1 97 ks = xks * SQRT(kf) 98 t = SQRT(grho) / (2._DP * ks * rho) 99 ! 100 expe = EXP( - 2._DP * al * ec / (be * be) ) 101 af = 2._DP * al / be * (1._DP / (expe-1._DP) ) 102 bf = expe * (vc - ec) 103 ! 104 y = af * t * t 105 xy = (1._DP + y) / (1._DP + y + y * y) 106 qy = y * y * (2._DP + y) / (1._DP + y + y * y)**2 107 s1 = 1._DP + 2._DP * al / be * t * t * xy 108 ! 109 h0 = be * be / (2._DP * al) * LOG(s1) 110 dh0 = be * t * t / s1 * ( - 7._DP / 3._DP * xy - qy * (af * bf / & 111 be-7._DP / 3._DP) ) 112 ddh0 = be / (2._DP * ks * ks * rho) * (xy - qy) / s1 113 ! 114 ee = - 100._DP * (ks / kf * t)**2 115 ! 116 cna = cxc0 + pa * rs1 + pb * rs2 117 dcna = pa * rs1 + 2._DP * pb * rs2 118 cnb = 1._DP + pc * rs1 + pd * rs2 + 1.d4 * pb * rs3 119 dcnb = pc * rs1 + 2._DP * pd * rs2 + 3.d4 * pb * rs3 120 cn = cna / cnb - cx 121 dcn = dcna / cnb - cna * dcnb / (cnb * cnb) 122 ! 123 h1 = nu * (cn - cc0 - 3._DP / 7._DP * cx) * t * t * EXP(ee) 124 dh1 = - third * ( h1 * (7._DP + 8._DP * ee) + nu * t * t * EXP(ee) & 125 * dcn ) 126 ddh1 = 2._DP * h1 * (1._DP + ee) * rho / grho 127 ! 128 sc = rho * (h0 + h1) 129 v1c = h0 + h1 + dh0 + dh1 130 v2c = ddh0 + ddh1 131 ! 132 RETURN 133 ! 134END SUBROUTINE ggac 135! 136! 137!----------------------------------------------------------------------- 138SUBROUTINE glyp( rho, grho, sc, v1c, v2c ) !<GPU:DEVICE> 139 !----------------------------------------------------------------------- 140 !! Lee Yang Parr: gradient correction part. 141 ! 142 USE kinds, ONLY: DP 143 ! 144 IMPLICIT NONE 145 ! 146 REAL(DP), INTENT(IN) :: rho, grho 147 REAL(DP), INTENT(OUT) :: sc, v1c, v2c 148 ! 149 ! ... local varibles 150 ! 151 REAL(DP), PARAMETER :: a=0.04918_DP, b=0.132_DP, c=0.2533_DP, & 152 d=0.349_DP 153 REAL(DP) :: rhom13, rhom43, rhom53, om, xl, ff, dom, dxl 154 ! 155 rhom13 = rho**(-1._DP/3._DP) 156 om = EXP(-c*rhom13) / (1._DP+d*rhom13) 157 xl = 1._DP + (7._DP/3._DP) * ( c*rhom13 + d * rhom13 / (1._DP + & 158 d * rhom13) ) 159 ff = a * b * grho / 24._DP 160 rhom53 = rhom13**5 161 ! 162 sc = ff * rhom53 * om * xl 163 ! 164 dom = - om * (c + d+c * d * rhom13) / (1.d0 + d * rhom13) 165 dxl = (7.d0 / 3.d0) * (c + d+2.d0 * c * d * rhom13 + c * d * d * & 166 rhom13**2) / (1.d0 + d * rhom13) **2 167 rhom43 = rhom13**4 168 ! 169 v1c = - ff * rhom43 / 3.d0 * ( 5.d0 * rhom43 * om * xl + rhom53 * & 170 dom * xl + rhom53 * om * dxl ) 171 v2c = 2.d0 * sc / grho 172 ! 173 RETURN 174 ! 175END SUBROUTINE glyp 176! 177! 178!--------------------------------------------------------------- 179SUBROUTINE pbec( rho, grho, iflag, sc, v1c, v2c ) !<GPU:DEVICE> 180 !--------------------------------------------------------------- 181 !! PBE correlation (without LDA part) 182 ! 183 !! * iflag=1: J.P.Perdew, K.Burke, M.Ernzerhof, PRL 77, 3865 (1996). 184 !! * iflag=2: J.P.Perdew et al., PRL 100, 136406 (2008). 185 !! * iflag=3: L. Chiodo et al, PRL 108, 126402 (2012) (PBEQ2D) 186 ! 187 USE kinds, ONLY: DP 188 ! 189 IMPLICIT NONE 190 ! 191 INTEGER, INTENT(IN) :: iflag !<GPU:VALUE> 192 REAL(DP), INTENT(IN) :: rho, grho 193 ! input: charge and squared gradient 194 REAL(DP), INTENT(OUT) :: sc, v1c, v2c 195 ! output: energy, potential 196 REAL(DP), PARAMETER :: ga = 0.0310906908696548950_DP 197 REAL(DP) :: be(3) 198 ! pbe pbesol pbeq2d 199 DATA be / 0.06672455060314922_DP, 0.046_DP, 0.06672455060314922_DP/ 200 REAL(DP), PARAMETER :: third = 1.d0 / 3.d0, pi34 = 0.6203504908994d0 201 REAL(DP), PARAMETER :: xkf = 1.919158292677513d0, xks = 1.128379167095513d0 202 ! pi34=(3/4pi)^(1/3), xkf=(9 pi/4)^(1/3), xks= sqrt(4/pi) 203 ! 204 REAL(DP) :: rs, ec, vc 205 ! 206 REAL(DP) :: kf, ks, t, expe, af, bf, y, xy, qy 207 REAL(DP) :: s1, h0, dh0, ddh0, sc2D, v1c2D, v2c2D 208 ! 209 rs = pi34 / rho**third 210 ! 211 CALL pw( rs, 1, ec, vc ) !<GPU:pw=>pw_d> 212 ! 213 kf = xkf / rs 214 ks = xks * SQRT(kf) 215 t = SQRT(grho) / (2._DP * ks * rho) 216 expe = EXP( - ec / ga ) 217 af = be(iflag) / ga * (1._DP / (expe-1._DP)) 218 bf = expe * (vc - ec) 219 y = af * t * t 220 xy = (1._DP + y) / (1._DP + y + y * y) 221 qy = y * y * (2._DP + y) / (1._DP + y + y * y)**2 222 s1 = 1._DP + be(iflag) / ga * t * t * xy 223 h0 = ga * LOG(s1) 224 dh0 = be(iflag) * t * t / s1 * ( - 7._DP / 3._DP * xy - qy * (af * bf / & 225 be(iflag)-7._DP / 3._DP) ) 226 ddh0 = be(iflag) / (2._DP * ks * ks * rho) * (xy - qy) / s1 227 ! 228 sc = rho * h0 229 v1c = h0 + dh0 230 v2c = ddh0 231 ! q2D 232 IF (iflag == 3) THEN 233 CALL cpbe2d( rho, grho, sc2D, v1c2D, v2c2D ) !<GPU:cpbe2d=>cpbe2d_d> 234 sc = sc + sc2D 235 v1c = v1c + v1c2D 236 v2c = v2c + v2c2D 237 ENDIF 238 ! 239 RETURN 240 ! 241END SUBROUTINE pbec 242! 243! 244! ===========> SPIN <=========== 245! 246!----------------------------------------------------------------------- 247SUBROUTINE perdew86_spin( rho, zeta, grho, sc, v1c_up, v1c_dw, v2c ) !<GPU:DEVICE> 248 !--------------------------------------------------------------------- 249 !! Perdew gradient correction on correlation: PRB 33, 8822 (1986) 250 !! spin-polarized case. 251 ! 252 USE kinds, ONLY: DP 253 ! 254 IMPLICIT NONE 255 ! 256 REAL(DP), INTENT(IN) :: rho 257 !! the total charge density 258 REAL(DP), INTENT(IN) :: zeta 259 !! the magnetization 260 REAL(DP), INTENT(IN) :: grho 261 !! the gradient of the charge squared 262 REAL(DP), INTENT(OUT) :: sc 263 !! correlation energies 264 REAL(DP), INTENT(OUT) :: v1c_up, v1c_dw !v1c(2) 265 !! derivative of correlation wr. rho 266 REAL(DP), INTENT(OUT) :: v2c 267 !! derivatives of correlation wr. grho 268 ! 269 ! ... local variables 270 ! 271 REAL(DP), PARAMETER :: p1=0.023266_DP, p2=7.389D-6, p3=8.723_DP, & 272 p4=0.472_DP 273 REAL(DP), PARAMETER :: pc1=0.001667_DP, pc2 = 0.002568_DP, pci=pc1+pc2 274 REAL(DP), PARAMETER :: third=1._DP/3._DP, pi34=0.6203504908994_DP 275 ! pi34=(3/4pi)^(1/3) 276 ! 277 REAL(DP) :: rho13, rho43, rs, rs2, rs3, cna, cnb, cn, drs 278 REAL(DP) :: dcna, dcnb, dcn, phi, ephi, dd, ddd 279 ! 280 rho13 = rho**third 281 rho43 = rho13**4 282 rs = pi34 / rho13 283 rs2 = rs * rs 284 rs3 = rs * rs2 285 ! 286 cna = pc2 + p1 * rs + p2 * rs2 287 cnb = 1._DP + p3 * rs + p4 * rs2 + 1.D4 * p2 * rs3 288 cn = pc1 + cna / cnb 289 ! 290 drs = -third * pi34 / rho43 291 dcna = (p1 + 2._DP * p2 * rs) * drs 292 dcnb = (p3 + 2._DP * p4 * rs + 3.D4 * p2 * rs2) * drs 293 dcn = dcna / cnb - cna / (cnb * cnb) * dcnb 294 phi = 0.192_DP * pci / cn * SQRT(grho) * rho**(-7._DP/6._DP) 295 !SdG: in the original paper 1.745*0.11=0.19195 is used 296 ! 297 dd = (2._DP)**third * SQRT( ( (1.0_DP + zeta) * 0.5_DP)**(5.0_DP/ & 298 3._DP) + ( (1.0_DP - zeta) * 0.5d0)**(5._DP/3._DP) ) 299 ! 300 ddd = (2._DP)**(-4.0_DP/3._DP) * 5._DP * ( ((1._DP + zeta) & 301 * 0.5_DP)**(2._DP/3._DP) - ((1._DP - zeta)*0.5d0)**(2._DP/ & 302 3._DP))/(3._DP*dd) 303 ! 304 ephi = EXP(-phi) 305 sc = grho / rho43 * cn * ephi / dd 306 ! 307 v1c_up = sc * ( (1._DP + phi) * dcn / cn - ( (4._DP / 3._DP) - & 308 (7._DP/6._DP)*phi)/rho) - sc * ddd / dd * (1._DP - zeta)/rho 309 v1c_dw = sc * ( (1._DP + phi) * dcn / cn - ( (4._DP / 3._DP) - & 310 (7._DP/6._DP)*phi)/rho) + sc * ddd / dd * (1._DP + zeta)/rho 311 v2c = cn * ephi / rho43 * (2._DP - phi) / dd 312 ! 313 RETURN 314 ! 315END SUBROUTINE perdew86_spin 316! 317! 318!----------------------------------------------------------------------- 319SUBROUTINE ggac_spin( rho, zeta, grho, sc, v1c_up, v1c_dw, v2c ) !<GPU:DEVICE> 320 !--------------------------------------------------------------------- 321 !! Perdew-Wang GGA (PW91) correlation part - spin-polarized 322 ! 323 USE kinds, ONLY: DP 324 ! 325 IMPLICIT NONE 326 ! 327 REAL(DP), INTENT(IN) :: rho 328 !! the total charge density 329 REAL(DP), INTENT(IN) :: zeta 330 !! the magnetization 331 REAL(DP), INTENT(IN) :: grho 332 !! the gradient of the charge squared 333 REAL(DP), INTENT(OUT) :: sc 334 !! correlation energies 335 REAL(DP), INTENT(OUT) :: v1c_up, v1c_dw 336 !! derivative of correlation wr. rho 337 REAL(DP), INTENT(OUT) :: v2c 338 !! derivatives of correlation wr. grho 339 ! 340 ! ... local variables 341 ! 342 REAL(DP), PARAMETER :: al=0.09_DP, pa=0.023266_DP, pb=7.389D-6, & 343 pc=8.723_DP, pd=0.472_DP 344 REAL(DP), PARAMETER :: cx=-0.001667_DP, cxc0=0.002568_DP, cc0=-cx+cxc0 345 REAL(DP), PARAMETER :: third=1._DP/3._DP, pi34=0.6203504908994_DP 346 ! pi34=(3/4pi)^(1/3) 347 REAL(DP), PARAMETER :: nu=15.755920349483144_DP, be=nu*cc0 348 ! nu=(16/pi)*(3 pi^2)^(1/3) 349 REAL(DP), PARAMETER :: xkf=1.919158292677513_DP, xks=1.128379167095513_DP 350 ! xkf=(9 pi/4)^(1/3) , xks=sqrt(4/pi) 351 ! 352 REAL(DP) :: rs, ec, vc_up, vc_dn 353 REAL(DP) :: kf, ks, rs2, rs3, t, expe, af, y, xy, qy, s1, h0, ddh0, ee, & 354 cn, dcn, cna, dcna, cnb, dcnb, h1, dh1, ddh1, fz, fz2, fz3, & 355 fz4, dfz, bfup, bfdw, dh0up, dh0dw, dh0zup, dh0zdw, dh1zup, & 356 dh1zdw 357 ! 358 rs = pi34 / rho**third 359 ! 360 CALL pw_spin( rs, zeta, ec, vc_up, vc_dn ) !<GPU:pw_spin=>pw_spin_d> 361 ! 362 rs2 = rs * rs 363 rs3 = rs * rs2 364 kf = xkf / rs 365 ks = xks * SQRT(kf) 366 ! 367 fz = 0.5_DP * ( (1._DP+zeta)**(2._DP/3._DP) + (1._DP-zeta)**(2._DP/3._DP) ) 368 fz2 = fz * fz 369 fz3 = fz2 * fz 370 fz4 = fz3 * fz 371 dfz = ( (1._DP+zeta)**(-1._DP/3._DP) - (1._DP-zeta)**(-1._DP/3._DP) )/3._DP 372 ! 373 t = SQRT(grho) / (2._DP * fz * ks * rho) 374 expe = EXP( - 2._DP * al * ec / (fz3 * be * be) ) 375 af = 2._DP * al / be * (1._DP / (expe-1._DP) ) 376 bfup = expe * (vc_up - ec) / fz3 377 bfdw = expe * (vc_dn - ec) / fz3 378 ! 379 y = af * t * t 380 xy = (1._DP + y) / (1._DP + y + y * y) 381 qy = y * y * (2._DP + y) / (1._DP + y + y * y)**2 382 s1 = 1._DP + 2._DP * al / be * t * t * xy 383 ! 384 h0 = fz3 * be * be / (2._DP * al) * LOG(s1) 385 dh0up = be * t * t * fz3 / s1 * ( - 7._DP / 3._DP * xy - qy * & 386 (af * bfup / be-7._DP / 3._DP) ) 387 dh0dw = be * t * t * fz3 / s1 * ( - 7._DP / 3._DP * xy - qy * & 388 (af * bfdw / be-7._DP / 3._DP) ) 389 dh0zup = (3._DP * h0 / fz - be * t * t * fz2 / s1 * (2._DP *xy - & 390 qy * (3._DP * af * expe * ec / fz3 / be+2._DP) ) ) * & 391 dfz * (1._DP - zeta ) 392 dh0zdw = -(3._DP * h0 / fz - be * t * t * fz3 / s1 * (2._DP*xy - & 393 qy * (3._DP * af * expe * ec / fz3 / be + 2._DP) ) ) * & 394 dfz * (1._DP + zeta ) 395 ddh0 = be * fz / (2._DP * ks * ks * rho) * (xy - qy) / s1 396 ! 397 ee = -100._DP * fz4 * (ks / kf * t)**2 398 cna = cxc0 + pa * rs + pb * rs2 399 dcna = pa * rs + 2._DP * pb * rs2 400 cnb = 1._DP + pc * rs + pd * rs2 + 1.D4 * pb * rs3 401 dcnb = pc * rs + 2._DP * pd * rs2 + 3.D4 * pb * rs3 402 cn = cna / cnb - cx 403 dcn = dcna / cnb - cna * dcnb / (cnb * cnb) 404 h1 = nu * (cn - cc0 - 3._DP / 7._DP * cx) * fz3 * t * t * EXP(ee) 405 dh1 = - third * (h1 * (7._DP + 8._DP * ee) + fz3 * nu * t * t * & 406 EXP(ee) * dcn) 407 ! 408 ddh1 = 2._DP * h1 * (1._DP + ee) * rho / grho 409 dh1zup = (1._DP - zeta) * dfz * h1 * (1._DP + 2._DP * ee / fz) 410 dh1zdw = -(1._DP + zeta) * dfz * h1 * (1._DP + 2._DP * ee / fz) 411 ! 412 sc = rho * (h0 + h1) 413 v1c_up = h0 + h1 + dh0up + dh1 + dh0zup + dh1zup 414 v1c_dw = h0 + h1 + dh0up + dh1 + dh0zdw + dh1zdw 415 v2c = ddh0 + ddh1 416 ! 417 ! 418 RETURN 419 ! 420END SUBROUTINE ggac_spin 421! 422! 423!------------------------------------------------------------------- 424SUBROUTINE pbec_spin( rho, zeta, grho, iflag, sc, v1c_up, v1c_dw, v2c ) !<GPU:DEVICE> 425 !----------------------------------------------------------------- 426 !! PBE correlation (without LDA part) - spin-polarized. 427 ! 428 !! * iflag = 1: J.P.Perdew, K.Burke, M.Ernzerhof, PRL 77, 3865 (1996); 429 !! * iflag = 2: J.P.Perdew et al., PRL 100, 136406 (2008). 430 ! 431 USE kinds, ONLY : DP 432 ! 433 IMPLICIT NONE 434 ! 435 INTEGER, INTENT(IN) :: iflag !<GPU:VALUE> 436 !! see main comments 437 REAL(DP), INTENT(IN) :: rho 438 !! the total charge density 439 REAL(DP), INTENT(IN) :: zeta 440 !! the magnetization 441 REAL(DP), INTENT(IN) :: grho 442 !! the gradient of the charge squared 443 REAL(DP), INTENT(OUT) :: sc 444 !! correlation energies 445 REAL(DP), INTENT(OUT) :: v1c_up, v1c_dw 446 !! derivative of correlation wr. rho 447 REAL(DP), INTENT(OUT) :: v2c 448 !! derivatives of correlation wr. grho 449 ! 450 ! ... local variables 451 ! 452 REAL(DP) :: rs, ec, vc_up, vc_dn 453 ! 454 REAL(DP), PARAMETER :: ga=0.031091_DP 455 REAL(DP) :: be(2) 456 DATA be / 0.06672455060314922_DP, 0.046_DP / 457 REAL(DP), PARAMETER :: third=1.D0/3.D0, pi34=0.6203504908994_DP 458 ! pi34=(3/4pi)^(1/3) 459 REAL(DP), PARAMETER :: xkf=1.919158292677513_DP, xks=1.128379167095513_DP 460 ! xkf=(9 pi/4)^(1/3) , xks=sqrt(4/pi) 461 REAL(DP) :: kf, ks, t, expe, af, y, xy, qy, s1, h0, ddh0 462 REAL(DP) :: fz, fz2, fz3, fz4, dfz, bfup, bfdw, dh0up, dh0dw, & 463 dh0zup, dh0zdw 464 ! 465 ! 466 rs = pi34 / rho**third 467 ! 468 CALL pw_spin( rs, zeta, ec, vc_up, vc_dn ) !<GPU:pw_spin=>pw_spin_d> 469 ! 470 kf = xkf / rs 471 ks = xks * SQRT(kf) 472 ! 473 fz = 0.5_DP*( (1._DP+zeta)**(2._DP/3._DP) + (1._DP-zeta)**(2._DP/3._DP) ) 474 fz2 = fz * fz 475 fz3 = fz2 * fz 476 fz4 = fz3 * fz 477 dfz = ( (1._DP+zeta)**(-1._DP/3._DP) - (1._DP - zeta)**(-1._DP/3._DP) ) & 478 / 3._DP 479 ! 480 t = SQRT(grho) / (2._DP * fz * ks * rho) 481 expe = EXP( - ec / (fz3 * ga) ) 482 af = be(iflag) / ga * (1._DP / (expe-1._DP) ) 483 bfup = expe * (vc_up - ec) / fz3 484 bfdw = expe * (vc_dn - ec) / fz3 485 ! 486 y = af * t * t 487 xy = (1._DP + y) / (1._DP + y + y * y) 488 qy = y * y * (2._DP + y) / (1._DP + y + y * y)**2 489 s1 = 1._DP + be(iflag) / ga * t * t * xy 490 ! 491 h0 = fz3 * ga * LOG(s1) 492 ! 493 dh0up = be(iflag) * t * t * fz3 / s1 * ( -7._DP/3._DP * xy - qy * & 494 (af * bfup / be(iflag)-7._DP/3._DP) ) 495 ! 496 dh0dw = be(iflag) * t * t * fz3 / s1 * ( -7._DP/3._DP * xy - qy * & 497 (af * bfdw / be(iflag)-7._DP/3._DP) ) 498 ! 499 dh0zup = (3._DP * h0 / fz - be(iflag) * t * t * fz2 / s1 * & 500 (2._DP * xy - qy * (3._DP * af * expe * ec / fz3 / & 501 be(iflag)+2._DP) ) ) * dfz * (1._DP - zeta) 502 ! 503 dh0zdw = - (3._DP * h0 / fz - be(iflag) * t * t * fz2 / s1 * & 504 (2._DP * xy - qy * (3._DP * af * expe * ec / fz3 / & 505 be(iflag)+2._DP) ) ) * dfz * (1._DP + zeta) 506 ! 507 ddh0 = be(iflag) * fz / (2._DP * ks * ks * rho) * (xy - qy) / s1 508 ! 509 sc = rho * h0 510 v1c_up = h0 + dh0up + dh0zup 511 v1c_dw = h0 + dh0dw + dh0zdw 512 v2c = ddh0 513 ! 514 ! 515 RETURN 516 ! 517END SUBROUTINE pbec_spin 518! 519! 520!------------------------------------------------------------------------ 521SUBROUTINE lsd_glyp( rho_in_up, rho_in_dw, grho_up, grho_dw, grho_ud, sc, v1c_up, v1c_dw, v2c_up, v2c_dw, v2c_ud ) !<GPU:DEVICE> 522 !---------------------------------------------------------------------- 523 !! Lee, Yang, Parr: gradient correction part. 524 ! 525 USE kinds, ONLY: DP 526 ! 527 IMPLICIT NONE 528 ! 529 REAL(DP), INTENT(IN) :: rho_in_up, rho_in_dw 530 !! the total charge density 531 REAL(DP), INTENT(IN) :: grho_up, grho_dw 532 !! the gradient of the charge squared 533 REAL(DP), INTENT(IN) :: grho_ud 534 !! gradient off-diagonal term up-down 535 REAL(DP), INTENT(OUT) :: sc 536 !! correlation energy 537 REAL(DP), INTENT(OUT) :: v1c_up, v1c_dw 538 !! derivative of correlation wr. rho 539 REAL(DP), INTENT(OUT) :: v2c_up, v2c_dw 540 !! derivatives of correlation wr. grho 541 REAL(DP), INTENT(OUT) :: v2c_ud 542 !! derivative of correlation wr. grho, off-diag. term 543 ! 544 ! ... local variables 545 ! 546 REAL(DP) :: ra, rb, rho, grhoaa, grhoab, grhobb 547 REAL(DP) :: rm3, dr, or, dor, der, dder 548 REAL(DP) :: dlaa, dlab, dlbb, dlaaa, dlaab, dlaba, dlabb, dlbba, dlbbb 549 REAL(DP), PARAMETER :: a=0.04918_DP, b=0.132_DP, c=0.2533_DP, d=0.349_DP 550 ! 551 ra = rho_in_up 552 rb = rho_in_dw 553 rho = ra + rb 554 rm3 = rho**(-1._DP/3._DP) 555 ! 556 dr = ( 1._DP + d*rm3 ) 557 or = EXP(-c*rm3) / dr * rm3**11._DP 558 dor = -1._DP/3._DP * rm3**4 * or * (11._DP/rm3-c-d/dr) 559 ! 560 der = c*rm3 + d*rm3/dr 561 dder = 1._DP/3._DP * (d*d*rm3**5/dr/dr - der/rho) 562 ! 563 dlaa = -a*b*or*( ra*rb/9._DP*(1._DP-3*der-(der-11._DP)*ra/rho)-rb*rb ) 564 dlab = -a*b*or*( ra*rb/9._DP*(47._DP-7._DP*der)-4._DP/3._DP*rho*rho ) 565 dlbb = -a*b*or*( ra*rb/9._DP*(1._DP-3*der-(der-11._DP)*rb/rho)-ra*ra ) 566 ! 567 dlaaa = dor/or*dlaa-a*b*or*(rb/9._DP*(1._DP-3*der-(der-11._DP)*ra/rho)- & 568 ra*rb/9._DP*((3._DP+ra/rho)*dder+(der-11._DP)*rb/rho/rho)) 569 dlaab = dor/or*dlaa-a*b*or*(ra/9._DP*(1._DP-3._DP*der-(der-11._DP)*ra/rho)- & 570 ra*rb/9._DP*((3._DP+ra/rho)*dder-(der-11._DP)*ra/rho/rho)-2._DP*rb) 571 dlaba = dor/or*dlab-a*b*or*(rb/9._DP*(47._DP-7._DP*der)-7._DP/9._DP*ra*rb*dder- & 572 8._DP/3._DP*rho) 573 dlabb = dor/or*dlab-a*b*or*(ra/9._DP*(47._DP-7._DP*der)-7._DP/9._DP*ra*rb*dder- & 574 8._DP/3._DP*rho) 575 dlbba = dor/or*dlbb-a*b*or*(rb/9._DP*(1._DP-3._DP*der-(der-11._DP)*rb/rho)- & 576 ra*rb/9._DP*((3._DP+rb/rho)*dder-(der-11._DP)*rb/rho/rho)-2._DP*ra) 577 dlbbb = dor/or*dlbb-a*b*or*(ra/9._DP*(1._DP-3*der-(der-11._DP)*rb/rho)- & 578 ra*rb/9._DP*((3._DP+rb/rho)*dder+(der-11._DP)*ra/rho/rho)) 579 ! 580 grhoaa = grho_up 581 grhobb = grho_dw 582 grhoab = grho_ud 583 ! 584 sc = dlaa *grhoaa + dlab *grhoab + dlbb *grhobb 585 v1c_up = dlaaa*grhoaa + dlaba*grhoab + dlbba*grhobb 586 v1c_dw = dlaab*grhoaa + dlabb*grhoab + dlbbb*grhobb 587 v2c_up = 2._DP*dlaa 588 v2c_dw = 2._DP*dlbb 589 v2c_ud = dlab 590 ! 591 ! 592 RETURN 593 ! 594END SUBROUTINE lsd_glyp 595! 596! 597!--------------------------------------------------------------- 598SUBROUTINE cpbe2d( rho, grho, sc, v1c, v2c ) !<GPU:DEVICE> 599 !--------------------------------------------------------------- 600 !! 2D correction (last term of Eq. 5, PRL 108, 126402 (2012)) 601 ! 602 USE kinds, ONLY: DP 603 ! 604 IMPLICIT NONE 605 ! 606 REAL(DP), INTENT(IN) :: rho, grho 607 REAL(DP), INTENT(OUT) :: sc, v1c, v2c 608 ! 609 ! ... local variables 610 ! 611 REAL(8), PARAMETER :: pi=3.14159265358979323846d0 612 REAL(DP), PARAMETER :: ex1=0.333333333333333333_DP, ex2=1.166666666666667_DP 613 REAL(DP), PARAMETER :: ex3=ex2+1.0_DP 614 REAL(DP) :: fac1, fac2, zeta, phi, gr, rs, drsdn, akf, aks, t, dtdn, dtdgr 615 REAL(DP) :: p, a, g, alpha1, beta1,beta2,beta3,beta4, dgdrs, epsc, depscdrs 616 REAL(DP) :: c, gamma1, beta, aa, cg, adddepsc, h, dhdaa, dhdt, dhdrs 617 REAL(DP) :: epscpbe, depscpbedrs, depscpbedt, a0,a1,a2, b0,b1,b2, c0,c1,c2 618 REAL(DP) :: e0,e1,e2, f0,f1,f2, g0,g1,g2, h0,h1,h2, d0,d1,d2, ff, dffdt 619 REAL(DP) :: rs3d, rs2d, drs2ddrs3d, eps2d, deps2ddrs2, depsGGAdrs, depsGGAdt 620 REAL(DP) :: drs2ddt, rs2, ec, decdn, decdgr, daadepsc 621 ! 622 fac1 = (3.d0*pi*pi)**ex1 623 fac2 = SQRT(4.d0*fac1/pi) 624 ! 625 zeta = 0.d0 626 phi = 1.d0 627 ! 628 gr = SQRT(grho) 629 ! 630 rs = (3.d0/4.d0/pi/rho)**ex1 631 drsdn = -DBLE(3 ** (0.1D1 / 0.3D1)) * DBLE(2 ** (0.1D1 / 0.3D1)) * & 632 0.3141592654D1 ** (-0.1D1 / 0.3D1) * (0.1D1 / rho) ** (-0.2D1 / & 633 0.3D1) / rho ** 2 / 0.6D1 634 ! 635 akf = (3.d0*pi*pi*rho)**(1.d0/3.d0) 636 aks = DSQRT(4.d0*akf/pi) 637 t = gr/2.d0 / phi / aks / rho 638 dtdn = -7.d0/6.d0 * gr/2.d0 / phi/DSQRT(4.d0/pi)/ & 639 ((3.d0*pi*pi)**(1.d0/6.d0)) / (rho**(13.d0/6.d0)) 640 dtdgr = 1.d0/2.d0/phi/aks/rho 641 ! 642 ! for the LDA correlation 643 p = 1.d0 644 A = 0.0310906908696548950_DP 645 alpha1 = 0.21370d0 646 beta1 = 7.5957d0 647 beta2 = 3.5876d0 648 beta3 = 1.6382d0 649 beta4 = 0.49294d0 650 G = -0.2D1 * A * DBLE(1 + alpha1 * rs) * LOG(0.1D1 + 0.1D1 / A / ( & 651 beta1 * SQRT(DBLE(rs)) + DBLE(beta2 * rs) + DBLE(beta3 * rs ** ( & 652 0.3D1 / 0.2D1)) + DBLE(beta4 * rs ** (p + 1))) / 0.2D1) 653 ! 654 dGdrs = -0.2D1 * A * alpha1 * LOG(0.1D1 + 0.1D1 / A / (beta1 * SQRT(rs) & 655 + beta2 * rs + beta3 * rs ** (0.3D1 / 0.2D1) + beta4 * rs ** & 656 (p + 1)) / 0.2D1) + (0.1D1 + alpha1 * rs) / (beta1 * SQRT(rs) + & 657 beta2 * rs + beta3 * rs ** (0.3D1 / 0.2D1) + beta4 * rs ** (p + 1)) & 658 ** 2 * (beta1 * rs ** (-0.1D1 / 0.2D1) / 0.2D1 + beta2 + 0.3D1 / & 659 0.2D1 * beta3 * SQRT(rs) + beta4 * rs ** (p + 1) * DBLE(p + 1) / & 660 rs) / (0.1D1 + 0.1D1 / A / (beta1 * SQRT(rs) + beta2 * rs + beta3 * & 661 rs ** (0.3D1 / 0.2D1) + beta4 * rs ** (p + 1)) / 0.2D1) 662 ! 663 epsc = G 664 depscdrs = dGdrs 665 ! 666 ! PBE 667 c = 1.d0 668 gamma1 = 0.0310906908696548950_dp 669 beta = 0.06672455060314922_dp 670 ! 671 AA = beta / gamma1 / (EXP(-epsc / gamma1 / phi ** 3) - 0.1D1) 672 cg = beta / gamma1 ** 2 / (EXP(-epsc/ gamma1 / phi ** 3) - 0.1D1) & 673 ** 2 / phi ** 3 * EXP(-epsc / gamma1 / phi ** 3) 674 dAAdepsc = cg 675 ! 676 IF (t <= 10.d0) THEN 677 ! 678 H = DBLE(gamma1) * phi ** 3 * LOG(DBLE(1 + beta / gamma1 * t ** 2 & 679 * (1 + AA * t ** 2) / (1 + c * AA * t ** 2 + AA ** 2 * t ** 4))) 680 ! 681 dHdAA = gamma1 * phi ** 3 * (beta / gamma1 * t ** 4 / (1 + c * AA & 682 * t ** 2 + AA ** 2 * t ** 4) - beta / gamma1 * t ** 2 * (1 + AA * & 683 t ** 2) / (1 + c * AA * t ** 2 + AA ** 2 * t ** 4) ** 2 * (c * t **& 684 2 + 2 * AA * t ** 4)) / (1 + beta / gamma1 * t ** 2 * (1 + AA * & 685 t ** 2) / (1 + c * AA * t ** 2 + AA ** 2 * t ** 4)) 686 ! 687 dHdt = gamma1 * phi ** 3 * (2 * beta / gamma1 * t * (1 + AA * t ** & 688 2) / (1 + c * AA * t ** 2 + AA ** 2 * t ** 4) + 2 * beta / gamma1 & 689 * t ** 3 * AA / (1 + c * AA * t ** 2 + AA ** 2 * t ** 4) - beta / & 690 gamma1 * t ** 2 * (1 + AA * t ** 2) / (1 + c * AA * t ** 2 + AA ** & 691 2 * t ** 4) ** 2 * (2 * c * AA * t + 4 * AA ** 2 * t ** 3)) / (1 & 692 + beta / gamma1 * t ** 2 * (1 + AA * t ** 2) / (1 + c * AA * t ** & 693 2 + AA ** 2 * t ** 4)) 694 ! 695 ELSE 696 ! 697 H = gamma1 * (phi**3) * DLOG(1.d0+(beta/gamma1)*(1.d0/AA)) 698 ! 699 dHdAA = gamma1 * (phi**3) * 1.d0/(1.d0+(beta/gamma1)*(1.d0/AA))* & 700 (beta/gamma1) * (-1.d0/AA/AA) 701 ! 702 dHdt = 0.d0 703 ! 704 ENDIF 705 ! 706 dHdrs = dHdAA*dAAdepsc*depscdrs 707 ! 708 epscPBE = epsc+H 709 depscPBEdrs = depscdrs+dHdrs 710 depscPBEdt = dHdt 711 ! 712 ! START THE 2D CORRECTION 713 ! 714 beta = 1.3386d0 715 a0 = -0.1925d0 716 a1 = 0.117331d0 717 a2 = 0.0234188d0 718 b0 = 0.0863136d0 719 b1 = -0.03394d0 720 b2 = -0.037093d0 721 c0 = 0.057234d0 722 c1 = -0.00766765d0 723 c2 = 0.0163618d0 724 e0 = 1.0022d0 725 e1 = 0.4133d0 726 e2 = 1.424301d0 727 f0 = -0.02069d0 728 f1 = 0.d0 729 f2 = 0.d0 730 g0 = 0.340d0 731 g1 = 0.0668467d0 732 g2 = 0.d0 733 h0 = 0.01747d0 734 h1 = 0.0007799d0 735 h2 = 1.163099d0 736 d0 = -a0*h0 737 d1 = -a1*h1 738 d2 = -a2*h2 739 ! 740 ff = t ** 4 * (1 + t ** 2) / (1000000 + t ** 6) 741 dffdt = 4 * t ** 3 * (1 + t ** 2) / (1000000 + t ** 6) + 2 * t ** & 742 5 / (1000000 + t ** 6) - 6 * t ** 9 * (1 + t ** 2) / (1000000 + t & 743 ** 6) ** 2 744! 745 rs3d=rs 746 rs2d = 0.4552100000D0 * DBLE(3 ** (0.7D1 / 0.12D2)) * DBLE(4 ** ( & 747 0.5D1 / 0.12D2)) * (0.1D1 / pi) ** (-0.5D1 / 0.12D2) * rs3d ** ( & 748 0.5D1 / 0.4D1) * SQRT(t) 749 750 cg = 0.5690125000D0 * DBLE(3 ** (0.7D1 / 0.12D2)) * DBLE(4 ** ( & 751 0.5D1 / 0.12D2)) * (0.1D1 / pi) ** (-0.5D1 / 0.12D2) * rs3d ** (0.1D1 & 752 / 0.4D1) * SQRT(t) 753 drs2ddrs3d=cg 754 755 cg = 0.2276050000D0 * DBLE(3 ** (0.7D1 / 0.12D2)) * DBLE(4 ** ( & 756 0.5D1 / 0.12D2)) * DBLE((1 / pi) ** (-0.5D1 / 0.12D2)) * DBLE(rs3d ** & 757 (0.5D1 / 0.4D1)) * DBLE(t ** (-0.1D1 / 0.2D1)) 758 drs2ddt=cg 759 rs2=rs2d 760 ! 761 eps2d = (EXP(-beta * rs2) - 0.1D1) * (-0.2D1 / 0.3D1 * SQRT(0.2D1) & 762 * DBLE((1 + zeta) ** (0.3D1 / 0.2D1) + (1 - zeta) ** (0.3D1 / & 763 0.2D1)) / pi / rs2 + 0.4D1 / 0.3D1 * (0.1D1 + 0.3D1 / 0.8D1 * DBLE( & 764 zeta ** 2) + 0.3D1 / 0.128D3 * DBLE(zeta ** 4)) * SQRT(0.2D1) / pi / & 765 rs2) + a0 + (b0 * rs2 + c0 * rs2 ** 2 + d0 * rs2 ** 3) * LOG(0.1D1 & 766 + 0.1D1 / (e0 * rs2 + f0 * rs2 ** (0.3D1 / 0.2D1) + g0 * rs2 ** & 767 2 + h0 * rs2 ** 3)) + (a1 + (b1 * rs2 + c1 * rs2 ** 2 + d1 * rs2 ** & 768 3) * LOG(0.1D1 + 0.1D1 / (e1 * rs2 + f1 * rs2 ** (0.3D1 / 0.2D1) & 769 + g1 * rs2 ** 2 + h1 * rs2 ** 3))) * DBLE(zeta ** 2) + (a2 + (b2 & 770 * rs2 + c2 * rs2 ** 2 + d2 * rs2 ** 3) * LOG(0.1D1 + 0.1D1 / (e2 * & 771 rs2 + f2 * rs2 ** (0.3D1 / 0.2D1) + g2 * rs2 ** 2 + h2 * rs2 ** 3 & 772 ))) * DBLE(zeta ** 4) 773 ! 774 cg = -beta * EXP(-beta * rs2) * (-0.2D1 / 0.3D1 * SQRT(0.2D1) * & 775 DBLE((1 + zeta) ** (0.3D1 / 0.2D1) + (1 - zeta) ** (0.3D1 / 0.2D1)) & 776 / pi / rs2 + 0.4D1 / 0.3D1 * (0.1D1 + 0.3D1 / 0.8D1 * DBLE(zeta ** & 777 2) + 0.3D1 / 0.128D3 * DBLE(zeta ** 4)) * SQRT(0.2D1) / pi / rs2) & 778 + (EXP(-beta * rs2) - 0.1D1) * (0.2D1 / 0.3D1 * SQRT(0.2D1) * DBLE & 779 ((1 + zeta) ** (0.3D1 / 0.2D1) + (1 - zeta) ** (0.3D1 / 0.2D1)) / & 780 pi / rs2 ** 2 - 0.4D1 / 0.3D1 * (0.1D1 + 0.3D1 / 0.8D1 * DBLE(zeta & 781 ** 2) + 0.3D1 / 0.128D3 * DBLE(zeta ** 4)) * SQRT(0.2D1) / pi / & 782 rs2 ** 2) + (b0 + 0.2D1 * c0 * rs2 + 0.3D1 * d0 * rs2 ** 2) * LOG( & 783 0.1D1 + 0.1D1 / (e0 * rs2 + f0 * rs2 ** (0.3D1 / 0.2D1) + g0 * rs2 & 784 ** 2 + h0 * rs2 ** 3)) - (b0 * rs2 + c0 * rs2 ** 2 + d0 * rs2 ** & 785 3) / (e0 * rs2 + f0 * rs2 ** (0.3D1 / 0.2D1) + g0 * rs2 ** 2 + h0 & 786 * rs2 ** 3) ** 2 * (e0 + 0.3D1 / 0.2D1 * f0 * SQRT(rs2) + 0.2D1 * & 787 g0 * rs2 + 0.3D1 * h0 * rs2 ** 2) / (0.1D1 + 0.1D1 / (e0 * rs2 + f0 & 788 * rs2 ** (0.3D1 / 0.2D1) + g0 * rs2 ** 2 + h0 * rs2 ** 3)) + (( & 789 b1 + 0.2D1 * c1 * rs2 + 0.3D1 * d1 * rs2 ** 2) * LOG(0.1D1 + 0.1D1 & 790 / (e1 * rs2 + f1 * rs2 ** (0.3D1 / 0.2D1) + g1 * rs2 ** 2 + h1 * & 791 rs2 ** 3)) - (b1 * rs2 + c1 * rs2 ** 2 + d1 * rs2 ** 3) / (e1 * rs2 & 792 + f1 * rs2 ** (0.3D1 / 0.2D1) + g1 * rs2 ** 2 + h1 * rs2 ** 3) ** & 793 2 * (e1 + 0.3D1 / 0.2D1 * f1 * SQRT(rs2) + 0.2D1 * g1 * rs2 + & 794 0.3D1 * h1 * rs2 ** 2) / (0.1D1 + 0.1D1 / (e1 * rs2 + f1 * rs2 ** ( & 795 0.3D1 / 0.2D1) + g1 * rs2 ** 2 + h1 * rs2 ** 3))) * DBLE(zeta ** 2) & 796 + ((b2 + 0.2D1 * c2 * rs2 + 0.3D1 * d2 * rs2 ** 2) * LOG(0.1D1 + & 797 0.1D1 / (e2 * rs2 + f2 * rs2 ** (0.3D1 / 0.2D1) + g2 * rs2 ** 2 + h2 & 798 * rs2 ** 3)) - (b2 * rs2 + c2 * rs2 ** 2 + d2 * rs2 ** 3) / (e2 & 799 * rs2 + f2 * rs2 ** (0.3D1 / 0.2D1) + g2 * rs2 ** 2 + h2 * rs2 ** & 800 3) ** 2 * (e2 + 0.3D1 / 0.2D1 * f2 * SQRT(rs2) + 0.2D1 * g2 * rs2 & 801 + 0.3D1 * h2 * rs2 ** 2) / (0.1D1 + 0.1D1 / (e2 * rs2 + f2 * rs2 ** & 802 (0.3D1 / 0.2D1) + g2 * rs2 ** 2 + h2 * rs2 ** 3))) * DBLE(zeta ** & 803 4) 804 ! 805 deps2ddrs2=cg 806 ! 807 ! GGA-2D 808 ! 809 depsGGAdrs = ff*(-depscPBEdrs+deps2ddrs2*drs2ddrs3d) 810 depsGGAdt = dffdt*(-epscPBE+eps2d)+ff* & 811 (-depscPBEdt+deps2ddrs2*drs2ddt) 812 ! 813 ec = rho*(ff*(-epscPBE+eps2d)) 814 ! 815 decdn = ff*(-epscPBE+eps2d)+rho*depsGGAdrs*drsdn+ & 816 rho*depsGGAdt*dtdn 817 ! 818 decdgr = rho*depsGGAdt*dtdgr 819 ! 820 sc = ec 821 v1c = decdn 822 v2c = decdgr/gr 823 ! 824 RETURN 825 ! 826END SUBROUTINE cpbe2d 827! ! 828END MODULE 829 830