1! 2! Copyright (C) 2001-2014 Quantum ESPRESSO group 3! This file is distributed under the terms of the 4! GNU General Public License. See the file `License' 5! in the root directory of the present distribution, 6! or http://www.gnu.org/copyleft/gpl.txt . 7! 8 9MODULE corr_lda !<GPU:corr_lda=>corr_lda_gpu> 10 11CONTAINS 12 13!------------------------------------------------------------------------- 14SUBROUTINE pz( rs, iflag, ec, vc ) !<GPU:DEVICE> 15 !----------------------------------------------------------------------- 16 !! LDA parametrization from Monte Carlo DATA: 17 ! 18 !! * iflag=1: J.P. Perdew and A. Zunger, PRB 23, 5048 (1981); 19 !! * iflag=2: G. Ortiz and P. Ballone, PRB 50, 1391 (1994). 20 ! 21 USE kinds, ONLY: DP 22 ! 23 IMPLICIT NONE 24 ! 25 INTEGER, INTENT(IN) :: iflag !<GPU:VALUE> 26 !! see routine comments 27 REAL(DP), INTENT(IN) :: rs 28 !! Wigner-Seitz radius 29 REAL(DP), INTENT(OUT) :: ec 30 !! correlation energy 31 REAL(DP), INTENT(OUT) :: vc 32 !! correlation potential 33 ! 34 ! ... local variables 35 ! 36 REAL(DP) :: a(2), b(2), c(2), d(2), gc(2), b1(2), b2(2) 37 REAL(DP) :: lnrs, rs12, ox, dox 38 ! 39 DATA a / 0.0311d0, 0.031091d0 /, b / -0.048d0, -0.046644d0 /, & 40 c / 0.0020d0, 0.00419d0 /, d / -0.0116d0, -0.00983d0 / 41 ! 42 DATA gc / -0.1423d0, -0.103756d0 /, b1 / 1.0529d0, 0.56371d0 /, & 43 b2 / 0.3334d0, 0.27358d0 / 44 ! 45 IF ( rs < 1.0d0 ) THEN 46 ! 47 ! high density formula 48 lnrs = LOG(rs) 49 ! 50 ec = a(iflag)*lnrs + b(iflag) + c(iflag)*rs*lnrs + d(iflag)*rs 51 ! 52 vc = a(iflag)*lnrs + ( b(iflag) - a(iflag)/3.d0 ) + 2.d0/3.d0 * & 53 c(iflag)*rs*lnrs + ( 2.d0*d(iflag) - c(iflag) )/3.d0*rs 54 ELSE 55 ! 56 ! interpolation formula 57 rs12 = SQRT(rs) 58 ! 59 ox = 1.d0 + b1(iflag)*rs12 + b2(iflag)*rs 60 dox = 1.d0 + 7.d0/6.d0*b1(iflag)*rs12 + 4.d0/3.d0*b2(iflag)*rs 61 ! 62 ec = gc(iflag)/ox 63 vc = ec*dox/ox 64 ! 65 ENDIF 66 ! 67 RETURN 68 ! 69END SUBROUTINE pz 70! 71! 72!----------------------------------------------------------------------- 73SUBROUTINE pzKZK( rs, ec, vc, vol ) !<GPU:DEVICE> 74 !----------------------------------------------------------------------- 75 !! LDA parametrization from Monte Carlo DATA: 76 ! 77 !! * iflag=1: J.P. Perdew and A. Zunger, PRB 23, 5048 (1981) 78 !! * iflag=2: G. Ortiz and P. Ballone, PRB 50, 1391 (1994) 79 ! 80 USE kinds, ONLY: DP 81 ! 82 IMPLICIT NONE 83 ! 84 REAL(DP), INTENT(IN) :: rs 85 !! Wigner-Seitz radius 86 REAL(DP), INTENT(OUT) :: ec 87 !! correlation energy 88 REAL(DP), INTENT(OUT) :: vc 89 !! correlation potential 90 REAL(DP) :: vol !<GPU:VALUE> 91 !! volume element 92 ! 93 ! ... local variables 94 ! 95 INTEGER :: iflag, kr 96 REAL(DP) :: ec0(2), vc0(2), ec0p 97 REAL(DP) :: a(2), b(2), c(2), d(2), gc(2), b1(2), b2(2) 98 REAL(DP) :: lnrs, rs12, ox, dox, lnrsk, rsk 99 REAL(DP) :: a1, grs, g1, g2, g3, g4, dL, gh, gl, grsp 100 REAL(DP) :: f3, f2, f1, f0, pi 101 REAL(DP) :: D1, D2, D3, P1, P2, ry2h 102 ! 103 DATA a / 0.0311_DP, 0.031091_DP /, b / -0.048_DP, -0.046644_DP /, & 104 c / 0.0020_DP, 0.00419_DP /, d / -0.0116_DP, -0.00983_DP / 105 ! 106 DATA gc / -0.1423_DP, -0.103756_DP /, b1 / 1.0529_DP, 0.56371_DP /, & 107 b2 / 0.3334_DP, 0.27358_DP / 108 ! 109 DATA a1 / -2.2037_DP /, g1 / 0.1182_DP/, g2 / 1.1656_DP/, g3 / -5.2884_DP/, & 110 g4 / -1.1233_DP / 111 ! 112 DATA ry2h / 0.5_DP / 113 ! 114 iflag = 1 115 pi = 4.d0 * ATAN(1.d0) 116 dL = vol**(1.d0/3.d0) 117 gh = 0.5d0 * dL / (2.d0 * pi)**(1.d0/3.d0) 118 gl = dL * (3.d0 / 2.d0 / pi)**(1.d0/3.d0) 119 ! 120 rsk = gh 121 ! 122 DO kr = 1, 2 123 ! 124 lnrsk = LOG(rsk) 125 IF (rsk < 1.0d0) THEN 126 ! high density formula 127 ec0(kr) = a(iflag) *lnrsk + b(iflag) + c(iflag) * rsk * lnrsk + d(iflag) * rsk 128 vc0(kr) = a(iflag) * lnrsk + (b(iflag) - a(iflag) / 3.d0) + 2.d0 / & 129 3.d0 * c (iflag) * rsk * lnrsk + (2.d0 * d (iflag) - c (iflag) ) & 130 / 3.d0 * rsk 131 ! 132 ELSE 133 ! interpolation formula 134 rs12 = SQRT(rsk) 135 ox = 1.d0 + b1 (iflag) * rs12 + b2 (iflag) * rsk 136 dox = 1.d0 + 7.d0 / 6.d0 * b1 (iflag) * rs12 + 4.d0 / 3.d0 * & 137 b2 (iflag) * rsk 138 ec0(kr) = gc (iflag) / ox 139 vc0(kr) = ec0(kr) * dox / ox 140 ! 141 ENDIF 142 ! 143 grs = g1 * rsk * lnrsk + g2 * rsk + g3 * rsk**1.5d0 + g4 * rsk**2.d0 144 grsp = g1 * lnrsk + g1 + g2 + 1.5d0 * g3 * rsk**0.5d0 + 2.d0 * g4 * rsk 145 ! 146 ec0(kr) = ec0(kr) + (-a1 * rsk / dL**2.d0 + grs / dL**3.d0) * ry2h 147 vc0(kr) = vc0(kr) + (-2.d0 * a1 * rsk / dL**2.d0 / 3.d0 + & 148 grs / dL**3.d0 - grsp * rsk / 3.d0 / dL**3.d0) * ry2h 149 ! 150 rsk = rs 151 ! 152 ENDDO 153 ! 154 lnrs = LOG(rs) 155 ! 156 IF (rs <= gh) THEN 157 ec = ec0(2) 158 vc = vc0(2) 159 ELSE 160 IF ( rs <= gl ) THEN 161 ec0p = 3.d0 * (ec0(1) - vc0(1)) / gh 162 P1 = 3.d0 * ec0(1) - gh * ec0p 163 P2 = ec0p 164 D1 = gl - gh 165 D2 = gl**2.d0 - gh**2.d0 166 D3 = gl**3.d0 - gh**3.d0 167 f2 = 2.d0 * gl**2.d0 * P2 * D1 + D2 * P1 168 f2 = f2 / (-(2.d0*gl*D1)**2.d0 + 4.d0*gl*D1*D2 - D2**2.d0 ) 169 f3 = - (P2 + 2.d0*D1*f2) / (3.d0 * D2) 170 f1 = - (P1 + D2 * f2) / (2.d0 * D1) 171 f0 = - gl * (gl * f2 + 2.d0 * f1) / 3.d0 172 ! 173 ec = f3 * rs**3.d0 + f2 * rs**2.d0 + f1 * rs + f0 174 vc = f2 * rs**2.d0 / 3.d0 + f1 * 2.d0 * rs / 3.d0 + f0 175 ELSE 176 ec = 0.d0 177 vc = 0.d0 178 ENDIF 179 ENDIF 180 ! 181 ! 182 RETURN 183 ! 184END SUBROUTINE pzKZK 185! 186! 187!----------------------------------------------------------------------- 188SUBROUTINE vwn( rs, ec, vc ) !<GPU:DEVICE> 189 !----------------------------------------------------------------------- 190 !! S.H. Vosko, L. Wilk, and M. Nusair, Can. J. Phys. 58, 1200 (1980). 191 ! 192 USE kinds, ONLY: DP 193 ! 194 IMPLICIT NONE 195 ! 196 REAL(DP), INTENT(IN) :: rs 197 !! Wigner-Seitz radius 198 REAL(DP), INTENT(OUT) :: ec 199 !! correlation energy 200 REAL(DP), INTENT(OUT) :: vc 201 !! correlation potential 202 ! 203 ! ... local variables 204 ! 205 REAL(DP), PARAMETER :: a = 0.0310907d0, b = 3.72744d0, & 206 c = 12.9352d0, x0 = -0.10498d0 207 REAL(DP) :: q, f1, f2, f3, rs12, fx, qx, tx, tt 208 ! 209 ! 210 q = SQRT( 4.d0*c - b*b ) 211 f1 = 2.d0*b/q 212 f2 = b*x0 / ( x0*x0 + b*x0 + c ) 213 f3 = 2.d0*( 2.d0 * x0 + b ) / q 214 ! 215 rs12 = SQRT(rs) 216 fx = rs + b*rs12 + c 217 qx = ATAN( q / (2.d0*rs12 + b) ) 218 ! 219 ec = a * ( LOG( rs/fx ) + f1*qx - f2*( LOG( (rs12 - x0)**2 / fx) & 220 + f3*qx) ) 221 ! 222 tx = 2.d0*rs12 + b 223 tt = tx*tx + q*q 224 ! 225 vc = ec - rs12*a / 6.d0*( 2.d0 / rs12 - tx/fx - 4.d0*b/tt - f2 * & 226 (2.d0 / (rs12 - x0) - tx/fx - 4.d0*(2.d0 * x0 + b)/tt) ) 227 ! 228 RETURN 229 ! 230END SUBROUTINE vwn 231! 232! 233!----------------------------------------------------------------------- 234SUBROUTINE vwn1_rpa( rs, ec, vc ) !<GPU:DEVICE> 235 !----------------------------------------------------------------------- 236 !! S.H. Vosko, L. Wilk, and M. Nusair, Can. J. Phys. 58, 1200 (1980). 237 ! 238 USE kinds, ONLY: DP 239 ! 240 IMPLICIT NONE 241 ! 242 REAL(DP), INTENT(IN) :: rs 243 !! Wigner-Seitz radius 244 REAL(DP), INTENT(OUT) :: ec 245 !! correlation energy 246 REAL(DP), INTENT(OUT) :: vc 247 !! correlation potential 248 ! 249 ! ... local variables 250 ! 251 REAL(DP), PARAMETER :: a = 0.0310907_DP, b = 13.0720_DP, & 252 c = 42.7198_DP, x0 = -0.409286_DP 253 REAL(DP) :: q, f1, f2, f3, rs12, fx, qx, tx, tt 254 ! 255 q = SQRT(4.d0*c - b*b) 256 f1 = 2.d0*b/q 257 f2 = b*x0 / (x0*x0 + b*x0 + c) 258 f3 = 2.d0 * (2.d0 * x0 + b) / q 259 ! 260 rs12 = SQRT(rs) 261 fx = rs + b * rs12 + c 262 qx = ATAN(q / (2.d0 * rs12 + b) ) 263 ! 264 ec = a*( LOG(rs/fx) + f1*qx - f2*(LOG((rs12 - x0)**2/fx) + f3*qx) ) 265 ! 266 tx = 2.d0*rs12 + b 267 tt = tx*tx + q*q 268 ! 269 vc = ec - rs12*a/6.d0*( 2.d0/rs12 - tx/fx - 4.d0*b/tt - f2 * & 270 (2.d0/(rs12 - x0) - tx/fx - 4.d0*(2.d0*x0 + b)/tt) ) 271 ! 272 ! 273 RETURN 274 ! 275END SUBROUTINE vwn1_rpa 276! 277! 278!----------------------------------------------------------------------- 279SUBROUTINE lyp( rs, ec, vc ) !<GPU:DEVICE> 280 !----------------------------------------------------------------------- 281 !! C. Lee, W. Yang, and R.G. Parr, PRB 37, 785 (1988). 282 !! LDA part only. 283 ! 284 USE kinds, ONLY: DP 285 ! 286 IMPLICIT NONE 287 ! 288 REAL(DP), INTENT(IN) :: rs 289 !! Wigner-Seitz radius 290 REAL(DP), INTENT(OUT) :: ec 291 !! correlation energy 292 REAL(DP), INTENT(OUT) :: vc 293 !! correlation potential 294 ! 295 ! ... local variables 296 ! 297 REAL(DP), PARAMETER :: a=0.04918d0, b=0.132d0*2.87123400018819108d0 298 REAL(DP), PARAMETER :: pi43=1.61199195401647d0 299 ! pi43=(4pi/3)^(1/3) 300 REAL(DP), PARAMETER :: c=0.2533d0*pi43, d=0.349d0*pi43 301 REAL(DP) :: ecrs, ox 302 ! 303 ! 304 ecrs = b*EXP( -c*rs ) 305 ox = 1.d0 / (1.d0 + d*rs) 306 ! 307 ec = - a*ox*(1.d0 + ecrs) 308 vc = ec - rs/3.d0*a*ox*( d*ox + ecrs*(d*ox + c) ) 309 ! 310 RETURN 311 ! 312END SUBROUTINE lyp 313! 314! 315!----------------------------------------------------------------------- 316SUBROUTINE pw( rs, iflag, ec, vc ) !<GPU:DEVICE> 317 !----------------------------------------------------------------------- 318 !! * iflag=1: J.P. Perdew and Y. Wang, PRB 45, 13244 (1992) 319 !! * iflag=2: G. Ortiz and P. Ballone, PRB 50, 1391 (1994) 320 ! 321 USE kinds, ONLY: DP 322 ! 323 IMPLICIT NONE 324 ! 325 REAL(DP), INTENT(IN) :: rs 326 !! Wigner-Seitz radius 327 INTEGER, INTENT(IN) :: iflag !<GPU:VALUE> 328 !! see routine comments 329 REAL(DP), INTENT(OUT) :: ec 330 !! correlation energy 331 REAL(DP), INTENT(OUT) :: vc 332 !! correlation potential 333 ! 334 ! ... local variables 335 ! 336 REAL(DP), PARAMETER :: a=0.031091d0, b1=7.5957d0, b2=3.5876d0, c0=a, & 337 c1=0.046644d0, c2=0.00664d0, c3=0.01043d0, d0=0.4335d0, & 338 d1=1.4408d0 339 REAL(DP) :: lnrs, rs12, rs32, rs2, om, dom, olog 340 REAL(DP) :: a1(2), b3(2), b4(2) 341 DATA a1 / 0.21370d0, 0.026481d0 /, b3 / 1.6382d0, -0.46647d0 /, & 342 b4 / 0.49294d0, 0.13354d0 / 343 ! 344 ! high- and low-density formulae implemented but not used in PW case 345 ! (reason: inconsistencies in PBE/PW91 functionals). 346 ! 347 IF ( rs < 1d0 .AND. iflag == 2 ) THEN 348 ! 349 ! high density formula 350 lnrs = LOG(rs) 351 ec = c0 * lnrs - c1 + c2 * rs * lnrs - c3 * rs 352 vc = c0 * lnrs - (c1 + c0 / 3.d0) + 2.d0 / 3.d0 * c2 * rs * & 353 lnrs - (2.d0 * c3 + c2) / 3.d0 * rs 354 ! 355 ELSEIF ( rs > 100.d0 .AND. iflag == 2 ) THEN 356 ! 357 ! low density formula 358 ec = - d0 / rs + d1 / rs**1.5d0 359 vc = - 4.d0 / 3.d0 * d0 / rs + 1.5d0 * d1 / rs**1.5d0 360 ! 361 ELSE 362 ! 363 ! interpolation formula 364 rs12 = SQRT(rs) 365 rs32 = rs * rs12 366 rs2 = rs**2 367 om = 2.d0*a*( b1*rs12 + b2*rs + b3(iflag) * rs32 + b4(iflag)*rs2 ) 368 dom = 2.d0*a*( 0.5d0 * b1 * rs12 + b2 * rs + 1.5d0 * b3(iflag) * & 369 rs32 + 2.d0 * b4(iflag) * rs2 ) 370 olog = LOG( 1.d0 + 1.0d0 / om ) 371 ! 372 ec = - 2.d0 * a * (1.d0 + a1(iflag) * rs) * olog 373 vc = - 2.d0 * a * (1.d0 + 2.d0 / 3.d0 * a1(iflag) * rs) & 374 * olog - 2.d0 / 3.d0 * a * (1.d0 + a1(iflag) * rs) * dom / & 375 (om * (om + 1.d0) ) 376 ! 377 ENDIF 378 ! 379 ! 380 RETURN 381 ! 382END SUBROUTINE pw 383! 384! 385!----------------------------------------------------------------------- 386SUBROUTINE wignerc( rs, ec, vc ) !<GPU:DEVICE> 387 !----------------------------------------------------------------------- 388 !! Wigner correlation. 389 ! 390 USE kinds, ONLY: DP 391 ! 392 IMPLICIT NONE 393 ! 394 REAL(DP), INTENT(IN) :: rs 395 !! Wigner-Seitz radius 396 REAL(DP), INTENT(OUT) :: ec 397 !! correlation energy 398 REAL(DP), INTENT(OUT) :: vc 399 !! correlation potential 400 ! 401 ! ... local variables 402 ! 403 REAL(DP) :: rho13 !rho13=rho^(1/3) 404 REAL(DP), PARAMETER :: pi34 = 0.6203504908994d0 405 ! pi34 = (3/4pi)^(1/3) 406 ! 407 ! 408 rho13 = pi34 / rs 409 vc = - rho13 * ( (0.943656d0 + 8.8963d0 * rho13) / (1.d0 + & 410 12.57d0 * rho13)**2 ) 411 ec = - 0.738d0 * rho13 * ( 0.959d0 / (1.d0 + 12.57d0 * rho13) ) 412 ! 413 RETURN 414 ! 415END SUBROUTINE wignerc 416! 417! 418!----------------------------------------------------------------------- 419SUBROUTINE hl( rs, ec, vc ) !<GPU:DEVICE> 420 !----------------------------------------------------------------------- 421 !! L. Hedin and B.I. Lundqvist, J. Phys. C 4, 2064 (1971). 422 ! 423 USE kinds, ONLY: DP 424 ! 425 IMPLICIT NONE 426 ! 427 REAL(DP), INTENT(IN) :: rs 428 !! Wigner-Seitz radius 429 REAL(DP), INTENT(OUT) :: ec 430 !! correlation energy 431 REAL(DP), INTENT(OUT) :: vc 432 !! correlation potential 433 ! 434 ! ... local variables 435 ! 436 REAL(DP) :: a, x 437 ! 438 ! 439 a = LOG(1.0d0 + 21.d0/rs) 440 x = rs / 21.0d0 441 ec = a + (x**3 * a - x*x) + x/2.d0 - 1.0d0/3.0d0 442 ec = - 0.0225d0 * ec 443 vc = - 0.0225d0 * a 444 ! 445 RETURN 446 ! 447END SUBROUTINE hl 448! 449! 450!----------------------------------------------------------------------- 451SUBROUTINE gl( rs, ec, vc ) !<GPU:DEVICE> 452 !--------------------------------------------------------------------- 453 !! O. Gunnarsson and B. I. Lundqvist, PRB 13, 4274 (1976). 454 ! 455 USE kinds, ONLY: DP 456 ! 457 IMPLICIT NONE 458 ! 459 REAL(DP), INTENT(IN) :: rs 460 !! Wigner-Seitz radius 461 REAL(DP), INTENT(OUT) :: ec 462 !! correlation energy 463 REAL(DP), INTENT(OUT) :: vc 464 !! correlation potential 465 ! 466 ! ... local variables 467 ! 468 REAL(DP) :: x 469 REAL(DP), PARAMETER :: c=0.0333d0, r=11.4d0 470 ! c=0.0203, r=15.9 for the paramagnetic case 471 ! 472 x = rs / r 473 vc = - c * LOG(1.d0 + 1.d0 / x) 474 ec = - c * ( (1.d0 + x**3) * LOG(1.d0 + 1.d0 / x) - 1.0d0 / & 475 3.0d0 + x * (0.5d0 - x) ) 476 ! 477 RETURN 478 ! 479END SUBROUTINE gl 480! 481! 482! ... LSDA 483! 484!----------------------------------------------------------------------- 485SUBROUTINE pz_polarized( rs, ec, vc ) !<GPU:DEVICE> 486 !----------------------------------------------------------------------- 487 !! J.P. Perdew and A. Zunger, PRB 23, 5048 (1981). 488 !! spin-polarized energy and potential. 489 ! 490 USE kinds, ONLY : DP 491 ! 492 IMPLICIT NONE 493 ! 494 REAL(DP), INTENT(IN) :: rs 495 !! Wigner-Seitz radius 496 REAL(DP), INTENT(OUT) :: ec 497 !! correlation energy 498 REAL(DP), INTENT(OUT) :: vc 499 !! correlation potential 500 ! 501 ! ... local variables 502 ! 503 REAL(DP), PARAMETER :: a=0.01555d0, b=-0.0269d0, c=0.0007d0, & 504 d=-0.0048d0, gc=-0.0843d0, b1=1.3981d0, b2=0.2611d0 505 REAL(DP) :: lnrs, rs12, ox, dox 506 REAL(DP), PARAMETER :: xcprefact=0.022575584d0, pi34=0.6203504908994d0 507 ! REAL(DP) :: betha, etha, csi, prefact 508 ! 509 ! 510 IF ( rs < 1.0d0 ) THEN 511 ! high density formula 512 lnrs = LOG(rs) 513 ec = a * lnrs + b + c * rs * lnrs + d * rs 514 vc = a * lnrs + (b - a / 3.d0) + 2.d0 / 3.d0 * c * rs * lnrs + & 515 (2.d0 * d-c) / 3.d0 * rs 516 ELSE 517 ! interpolation formula 518 rs12 = SQRT(rs) 519 ox = 1.d0 + b1 * rs12 + b2 * rs 520 dox = 1.d0 + 7.d0 / 6.d0 * b1 * rs12 + 4.d0 / 3.d0 * b2 * rs 521 ec = gc / ox 522 vc = ec * dox / ox 523 ENDIF 524 ! 525 ! IF ( lxc_rel ) THEN 526 ! betha = prefact * pi34 / rs 527 ! etha = DSQRT( 1 + betha**2 ) 528 ! csi = betha + etha 529 ! prefact = 1.0D0 - (3.0D0/2.0D0) * ( (betha*etha - LOG(csi))/betha**2 )**2 530 ! ec = ec * prefact 531 ! vc = vc * prefact 532 ! ENDIF 533 ! 534 RETURN 535 ! 536END SUBROUTINE pz_polarized 537! 538! 539!----------------------------------------------------------------------- 540SUBROUTINE pz_spin( rs, zeta, ec, vc_up, vc_dw ) !<GPU:DEVICE> 541 !----------------------------------------------------------------------- 542 !! Perdew and Zunger, PRB 23, 5048 (1981). Spin polarized case. 543 ! 544 USE kinds, ONLY : DP 545 ! 546 IMPLICIT NONE 547 ! 548 REAL(DP), INTENT(IN) :: rs 549 !! Wigner-Seitz radius 550 REAL(DP), INTENT(IN) :: zeta 551 !! zeta = (rho_up - rho_dw) / rho_tot 552 REAL(DP), INTENT(OUT) :: ec 553 !! correlation energy 554 REAL(DP), INTENT(OUT) :: vc_up, vc_dw 555 !! correlation potential (up, down) 556 ! 557 ! ... local variables 558 ! 559 REAL(DP) :: ecu, vcu, ecp, vcp 560 REAL(DP) :: fz, dfz 561 REAL(DP), PARAMETER :: p43=4.0d0/3.d0, third=1.d0/3.d0 562 ! 563 ! unpolarized part (Perdew-Zunger formula) 564 CALL pz( rs, 1, ecu, vcu ) !<GPU:pz=>pz_d> 565 ! 566 ! polarization contribution 567 CALL pz_polarized( rs, ecp, vcp ) !<GPU:pz_polarized=>pz_polarized_d> 568 ! 569 fz = ( (1.0d0 + zeta)**p43 + (1.d0 - zeta)**p43 - 2.d0) / & 570 (2.d0**p43 - 2.d0) 571 dfz = p43 * ( (1.0d0 + zeta)**third- (1.d0 - zeta)**third) & 572 / (2.d0**p43 - 2.d0) 573 ! 574 ec = ecu + fz * (ecp - ecu) 575 vc_up = vcu + fz * (vcp - vcu) + (ecp - ecu) * dfz * ( 1.d0 - zeta) 576 vc_dw = vcu + fz * (vcp - vcu) + (ecp - ecu) * dfz * (-1.d0 - zeta) 577 ! 578 RETURN 579 ! 580END SUBROUTINE pz_spin 581! 582!------------------------------------------------------------------------------- 583SUBROUTINE vwn_spin( rs, zeta, ec, vc_up, vc_dw ) !<GPU:DEVICE> 584 !------------------------------------------------------------------------------ 585 !! S.H. Vosko, L. Wilk, and M. Nusair. Spin polarized case. 586 ! 587 USE kinds, ONLY: DP 588 ! 589 IMPLICIT NONE 590 ! 591 REAL(DP), INTENT(IN) :: rs 592 !! Wigner-Seitz radius 593 REAL(DP), INTENT(IN) :: zeta 594 !! zeta = (rho_up - rho_dw) / rho_tot 595 REAL(DP), INTENT(OUT) :: ec 596 !! correlation energy 597 REAL(DP), INTENT(OUT) :: vc_up, vc_dw 598 !! correlation potential (up, down) 599 ! 600 ! ... local variables 601 ! 602 REAL(DP) :: zeta3, zeta4, trup, trdw, trup13, trdw13, fz, dfz, fzz4 603 REAL(DP) :: SQRTrs, ecP, ecF, ac, De, vcP, vcF, dac, dec1, dec2 604 REAL(DP) :: cfz, cfz1, cfz2, iddfz0 605 ! 606 ! parameters: e_c/para, e_c/ferro, alpha_c 607 REAL(DP), PARAMETER :: & 608 A(3) = (/ 0.0310907_DP, 0.01554535_DP, -0.01688686394039_DP /), & 609 x0(3) = (/ -0.10498_DP, -0.32500_DP, -0.0047584_DP /), & 610 b(3) = (/3.72744_DP, 7.06042_DP, 1.13107_DP /), & 611 c(3) = (/ 12.9352_DP, 18.0578_DP, 13.0045_DP /),& 612 Q(3) = (/ 6.15199081975908_DP, 4.73092690956011_DP, 7.12310891781812_DP /), & 613 tbQ(3) = (/ 1.21178334272806_DP, 2.98479352354082_DP, 0.31757762321188_DP /), & 614 fx0(3) = (/ 12.5549141492_DP, 15.8687885_DP, 12.99914055888256_DP /), & 615 bx0fx0(3) = (/ -0.03116760867894_DP, -0.14460061018521_DP, -0.00041403379428_DP /) 616 ! 617 ! N.B.: A is expressed in Hartree 618 ! Q = SQRT(4*c - b^2) 619 ! tbQ = 2*b/Q 620 ! fx0 = X(x_0) = x_0^2 + b*x_0 + c 621 ! bx0fx0 = b*x_0/X(x_0) 622 ! 623 ! 624 ! coefficients for f(z), df/dz, ddf/ddz(0) 625 cfz = 2.0_DP**(4.0_DP/3.0_DP) - 2.0_DP 626 cfz1 = 1.0_DP / cfz 627 cfz2 = 4.0_DP/3.0_DP * cfz1 628 iddfz0 = 9.0_DP / 8.0_DP *cfz 629 ! 630 SQRTrs = SQRT(rs) 631 zeta3 = zeta**3 632 zeta4 = zeta3*zeta 633 trup = 1.0_DP + zeta 634 trdw = 1.0_DP - zeta 635 trup13 = trup**(1.0_DP/3.0_DP) 636 trdw13 = trdw**(1.0_DP/3.0_DP) 637 fz = cfz1 * (trup13*trup + trdw13*trdw - 2.0_DP) ! f(zeta) 638 dfz = cfz2 * (trup13 - trdw13) ! df / dzeta 639 ! 640 CALL padefit_ParSet1( sqrtrs, 1, ecP, vcP ) ! ecF = e_c Paramagnetic !<GPU:padefit_ParSet1=>padefit_ParSet1_d> 641 CALL padefit_ParSet1( sqrtrs, 2, ecF, vcF ) ! ecP = e_c Ferromagnetic !<GPU:padefit_ParSet1=>padefit_ParSet1_d> 642 CALL padefit_ParSet1( sqrtrs, 3, ac, dac ) ! ac = "spin stiffness" !<GPU:padefit_ParSet1=>padefit_ParSet1_d> 643 ! 644 ac = ac * iddfz0 645 dac = dac * iddfz0 646 De = ecF - ecP - ac ! e_c[F]-e_c[P]-alpha_c/(ddf/ddz(z=0)) 647 fzz4 = fz * zeta4 648 ec = ecP + ac * fz + De * fzz4 649 ! 650 dec1 = vcP + dac*fz + (vcF - vcP - dac) * fzz4 ! e_c - (r_s/3)*(de_c/dr_s) 651 dec2 = ac*dfz + De*(4.0_DP*zeta3*fz + zeta4*dfz) ! de_c/dzeta 652 ! 653 vc_up = dec1 + (1.0_DP - zeta)*dec2 ! v_c[s] = e_c - (r_s/3)*(de_c/dr_s) 654 vc_dw = dec1 - (1.0_DP + zeta)*dec2 ! + [sign(s)-zeta]*(de_c/dzeta) 655 ! 656END SUBROUTINE 657! 658! 659!---- 660SUBROUTINE padefit_ParSet1( x, i, fit, dfit ) !<GPU:DEVICE> 661 !---- 662 !! It implements formula [4.4] in: 663 !! S.H. Vosko, L. Wilk, and M. Nusair, Can. J. Phys. 58, 1200 (1980) 664 ! 665 USE kinds, ONLY: DP 666 ! 667 IMPLICIT NONE 668 ! 669 REAL(DP) :: x 670 !! input: x is SQRT(r_s) 671 INTEGER :: i 672 !! input: i is the index of the fit 673 REAL(DP) :: fit 674 !! output: Pade fit calculated in x [eq. 4.4] 675 REAL(DP) :: dfit 676 !! output: dfit/drho = fit - (rs/3)*dfit/drs = ec - (x/6)*dfit/dx 677 ! 678 ! ... local variables 679 ! 680 REAL(DP), PARAMETER :: & 681 A(3) = (/ 0.0310907d0, 0.01554535d0, -0.01688686394039d0 /), & 682 x0(3) = (/ -0.10498d0, -0.32500d0, -0.0047584d0 /), & 683 b(3) = (/3.72744d0, 7.06042d0, 1.13107d0 /), & 684 c(3) = (/ 12.9352d0, 18.0578d0, 13.0045d0 /),& 685 Q(3) = (/ 6.15199081975908d0, 4.73092690956011d0, 7.12310891781812d0 /), & 686 tbQ(3) = (/ 1.21178334272806d0, 2.98479352354082d0, 0.31757762321188d0 /), & 687 fx0(3) = (/ 12.5549141492d0, 15.8687885d0, 12.99914055888256d0 /), & 688 bx0fx0(3) = (/ -0.03116760867894d0, -0.14460061018521d0, -0.00041403379428d0 /) 689 690 REAL(DP) :: sqx, xx0, Qtxb, atg, fx 691 REAL(DP) :: txb, txbfx, itxbQ 692 ! 693 sqx = x * x ! x^2 = r_s 694 xx0 = x - x0(i) ! x - x_0 695 Qtxb = Q(i) / (2.0_DP*x + b(i)) ! Q / (2x+b) 696 atg = ATAN(Qtxb) ! tan^-1(Q/(2x+b)) 697 fx = sqx + b(i)*x + c(i) ! X(x) = x^2 + b*x + c 698 ! 699 fit = A(i) * ( LOG(sqx/fx) + tbQ(i)*atg - & 700 bx0fx0(i) * ( LOG(xx0*xx0/fx) + (tbQ(i) + 4.0_DP*x0(i)/Q(i)) * atg ) ) 701 ! 702 txb = 2.0_DP*x + b(i) 703 txbfx = txb / fx 704 itxbQ = 1.0_DP / (txb*txb + Q(i)*Q(i)) 705 ! 706 dfit = fit - A(i) / 3.0_DP + A(i)*x/6.0_DP * ( txbfx + 4.0_DP*b(i)*itxbQ + & 707 bx0fx0(i) * ( 2.0_DP/xx0 - txbfx - 4.0_DP*(b(i)+2.0_DP*x0(i))*itxbQ ) ) 708 ! 709END SUBROUTINE 710SUBROUTINE padefit_ParSet2( x, i, fit, dfit ) !<GPU:DEVICE> 711 !---- 712 !! It implements formula [4.4] in: 713 !! S.H. Vosko, L. Wilk, and M. Nusair, Can. J. Phys. 58, 1200 (1980) 714 ! 715 USE kinds, ONLY: DP 716 ! 717 IMPLICIT NONE 718 ! 719 REAL(DP) :: x 720 !! input: x is SQRT(r_s) 721 INTEGER :: i 722 !! input: i is the index of the fit 723 REAL(DP) :: fit 724 !! output: Pade fit calculated in x [eq. 4.4] 725 REAL(DP) :: dfit 726 !! output: dfit/drho = fit - (rs/3)*dfit/drs = ec - (x/6)*dfit/dx 727 ! 728 ! ... local variables 729 ! 730 REAL(DP), PARAMETER :: & 731 A(3) = (/ 0.0310907_DP, 0.01554535_DP, -0.01688686394039_DP /), & 732 x0(3) = (/ -0.409286_DP, -0.743294_DP, -0.228344_DP /), & 733 b(3) = (/ 13.0720_DP, 20.1231_DP, 1.06835_DP /), & 734 c(3) = (/ 42.7198_DP, 101.578_DP, 11.4813_DP /),& 735 Q(3) = (/ 0.044899888641577_DP, 1.171685277708971_DP, 6.692072046645942_DP /), & 736 tbQ(3) = (/ 582.273159042780890_DP, 34.348984975465861_DP, 0.319288254087299_DP /), & 737 fx0(3) = (/ 37.537128437796000_DP, 87.173106479036008_DP, 11.289489669936000_DP /), & 738 bx0fx0(3) = (/ -0.142530524167984_DP, -0.171582499414508_DP, -0.021608710360898_DP /) 739 740 REAL(DP) :: sqx, xx0, Qtxb, atg, fx 741 REAL(DP) :: txb, txbfx, itxbQ 742 ! 743 sqx = x * x ! x^2 = r_s 744 xx0 = x - x0(i) ! x - x_0 745 Qtxb = Q(i) / (2.0_DP*x + b(i)) ! Q / (2x+b) 746 atg = ATAN(Qtxb) ! tan^-1(Q/(2x+b)) 747 fx = sqx + b(i)*x + c(i) ! X(x) = x^2 + b*x + c 748 ! 749 fit = A(i) * ( LOG(sqx/fx) + tbQ(i)*atg - & 750 bx0fx0(i) * ( LOG(xx0*xx0/fx) + (tbQ(i) + 4.0_DP*x0(i)/Q(i)) * atg ) ) 751 ! 752 txb = 2.0_DP*x + b(i) 753 txbfx = txb / fx 754 itxbQ = 1.0_DP / (txb*txb + Q(i)*Q(i)) 755 ! 756 dfit = fit - A(i) / 3.0_DP + A(i)*x/6.0_DP * ( txbfx + 4.0_DP*b(i)*itxbQ + & 757 bx0fx0(i) * ( 2.0_DP/xx0 - txbfx - 4.0_DP*(b(i)+2.0_DP*x0(i))*itxbQ ) ) 758 ! 759END SUBROUTINE 760! 761! 762!----------------------------------------------------------------------------------------- 763SUBROUTINE vwn1_rpa_spin( rs, zeta, ec, vc_up, vc_dw ) !<GPU:DEVICE> 764 !--------------------------------------------------------------------------------------- 765 ! 766 USE kinds, ONLY: DP 767 ! 768 IMPLICIT NONE 769 ! 770 REAL(DP), INTENT(IN) :: rs 771 !! Wigner-Seitz radius 772 REAL(DP), INTENT(IN) :: zeta 773 !! zeta = (rho_up - rho_dw)/rho_tot 774 REAL(DP), INTENT(OUT) :: ec 775 !! correlation energy 776 REAL(DP), INTENT(OUT) :: vc_up, vc_dw 777 !! correlation potential (up, down) 778 ! 779 ! ... local variables 780 ! 781 REAL(DP) :: zeta3, zeta4, trup, trdw, trup13, trdw13, fz, dfz, fzz4 782 REAL(DP) :: SQRTrs, ecP, ecF, ac, De, vcP, vcF, dac, dec1, dec2 783 REAL(DP) :: cfz, cfz1, cfz2, iddfz0 784 ! 785 ! PARAMETERs: e_c/para, e_c/ferro, alpha_c 786 REAL(DP), PARAMETER :: & 787 A(3) = (/ 0.0310907_DP, 0.01554535_DP, -0.01688686394039_DP /), & 788 x0(3) = (/ -0.409286_DP, -0.743294_DP, -0.228344_DP /), & 789 b(3) = (/ 13.0720_DP, 20.1231_DP, 1.06835_DP /), & 790 c(3) = (/ 42.7198_DP, 101.578_DP, 11.4813_DP /),& 791 Q(3) = (/ 0.044899888641577_DP, 1.171685277708971_DP, 6.692072046645942_DP /), & 792 tbQ(3) = (/ 582.273159042780890_DP, 34.348984975465861_DP, 0.319288254087299_DP /), & 793 fx0(3) = (/ 37.537128437796000_DP, 87.173106479036008_DP, 11.289489669936000_DP /), & 794 bx0fx0(3) = (/ -0.142530524167984_DP, -0.171582499414508_DP, -0.021608710360898_DP /) 795 ! N.B.: A is expressed in Hartree 796 ! Q = SQRT(4*c - b^2) 797 ! tbQ = 2*b/Q 798 ! fx0 = X(x_0) = x_0^2 + b*x_0 + c 799 ! bx0fx0 = b*x_0/X(x_0) 800 ! 801 ! 802 ! coefficients for f(z), df/dz, ddf/ddz(0) 803 cfz = 2.0_DP**(4.0_DP/3.0_DP) - 2.0_DP 804 cfz1 = 1.0_DP / cfz 805 cfz2 = 4.0_DP/3.0_DP * cfz1 806 iddfz0 = 9.0_DP / 8.0_DP *cfz 807 ! 808 SQRTrs = SQRT(rs) 809 zeta3 = zeta**3 810 zeta4 = zeta3*zeta 811 trup = 1.0_DP + zeta 812 trdw = 1.0_DP - zeta 813 trup13 = trup**(1.0_DP/3.0_DP) 814 trdw13 = trdw**(1.0_DP/3.0_DP) 815 fz = cfz1 * (trup13*trup + trdw13*trdw - 2.0_DP) ! f(zeta) 816 dfz = cfz2 * (trup13 - trdw13) ! df / dzeta 817 ! 818 CALL padefit_ParSet2( sqrtrs, 1, ecP, vcP ) ! ecF = e_c Paramagnetic !<GPU:padefit_ParSet2=>padefit_ParSet2_d> 819 CALL padefit_ParSet2( sqrtrs, 2, ecF, vcF ) ! ecP = e_c Ferromagnetic !<GPU:padefit_ParSet2=>padefit_ParSet2_d> 820 CALL padefit_ParSet2( sqrtrs, 3, ac, dac ) ! ac = "spin stiffness" !<GPU:padefit_ParSet2=>padefit_ParSet2_d> 821 ! 822 ac = ac * iddfz0 823 dac = dac * iddfz0 824 De = ecF - ecP - ac ! e_c[F]-e_c[P]-alpha_c/(ddf/ddz(z=0)) 825 fzz4 = fz * zeta4 826 ec = ecP + ac * fz + De * fzz4 827 ! 828 dec1 = vcP + dac*fz + (vcF - vcP - dac) * fzz4 ! e_c - (r_s/3)*(de_c/dr_s) 829 dec2 = ac*dfz + De*(4.0_DP*zeta3*fz + zeta4*dfz) ! de_c/dzeta 830 ! 831 vc_up = dec1 + (1.0_DP - zeta)*dec2 ! v_c[s] = e_c - (r_s/3)*(de_c/dr_s) 832 vc_dw = dec1 - (1.0_DP + zeta)*dec2 ! +[sign(s)-zeta]*(de_c/dzeta) 833 ! 834END SUBROUTINE 835! 836! 837!----------------------------------------------------------------------- 838SUBROUTINE pw_spin( rs, zeta, ec, vc_up, vc_dw ) !<GPU:DEVICE> 839 !----------------------------------------------------------------------- 840 !! J.P. Perdew and Y. Wang, PRB 45, 13244 (1992). 841 ! 842 USE kinds, ONLY : DP 843 ! 844 IMPLICIT NONE 845 ! 846 REAL(DP), INTENT(IN) :: rs 847 !! Wigner-Seitz radius 848 REAL(DP), INTENT(IN) :: zeta 849 !! zeta = (rho_up - rho_dw)/rho_tot 850 REAL(DP), INTENT(OUT) :: ec 851 !! correlation energy 852 REAL(DP), INTENT(OUT) :: vc_up, vc_dw 853 !! correlation potential (up, down) 854 ! 855 ! ... local variables 856 ! 857 REAL(DP) :: rs12, rs32, rs2, zeta2, zeta3, zeta4, fz, dfz 858 REAL(DP) :: om, dom, olog, epwc, vpwc 859 REAL(DP) :: omp, domp, ologp, epwcp, vpwcp 860 REAL(DP) :: oma, doma, ologa, alpha, vpwca 861 ! 862 ! xc parameters, unpolarised 863 REAL(DP), PARAMETER :: a = 0.031091d0, a1 = 0.21370d0, b1 = 7.5957d0, b2 = & 864 3.5876d0, b3 = 1.6382d0, b4 = 0.49294d0, c0 = a, c1 = 0.046644d0, & 865 c2 = 0.00664d0, c3 = 0.01043d0, d0 = 0.4335d0, d1 = 1.4408d0 866 ! xc parameters, polarised 867 REAL(DP), PARAMETER :: ap = 0.015545d0, a1p = 0.20548d0, b1p = 14.1189d0, b2p & 868 = 6.1977d0, b3p = 3.3662d0, b4p = 0.62517d0, c0p = ap, c1p = & 869 0.025599d0, c2p = 0.00319d0, c3p = 0.00384d0, d0p = 0.3287d0, d1p & 870 = 1.7697d0 871 ! xc PARAMETERs, antiferro 872 REAL(DP), PARAMETER :: aa = 0.016887d0, a1a = 0.11125d0, b1a = 10.357d0, b2a = & 873 3.6231d0, b3a = 0.88026d0, b4a = 0.49671d0, c0a = aa, c1a = & 874 0.035475d0, c2a = 0.00188d0, c3a = 0.00521d0, d0a = 0.2240d0, d1a & 875 = 0.3969d0 876 REAL(DP), PARAMETER :: fz0 = 1.709921d0 877 ! 878 ! if (rs < 0.5d0) then 879 ! high density formula (not implemented) 880 ! 881 ! elseif (rs > 100.d0) then 882 ! low density formula (not implemented) 883 ! 884 ! else 885 ! interpolation formula 886 ! 887 zeta2 = zeta * zeta 888 zeta3 = zeta2 * zeta 889 zeta4 = zeta3 * zeta 890 rs12 = SQRT(rs) 891 rs32 = rs * rs12 892 rs2 = rs**2 893 ! 894 ! unpolarised 895 om = 2.d0 * a * (b1 * rs12 + b2 * rs + b3 * rs32 + b4 * rs2) 896 dom = 2.d0 * a * (0.5d0 * b1 * rs12 + b2 * rs + 1.5d0 * b3 * rs32 & 897 + 2.d0 * b4 * rs2) 898 olog = LOG(1.d0 + 1.0d0 / om) 899 epwc = - 2.d0 * a * (1.d0 + a1 * rs) * olog 900 vpwc = - 2.d0 * a * (1.d0 + 2.d0 / 3.d0 * a1 * rs) * olog - 2.d0 / & 901 3.d0 * a * (1.d0 + a1 * rs) * dom / (om * (om + 1.d0) ) 902 ! 903 ! polarized 904 omp = 2.d0 * ap * (b1p * rs12 + b2p * rs + b3p * rs32 + b4p * rs2) 905 domp = 2.d0 * ap * (0.5d0 * b1p * rs12 + b2p * rs + 1.5d0 * b3p * & 906 rs32 + 2.d0 * b4p * rs2) 907 ologp = LOG(1.d0 + 1.0d0 / omp) 908 epwcp = - 2.d0 * ap * (1.d0 + a1p * rs) * ologp 909 vpwcp = - 2.d0 * ap * (1.d0 + 2.d0 / 3.d0 * a1p * rs) * ologp - & 910 2.d0 / 3.d0 * ap * (1.d0 + a1p * rs) * domp / (omp * (omp + 1.d0)) 911 ! 912 ! antiferro 913 oma = 2.d0 * aa * (b1a * rs12 + b2a * rs + b3a * rs32 + b4a * rs2) 914 doma = 2.d0 * aa * ( 0.5d0 * b1a * rs12 + b2a * rs + 1.5d0 * b3a * & 915 rs32 + 2.d0 * b4a * rs2 ) 916 ologa = LOG( 1.d0 + 1.0d0/oma ) 917 alpha = 2.d0 * aa * (1.d0 + a1a*rs) * ologa 918 vpwca = + 2.d0 * aa * (1.d0 + 2.d0/3.d0 * a1a * rs) * ologa + & 919 2.d0 / 3.d0 * aa * (1.d0 + a1a*rs) * doma / (oma * (oma + 1.d0)) 920 ! 921 ! 922 fz = ( (1.d0 + zeta)**(4.d0 / 3.d0) + (1.d0 - zeta)**(4.d0 / & 923 3.d0) - 2.d0) / (2.d0** (4.d0 / 3.d0) - 2.d0) 924 dfz = ( (1.d0 + zeta)**(1.d0 / 3.d0) - (1.d0 - zeta)**(1.d0 / & 925 3.d0) ) * 4.d0 / (3.d0 * (2.d0** (4.d0 / 3.d0) - 2.d0) ) 926 ! 927 ! 928 ec = epwc + alpha * fz * (1.d0 - zeta4) / fz0 + (epwcp - epwc) & 929 * fz * zeta4 930 ! 931 vc_up = vpwc + vpwca * fz * (1.d0 - zeta4) / fz0 + (vpwcp - vpwc) & 932 * fz * zeta4 + (alpha / fz0 * (dfz * (1.d0 - zeta4) - 4.d0 * fz * & 933 zeta3) + (epwcp - epwc) * (dfz * zeta4 + 4.d0 * fz * zeta3) ) & 934 * (1.d0 - zeta) 935 ! 936 vc_dw = vpwc + vpwca * fz * (1.d0 - zeta4) / fz0 + (vpwcp - vpwc) & 937 * fz * zeta4 - (alpha / fz0 * (dfz * (1.d0 - zeta4) - 4.d0 * fz * & 938 zeta3) + (epwcp - epwc) * (dfz * zeta4 + 4.d0 * fz * zeta3) ) & 939 * (1.d0 + zeta) 940 ! 941 RETURN 942 ! 943END SUBROUTINE pw_spin 944! 945! 946!----------------------------------------------------------------------------- 947SUBROUTINE lsd_lyp( rho, zeta, elyp, vlyp_up, vlyp_dw ) !<GPU:DEVICE> 948 !==--------------------------------------------------------------== 949 !== C. LEE, W. YANG, AND R.G. PARR, PRB 37, 785 (1988) == 950 !== THIS IS ONLY THE LDA PART == 951 !==--------------------------------------------------------------== 952 ! 953 USE kinds, ONLY: DP 954 ! 955 IMPLICIT NONE 956 ! 957 REAL(DP), INTENT(IN) :: rho 958 !! total charge density 959 REAL(DP), INTENT(IN) :: zeta 960 !! zeta = (rho_up - rho_dw)/rho_tot 961 REAL(DP), INTENT(OUT) :: elyp 962 !! correlation energy 963 REAL(DP), INTENT(OUT) :: vlyp_up, vlyp_dw 964 !! correlation potential (up, down) 965 ! 966 ! ... local variables 967 ! 968 REAL(DP) :: ra,rb,rm3,dr,e1,or,dor,e2,de1a,de1b,de2a,de2b 969 REAL(DP), PARAMETER :: small=1.D-24, a=0.04918D0, b=0.132D0, & 970 c=0.2533D0, d=0.349D0, cf=2.87123400018819108D0 971 !==--------------------------------------------------------------== 972 ra = rho*0.5D0*(1.D0+zeta) 973 ra = MAX(ra,small) 974 rb = rho*0.5D0*(1.D0-zeta) 975 rb = MAX(rb,small) 976 ! 977 rm3= rho**(-1.D0/3.D0) 978 dr = (1.D0+d*rm3) 979 ! 980 e1 = 4.D0*a*ra*rb/rho/dr 981 or = EXP(-c*rm3)/dr*rm3**11.D0 982 dor= -1.D0/3.D0*rm3**4*or*(11.D0/rm3-c-d/dr) 983 e2 = 2.D0**(11.D0/3.D0)*cf*a*b*or*ra*rb*(ra**(8.d0/3.d0)+ rb**(8.d0/3.d0)) 984 ! 985 elyp = (-e1-e2)/rho 986 ! 987 de1a = -e1*(1.D0/3.D0*d*rm3**4/dr+1./ra-1./rho) 988 de1b = -e1*(1.D0/3.D0*d*rm3**4/dr+1./rb-1./rho) 989 de2a = -2.D0**(11.D0/3.D0)*cf*a*b*(dor*ra*rb*(ra**(8.d0/3.d0)+ & 990 rb**(8.d0/3.d0))+or*rb*(11.d0/3.d0*ra**(8.d0/3.d0)+ & 991 rb**(8.d0/3.d0))) 992 de2b = -2.D0**(11.D0/3.D0)*cf*a*b*(dor*ra*rb*(ra**(8.d0/3.d0)+ & 993 rb**(8.d0/3.d0))+or*ra*(11.d0/3.d0*rb**(8.d0/3.d0)+ & 994 ra**(8.d0/3.d0))) 995 ! 996 vlyp_up = de1a + de2a 997 vlyp_dw = de1b + de2b 998 !==--------------------------------------------------------------== 999 ! 1000 RETURN 1001 ! 1002END SUBROUTINE lsd_lyp 1003 1004END MODULE 1005