1! 2! Copyright (C) 2001-2012 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! 9!------------------------------------------------------------------------- 10! 11! META-GGA FUNCTIONALS 12! 13! Available functionals : 14! - TPSS (Tao, Perdew, Staroverov & Scuseria) 15! - M06L 16! Other options are available through libxc library (must be specified 17! in input). 18! 19!========================================================================= 20! 21MODULE metagga !<GPU:metagga=>metagga_gpu> 22! 23USE exch_lda, ONLY : slater !<GPU:slater=>slater_d,exch_lda=>exch_lda_gpu> 24USE corr_lda, ONLY : pw, pw_spin !<GPU:pw=>pw_d,pw_spin=>pw_spin_d,corr_lda=>corr_lda_gpu> 25USE corr_gga, ONLY : pbec, pbec_spin !<GPU:pbec=>pbec_d,pbec_spin=>pbec_spin_d,corr_gga=>corr_gga_gpu> 26! 27 CONTAINS 28! TPSS 29! 30!------------------------------------------------------------------------- 31SUBROUTINE tpsscxc( rho, grho, tau, sx, sc, v1x, v2x, v3x, v1c, v2c, v3c ) !<GPU:DEVICE> 32 !----------------------------------------------------------------------- 33 !! TPSS metaGGA corrections for exchange and correlation - Hartree a.u. 34 ! 35 !! Definition: E_x = \int E_x(rho,grho) dr 36 ! 37 USE kinds, ONLY : DP 38 ! 39 IMPLICIT NONE 40 ! 41 REAL(DP), INTENT(IN) :: rho 42 !! the charge density 43 REAL(DP), INTENT(IN) :: grho 44 !! grho = |\nabla rho|^2 45 REAL(DP), INTENT(IN) :: tau 46 !! kinetic energy density 47 REAL(DP), INTENT(OUT) :: sx 48 !! sx = E_x(rho,grho) 49 REAL(DP), INTENT(OUT) :: sc 50 !! sc = E_c(rho,grho) 51 REAL(DP), INTENT(OUT) :: v1x 52 !! v1x = D(E_x)/D(rho) 53 REAL(DP), INTENT(OUT) :: v2x 54 !! v2x = D(E_x)/D( D rho/D r_alpha ) / |\nabla rho| 55 REAL(DP), INTENT(OUT) :: v3x 56 !! v3x = D(E_x)/D(tau) 57 REAL(DP), INTENT(OUT) :: v1c 58 !! v1c = D(E_c)/D(rho) 59 REAL(DP), INTENT(OUT) :: v2c 60 !! v2c = D(E_c)/D( D rho/D r_alpha ) / |\nabla rho| 61 REAL(DP), INTENT(OUT) :: v3c 62 !! v3c = D(E_c)/D(tau) 63 ! 64 ! ... local variables 65 ! 66 REAL(DP), PARAMETER :: small = 1.E-10_DP 67 ! 68 IF (rho <= small) THEN 69 sx = 0.0_DP 70 v1x = 0.0_DP 71 v2x = 0.0_DP 72 sc = 0.0_DP 73 v1c = 0.0_DP 74 v2c = 0.0_DP 75 v3x = 0.0_DP 76 v3c = 0.0_DP 77 RETURN 78 ENDIF 79 ! 80 ! exchange 81 CALL metax( rho, grho, tau, sx, v1x, v2x, v3x ) !<GPU:metax=>metax_d> 82 ! correlation 83 CALL metac( rho, grho, tau, sc, v1c, v2c, v3c ) !<GPU:metac=>metac_d> 84 ! 85 RETURN 86 ! 87END SUBROUTINE tpsscxc 88! 89! 90!------------------------------------------------------------------------- 91SUBROUTINE metax( rho, grho2, tau, ex, v1x, v2x, v3x ) !<GPU:DEVICE> 92 !-------------------------------------------------------------------- 93 !! TPSS meta-GGA exchange potential and energy. 94 ! 95 !! NOTE: E_x(rho,grho) = rho\epsilon_x(rho,grho) ; 96 !! ex = E_x(rho,grho) NOT \epsilon_x(rho,grho) ; 97 !! v1x= D(E_x)/D(rho) ; 98 !! v2x= D(E_x)/D( D rho/D r_alpha ) / |\nabla rho| ; 99 !! v3x= D(E_x)/D( tau ) ; 100 !! tau is the kinetic energy density ; 101 !! the same applies to correlation terms ; 102 !! input grho2 is |\nabla rho|^2 . 103 ! 104 USE kinds, ONLY : DP 105 ! 106 IMPLICIT NONE 107 ! 108 REAL(DP), INTENT(IN) :: rho 109 !! the charge density 110 REAL(DP), INTENT(IN) :: grho2 111 !! grho2 = |\nabla rho|^2 112 REAL(DP), INTENT(IN) :: tau 113 !! kinetic energy density 114 REAL(DP), INTENT(OUT) :: ex 115 !! ex = E_x(rho,grho) 116 REAL(DP), INTENT(OUT) :: v1x 117 !! v1x = D(E_x)/D(rho) 118 REAL(DP), INTENT(OUT) :: v2x 119 !! v2x = D(E_x)/D( D rho/D r_alpha ) / |\nabla rho| 120 REAL(DP), INTENT(OUT) :: v3x 121 !! v3x = D(E_x)/D(tau) 122 ! 123 ! ... local variables 124 ! 125 REAL(DP) :: rs, vx_unif, ex_unif 126 ! ex_unif: lda \epsilon_x(rho) 127 ! ec_unif: lda \epsilon_c(rho) 128 REAL(DP), PARAMETER :: small=1.E-10_DP 129 REAL(DP), PARAMETER :: pi34=0.6203504908994_DP, third=1.0_DP/3.0_DP 130 ! fx=Fx(p,z) 131 ! fxp=d Fx / d p 132 ! fxz=d Fx / d z 133 REAL(DP) :: fx, f1x, f2x, f3x 134 ! 135 IF (ABS(tau) < small) THEN 136 ex = 0.0_DP 137 v1x = 0.0_DP 138 v2x = 0.0_DP 139 v3x = 0.0_DP 140 RETURN 141 ENDIF 142 ! 143 rs = pi34/rho**third 144 CALL slater( rs, ex_unif, vx_unif ) !<GPU:slater=>slater_d> 145 CALL metaFX( rho, grho2, tau, fx, f1x, f2x, f3x ) !<GPU:metaFX=>metaFX_d> 146 ! 147 ex = rho*ex_unif 148 v1x = vx_unif*fx + ex*f1x 149 v2x = ex*f2x 150 v3x = ex*f3x 151 ex = ex*fx 152 ! 153 RETURN 154 ! 155END SUBROUTINE metax 156! 157! 158!------------------------------------------------------------------ 159SUBROUTINE metac( rho, grho2, tau, ec, v1c, v2c, v3c ) !<GPU:DEVICE> 160 !-------------------------------------------------------------- 161 !! TPSS meta-GGA correlation energy and potentials. 162 ! 163 USE kinds, ONLY : DP 164 ! 165 IMPLICIT NONE 166 ! 167 REAL(DP), INTENT(IN) :: rho 168 !! the charge density 169 REAL(DP), INTENT(IN) :: grho2 170 !! grho2 = |\nabla rho|^2 171 REAL(DP), INTENT(IN) :: tau 172 !! kinetic energy density 173 REAL(DP), INTENT(OUT) :: ec 174 !! ec = E_c(rho,grho) 175 REAL(DP), INTENT(OUT) :: v1c 176 !! v1c = D(E_c)/D(rho) 177 REAL(DP), INTENT(OUT) :: v2c 178 !! v2c = D(E_c)/D( D rho/D r_alpha ) / |\nabla rho| 179 REAL(DP), INTENT(OUT) :: v3c 180 !! v3c = D(E_c)/D(tau) 181 ! 182 ! ... local variables 183 ! 184 REAL(DP) :: z, z2, tauw, ec_rev 185 REAL(DP) :: d1rev, d2rev, d3rev 186 ! d1ec= D ec_rev / D rho 187 ! d2ec= D ec_rev / D |D rho/ D r| / |\nabla rho| 188 ! d3ec= D ec_rev / D tau 189 REAL(DP) :: cf1, cf2, cf3 190 REAL(DP) :: v1c_pbe, v2c_pbe, ec_pbe 191 REAL(DP) :: v1c_sum(2), v2c_sum, ec_sum 192 ! 193 REAL(DP) :: rs 194 REAL(DP) :: ec_unif, vc_unif 195 REAL(DP) :: vc_unif_s(2) 196 ! 197 REAL(DP) :: rhoup, grhoup 198 ! 199 REAL(DP), PARAMETER :: small=1.0E-10_DP 200 REAL(DP), PARAMETER :: pi34=0.75_DP/3.141592653589793_DP, & 201 third=1.0_DP/3.0_DP 202 REAL(DP), PARAMETER :: dd=2.80_DP !in unit of Hartree^-1 203 REAL(DP), PARAMETER :: cab=0.53_DP, cabone=1.0_DP+cab 204 ! 205 IF (ABS(tau) < small) THEN 206 ec = 0.0_DP 207 v1c = 0.0_DP 208 v2c = 0.0_DP 209 v3c = 0.0_DP 210 RETURN 211 ENDIF 212 ! 213 rhoup = 0.5_DP*rho 214 grhoup = 0.5_DP*SQRT(grho2) 215 ! 216 IF (rhoup > small) THEN 217 ! 218 rs = (pi34/rhoup)**third 219 CALL pw_spin( rs, 1.0_DP, ec_unif, vc_unif_s(1), vc_unif_s(2) ) !<GPU:pw_spin=>pw_spin_d> 220 ! 221 IF (ABS(grhoup) > small) THEN 222 ! zeta=1.0_DP-small to avoid pow_e of 0 in pbec_spin 223 CALL pbec_spin( rhoup, 1.0_DP-small, grhoup**2, 1, & !<GPU:pbec_spin=>pbec_spin_d> 224 ec_sum, v1c_sum(1), v1c_sum(2), v2c_sum ) 225 ELSE 226 ec_sum = 0.0_DP 227 v1c_sum = 0.0_DP 228 v2c_sum = 0.0_DP 229 ENDIF 230 ec_sum = ec_sum/rhoup + ec_unif 231 v1c_sum(1) = (v1c_sum(1) + vc_unif_s(1)-ec_sum)/rho !rho, not rhoup 232 v2c_sum = v2c_sum/(2.0_DP*rho) 233 ELSE 234 ec_sum = 0.0_DP 235 v1c_sum = 0.0_DP 236 v2c_sum = 0.0_DP 237 ENDIF 238 ! 239 rs = (pi34/rho)**third 240 CALL pw( rs, 1, ec_unif, vc_unif ) !<GPU:pw=>pw_d> 241 ! 242 ! ... PBE correlation energy and potential: 243 ! ec_pbe=rho*H, not rho*(epsion_c_uinf + H) 244 ! v1c_pbe=D (rho*H) /D rho 245 ! v2c_pbe= for rho, 2 for 246 CALL pbec( rho, grho2, 1, ec_pbe, v1c_pbe, v2c_pbe ) !<GPU:pbec=>pbec_d> 247 ! 248 ec_pbe = ec_pbe/rho+ec_unif 249 v1c_pbe = (v1c_pbe+vc_unif-ec_pbe)/rho 250 v2c_pbe = v2c_pbe/rho 251 ! 252 IF (ec_sum < ec_pbe) THEN 253 ec_sum = ec_pbe 254 v1c_sum(1) = v1c_pbe 255 v2c_sum = v2c_pbe 256 ENDIF 257 ! 258 tauw = 0.1250_DP * grho2/rho 259 z = tauw/tau 260 z2 = z*z 261 ! 262 ec_rev = ec_pbe*(1+cab*z2)-cabone*z2*ec_sum 263 ! 264 d1rev = v1c_pbe + (cab*v1c_pbe-cabone*v1c_sum(1))*z2 & 265 - (ec_pbe*cab - ec_sum*cabone)*2.0_DP*z2/rho 266 d2rev = v2c_pbe + (cab*v2c_pbe-cabone*v2c_sum)*z2 & 267 + (ec_pbe*cab - ec_sum*cabone)*4.0_DP*z2/grho2 268 d3rev = -(ec_pbe*cab - ec_sum*cabone)*2.0_DP*z2/tau 269 ! 270 cf1 = 1.0_DP + dd*ec_rev*z2*z 271 cf2 = rho*(1.0_DP+2.0_DP*z2*z*dd*ec_rev) 272 cf3 = ec_rev*ec_rev*3.0_DP*dd*z2*z 273 v1c = ec_rev*cf1 + cf2*d1rev-cf3 274 ! 275 cf3 = cf3*rho 276 v2c = cf2*d2rev + cf3*2.0_DP/grho2 277 v3c = cf2*d3rev - cf3/tau 278 ! 279 ec = rho*ec_rev*(1.0_DP+dd*ec_rev*z2*z) !-rho*ec_unif(1) 280 !v1c = v1c - vc_unif(1) 281 ! 282 RETURN 283 ! 284END SUBROUTINE metac 285! 286! 287!------------------------------------------------------------------------- 288SUBROUTINE metaFX( rho, grho2, tau, fx, f1x, f2x, f3x ) !<GPU:DEVICE> 289 !------------------------------------------------------------------------- 290 ! 291 USE kinds, ONLY : DP 292 ! 293 IMPLICIT NONE 294 ! 295 REAL(DP), INTENT(IN) :: rho 296 !! charge density 297 REAL(DP), INTENT(IN) :: grho2 298 !! square of gradient of rho 299 REAL(DP), INTENT(IN) :: tau 300 !! kinetic energy density 301 REAL(DP), INTENT(OUT) :: fx 302 !! fx = Fx(p,z) 303 REAL(DP), INTENT(OUT) :: f1x 304 !! f1x=D (Fx) / D rho 305 REAL(DP), INTENT(OUT) :: f2x 306 !! f2x=D (Fx) / D ( D rho/D r_alpha) /|nabla rho| 307 REAL(DP), INTENT(OUT) :: f3x 308 !! f3x=D (Fx) / D tau 309 ! 310 ! ... local variables 311 ! 312 REAL(DP) x, p, z, qb, al, localdp, dz 313 REAL(DP) dfdx, dxdp, dxdz, dqbdp, daldp, dqbdz, daldz 314 REAL(DP) fxp, fxz ! fxp =D fx /D p 315 REAL(DP) tauw, tau_unif 316 ! work variables 317 REAL(DP) xf1, xf2 318 REAL(DP) xfac1, xfac2, xfac3, xfac4, xfac5, xfac6, xfac7, z2 319 ! 320 REAL(DP), PARAMETER :: pi=3.141592653589793_DP 321 REAL(DP), PARAMETER :: THRD=0.3333333333333333_DP 322 REAL(DP), PARAMETER :: ee=1.537_DP 323 REAL(DP), PARAMETER :: cc=1.59096_DP 324 REAL(DP), PARAMETER :: kk=0.804_DP 325 REAL(DP), PARAMETER :: bb=0.40_DP 326 REAL(DP), PARAMETER :: miu=0.21951_DP 327 REAL(DP), PARAMETER :: fac1=9.57078000062731_DP !fac1=(3*pi^2)^(2/3) 328 REAL(DP), PARAMETER :: small=1.0E-6_DP 329 ! 330 tauw = 0.125_DP*grho2/rho 331 z = tauw/tau 332 ! 333 p = SQRT(grho2)/rho**THRD/rho 334 p = p*p/(fac1*4.0_DP) 335 tau_unif = 0.3_DP*fac1*rho**(5.0_DP/3.0_DP) 336 al = (tau-tauw)/tau_unif 337 al = ABS(al) !make sure al is always .gt. 0.0_DP 338 qb = 0.45_DP*(al-1.0_DP)/SQRT(1.0_DP+bb*al*(al-1.0_DP)) 339 qb = qb+2.0_DP*THRD*p 340 ! 341 ! calculate x(p,z) and fx 342 z2 = z*z 343 xf1 = 10.0_DP/81.0_DP 344 xfac1 = xf1+cc*z2/(1+z2)**2.0_DP 345 xfac2 = 146.0_DP/2025.0_DP 346 xfac3 = SQRT(0.5_DP*(0.36_DP*z2+p*p)) 347 xfac4 = xf1*xf1/kk 348 xfac5 = 2.0_DP*SQRT(ee)*xf1*0.36_DP 349 xfac6 = xfac1*p+xfac2*qb**2.0_DP-73.0_DP/405.0_DP*qb*xfac3 350 xfac6 = xfac6+xfac4*p**2.0_DP+xfac5*z2+ee*miu*p**3.0_DP 351 xfac7 = (1+SQRT(ee)*p) 352 x = xfac6/(xfac7*xfac7) 353 ! fx=kk-kk/(1.0_DP+x/kk) 354 fx = 1.0_DP + kk-kk/(1.0_DP+x/kk) 355 ! 356 ! calculate the derivatives of fx w.r.t p and z 357 dfdx = (kk/(kk+x))**2.0_DP 358 daldp = 5.0_DP*THRD*(tau/tauw-1.0_DP) 359 ! 360 ! daldz=-0.50_DP*THRD* 361 ! * (tau/(2.0_DP*fac1*rho**THRD*0.1250_DP*sqrt(grho2)))**2.0_DP 362 daldz = -5.0_DP*THRD*p/z2 363 dqbdz = 0.45_DP*(0.50_DP*bb*(al-1.0_DP)+1.0_DP) 364 dqbdz = dqbdz/(1.0_DP+bb*al*(al-1.0_DP))**1.5_DP 365 ! 366 dqbdp = dqbdz*daldp+2.0_DP*THRD 367 dqbdz = dqbdz*daldz 368 ! 369 ! calculate dx/dp 370 xf1 = 73.0_DP/405.0_DP/xfac3*0.50_DP*qb 371 xf2 = 2.0_DP*xfac2*qb-73.0_DP/405.0_DP*xfac3 372 ! 373 dxdp = -xf1*p 374 dxdp = dxdp+xfac1+xf2*dqbdp 375 dxdp = dxdp+2.0_DP*xfac4*p 376 dxdp = dxdp+3.0_DP*ee*miu*p*p 377 dxdp = dxdp/(xfac7*xfac7)-2.0_DP*x*SQRT(ee)/xfac7 378 ! 379 ! dx/dz 380 dxdz = -xf1*0.36_DP*z 381 xfac1 = cc*2.0_DP*z*(1-z2)/(1+z2)**3.0_DP 382 dxdz = dxdz+xfac1*p+xf2*dqbdz 383 dxdz = dxdz+xfac5*2.0_DP*z 384 dxdz = dxdz/(xfac7*xfac7) 385 ! 386 fxp = dfdx*dxdp 387 fxz = dfdx*dxdz 388 ! 389 ! calculate f1x 390 localdp = -8.0_DP*THRD*p/rho ! D p /D rho 391 dz = -z/rho ! D z /D rho 392 f1x = fxp*localdp+fxz*dz 393 ! 394 ! f2x 395 localdp = 2.0_DP/(fac1*4.0_DP*rho**(8.0_DP/3.0_DP)) 396 dz = 2.0_DP*0.125_DP/(rho*tau) 397 f2x = fxp*localdp + fxz*dz 398 ! 399 ! f3x 400 localdp = 0.0_DP 401 dz = -z/tau 402 f3x = fxz*dz 403 ! 404 RETURN 405 ! 406END SUBROUTINE metaFX 407! 408! 409!--------------------------------------------------------------------------- 410SUBROUTINE tpsscx_spin( rhoup, rhodw, grhoup2, grhodw2, tauup, taudw, sx, & !<GPU:DEVICE> 411 v1xup, v1xdw, v2xup, v2xdw, v3xup, v3xdw ) 412 !----------------------------------------------------------------------- 413 !! TPSS metaGGA for exchange - Hartree a.u. 414 ! 415 USE kinds, ONLY : DP 416 ! 417 IMPLICIT NONE 418 ! 419 REAL(DP), INTENT(IN) :: rhoup 420 !! up charge 421 REAL(DP), INTENT(IN) :: rhodw 422 !! down charge 423 REAL(DP), INTENT(IN) :: grhoup2 424 !! up gradient of the charge 425 REAL(DP), INTENT(IN) :: grhodw2 426 !! down gradient of the charge 427 REAL(DP), INTENT(IN) :: tauup 428 !! up kinetic energy density 429 REAL(DP), INTENT(IN) :: taudw 430 !! down kinetic energy density 431 REAL(DP), INTENT(OUT) :: sx 432 !! exchange energy 433 REAL(DP), INTENT(OUT) :: v1xup 434 !! derivatives of exchange wr. rho - up 435 REAL(DP), INTENT(OUT) :: v1xdw 436 !! derivatives of exchange wr. rho - down 437 REAL(DP), INTENT(OUT) :: v2xup 438 !! derivatives of exchange wr. grho - up 439 REAL(DP), INTENT(OUT) :: v2xdw 440 !! derivatives of exchange wr. grho - down 441 REAL(DP), INTENT(OUT) :: v3xup 442 !! derivatives of exchange wr. tau - up 443 REAL(DP), INTENT(OUT) :: v3xdw 444 !! derivatives of exchange wr. tau - down 445 ! 446 ! ... local variables 447 ! 448 REAL(DP), PARAMETER :: small = 1.E-10_DP 449 REAL(DP) :: rho, sxup, sxdw 450 ! 451 ! exchange 452 rho = rhoup + rhodw 453 IF (rhoup>small .AND. SQRT(ABS(grhoup2))>small & 454 .AND. ABS(tauup) > small) THEN 455 CALL metax( 2.0_DP*rhoup, 4.0_DP*grhoup2, & !<GPU:metax=>metax_d> 456 2.0_DP*tauup, sxup, v1xup, v2xup, v3xup ) 457 ELSE 458 sxup = 0.0_DP 459 v1xup = 0.0_DP 460 v2xup = 0.0_DP 461 v3xup = 0.0_DP 462 ENDIF 463 ! 464 IF (rhodw > small .AND. SQRT(ABS(grhodw2)) > small & 465 .AND. ABS(taudw) > small) THEN 466 CALL metax( 2.0_DP*rhodw, 4.0_DP*grhodw2, & !<GPU:metax=>metax_d> 467 2.0_DP*taudw, sxdw, v1xdw, v2xdw, v3xdw ) 468 ELSE 469 sxdw = 0.0_DP 470 v1xdw = 0.0_DP 471 v2xdw = 0.0_DP 472 v3xdw = 0.0_DP 473 ENDIF 474 ! 475 sx = 0.5_DP*(sxup+sxdw) 476 v2xup = 2.0_DP*v2xup 477 v2xdw = 2.0_DP*v2xdw 478 ! 479 RETURN 480 ! 481END SUBROUTINE tpsscx_spin 482! 483! 484!--------------------------------------------------------------------------- 485SUBROUTINE tpsscc_spin( rho, zeta, grhoup, grhodw, & !<GPU:DEVICE> 486 atau, sc, v1cup, v1cdw, v2cup, v2cdw, v3cup, v3cdw ) 487 !-------------------------------------------------------------------------- 488 !! TPSS metaGGA for correlations - Hartree a.u. 489 ! 490 USE kinds, ONLY : DP 491 ! 492 IMPLICIT NONE 493 ! 494 REAL(DP), INTENT(IN) :: rho 495 !! the total charge 496 REAL(DP), INTENT(IN) :: zeta 497 !! the magnetization 498 REAL(DP), INTENT(IN) :: atau 499 !! the total kinetic energy density 500 REAL(DP), INTENT(IN) :: grhoup(3) 501 !! the gradient of the charge - up 502 REAL(DP), INTENT(IN) :: grhodw(3) 503 !! the gradient of the charge - down 504 REAL(DP), INTENT(OUT) :: sc 505 !! correlation energy 506 REAL(DP), INTENT(OUT) :: v1cup 507 !! derivatives of correlation wr. rho - up 508 REAL(DP), INTENT(OUT) :: v1cdw 509 !! derivatives of correlation wr. rho - down 510 REAL(DP), INTENT(OUT) :: v2cup(3) 511 !! derivatives of correlation wr. grho - up 512 REAL(DP), INTENT(OUT) :: v2cdw(3) 513 !! derivatives of correlation wr. grho - down 514 REAL(DP), INTENT(OUT) :: v3cup 515 !! derivatives of correlation wr. tau - up 516 REAL(DP), INTENT(OUT) :: v3cdw 517 !! derivatives of correlation wr. tau - down 518 ! 519 ! ... local variables 520 ! 521 REAL(DP) :: grho_vec(3) 522 REAL(DP) :: v3c, grho !grho=grho2 523 INTEGER :: ipol 524 REAL(DP), PARAMETER :: small = 1.E-10_DP 525 ! 526 ! vector 527 grho_vec = grhoup + grhodw 528 grho = 0.0_DP 529 ! 530 DO ipol = 1, 3 531 grho = grho + grho_vec(ipol)**2 532 ENDDO 533 ! 534 IF (rho <= small .OR. ABS(zeta) > 1.0_DP .OR. SQRT(ABS(grho)) <= small & 535 .OR. ABS(atau) < small ) THEN 536 ! 537 sc = 0.0_DP 538 v1cup = 0.0_DP 539 v1cdw = 0.0_DP 540 ! 541 v2cup(:) = 0.0_DP 542 v2cdw(:) = 0.0_DP 543 ! 544 v3cup = 0.0_DP 545 v3cdw = 0.0_DP 546 ! 547 v3c = 0.0_DP 548 ELSE 549 CALL metac_spin( rho, zeta, grhoup, grhodw, & !<GPU:metac_spin=>metac_spin_d> 550 atau, sc, v1cup, v1cdw, v2cup, v2cdw, v3c ) 551 ENDIF 552 ! 553 v3cup = v3c 554 v3cdw = v3c 555 ! 556 RETURN 557 ! 558END SUBROUTINE tpsscc_spin 559! 560! 561!--------------------------------------------------------------- 562SUBROUTINE metac_spin( rho, zeta, grhoup, grhodw, & !<GPU:DEVICE> 563 tau, sc, v1up, v1dw, v2up, v2dw, v3 ) 564 !--------------------------------------------------------------- 565 USE kinds, ONLY : DP 566 ! 567 IMPLICIT NONE 568 ! 569 REAL(DP), INTENT(IN) :: rho 570 !! the total charge 571 REAL(DP), INTENT(IN) :: zeta 572 !! the magnetization 573 REAL(DP), INTENT(IN) :: grhoup(3) 574 !! the gradient of the charge - up 575 REAL(DP), INTENT(IN) :: grhodw(3) 576 !! the gradient of the charge - down 577 REAL(DP), INTENT(IN) :: tau 578 !! the total kinetic energy density 579 REAL(DP), INTENT(OUT) :: sc 580 !! correlation energy 581 REAL(DP), INTENT(OUT) :: v1up 582 !! derivatives of correlation wr. rho - up 583 REAL(DP), INTENT(OUT) :: v1dw 584 !! derivatives of correlation wr. rho - down 585 REAL(DP), INTENT(OUT) :: v2up(3) 586 !! derivatives of correlation wr. grho - up 587 REAL(DP), INTENT(OUT) :: v2dw(3) 588 !! derivatives of correlation wr. grho - down 589 REAL(DP), INTENT(OUT) :: v3 590 !! derivatives of correlation wr. tau 591 ! 592 ! ... local variables 593 ! 594 REAL(DP) :: rhoup, rhodw, tauw, grhovec(3), grho2, grho, & 595 grhoup2, grhodw2 596 ! 597 REAL(DP) :: rs, ec_u, vc_u(2), v1_0v(2), v1_pbe(2) 598 ! 599 !grhovec vector gradient of rho 600 !grho mod of gradient of rho 601 !REAL(DP) :: ec_u, vcup_u, vcdw_u 602 REAL(DP) :: ec_pbe, v1up_pbe, v1dw_pbe, v2up_pbe(3), v2dw_pbe(3) 603 REAL(DP) :: ecup_0, v1up_0, v2up_0(3), v2_tmp 604 REAL(DP) :: ecdw_0, v1dw_0, v2dw_0(3) 605 REAL(DP) :: ec_rev, cab, aa, bb, aa2 606 ! 607 REAL(DP) :: z2, z, ca0, dca0da, dcabda, dcabdb 608 REAL(DP) :: term(3), term1, term2, term3 609 ! 610 REAL(DP) :: drev1up, drev1dw, drev2up(3), drev2dw(3), drev3 611 REAL(DP) :: sum, dsum1up, dsum1dw, dsum2up(3), dsum2dw(3) 612 ! 613 REAL(DP) :: dcab1up, dcab1dw, dcab2up(3), dcab2dw(3) 614 REAL(DP) :: db1up, db1dw, db2up(3), db2dw(3) 615 REAL(DP) :: da1up, da1dw 616 REAL(DP) :: ecup_til, ecdw_til 617 ! 618 REAL(DP) :: v1up_uptil, v1up_dwtil, v2up_uptil(3), v2up_dwtil(3) 619 REAL(DP) :: v1dw_uptil, v1dw_dwtil, v2dw_uptil(3), v2dw_dwtil(3) 620 ! 621 REAL(DP), PARAMETER :: small=1.0E-10_DP, & 622 fac=3.09366772628013593097_DP**2 623 ! fac = (3*PI**2)**(2/3) 624 REAL(DP), PARAMETER :: pi34=0.75_DP/3.141592653589793_DP, & 625 p43=4.0_DP/3.0_DP, third=1.0_DP/3.0_DP 626 INTEGER :: ipol 627 ! 628 rhoup = (1+zeta)*0.5_DP*rho 629 rhodw = (1-zeta)*0.5_DP*rho 630 grho2 = 0.0_DP 631 grhoup2 = 0.0_DP 632 grhodw2 = 0.0_DP 633 ! 634 DO ipol = 1, 3 635 grhovec(ipol) = grhoup(ipol) + grhodw(ipol) 636 grho2 = grho2 + grhovec(ipol)**2 637 grhoup2 = grhoup2 + grhoup(ipol)**2 638 grhodw2 = grhodw2 + grhodw(ipol)**2 639 ENDDO 640 ! 641 grho = SQRT(grho2) 642 ! 643 IF (rho > small) THEN 644 ! 645 v2_tmp = 0.0_DP 646 rs = (pi34/rho)**third 647 CALL pw_spin( rs, zeta, ec_u, vc_u(1), vc_u(2) ) !<GPU:pw_spin=>pw_spin_d> 648 ! 649 IF ((ABS(grho) > small) .AND. (zeta <= 1.0_DP)) THEN 650 CALL pbec_spin( rho, zeta, grho2, 1, ec_pbe, v1_pbe(1), v1_pbe(2), v2_tmp ) !<GPU:pbec_spin=>pbec_spin_d> 651 v1up_pbe=v1_pbe(1) 652 v1dw_pbe=v1_pbe(2) 653 ELSE 654 ec_pbe = 0.0_DP 655 v1up_pbe = 0.0_DP 656 v1dw_pbe = 0.0_DP 657 v2up_pbe = 0.0_DP 658 ENDIF 659 ! 660 ec_pbe = ec_pbe/rho+ec_u 661 ! v1xx_pbe = D_epsilon_c/ D_rho_xx :xx= up, dw 662 v1up_pbe = (v1up_pbe+vc_u(1)-ec_pbe)/rho 663 v1dw_pbe = (v1dw_pbe+vc_u(2)-ec_pbe)/rho 664 ! v2xx_pbe = (D_Ec / D grho)/rho = (D_Ec/ D|grho|/|grho|)*grho/rho 665 v2up_pbe = v2_tmp/rho*grhovec 666 ! v2dw === v2up for PBE 667 v2dw_pbe = v2up_pbe 668 ELSE 669 ec_pbe = 0.0_DP 670 v1up_pbe = 0.0_DP 671 v1dw_pbe = 0.0_DP 672 v2up_pbe = 0.0_DP 673 v2dw_pbe = 0.0_DP 674 ENDIF 675 ! 676 ! ec_pbe(rhoup,0,grhoup,0) 677 IF (rhoup > small) THEN 678 v2_tmp = 0.0_DP 679 ! 680 rs = (pi34/rhoup)**third 681 CALL pw_spin( rs, 1.0_DP, ec_u, vc_u(1), vc_u(2) ) !<GPU:pw_spin=>pw_spin_d> 682 ! 683 IF (SQRT(grhoup2) > small) THEN 684 CALL pbec_spin( rhoup, 1.0_DP-small, grhoup2, 1, & !<GPU:pbec_spin=>pbec_spin_d> 685 ecup_0, v1_0v(1), v1_0v(2), v2_tmp ) 686 v1up_0 = v1_0v(1) 687 v1dw_0 = v1_0v(2) 688 ELSE 689 ecup_0 = 0.0_DP 690 v1up_0 = 0.0_DP 691 v2up_0 = 0.0_DP 692 ENDIF 693 ! 694 ecup_0 = ecup_0/rhoup + ec_u 695 v1up_0 = (v1up_0 + vc_u(1)-ecup_0)/rhoup 696 v2up_0 = v2_tmp/rhoup*grhoup 697 ELSE 698 ecup_0 = 0.0_DP 699 v1up_0 = 0.0_DP 700 v2up_0 = 0.0_DP 701 ENDIF 702 ! 703 IF (ecup_0 > ec_pbe) THEN 704 ecup_til = ecup_0 705 v1up_uptil = v1up_0 706 v2up_uptil = v2up_0 707 v1up_dwtil = 0.0_DP 708 v2up_dwtil = 0.0_DP 709 ELSE 710 ecup_til = ec_pbe 711 v1up_uptil = v1up_pbe 712 v1up_dwtil = v1dw_pbe 713 v2up_uptil = v2up_pbe 714 v2up_dwtil = v2up_pbe 715 ENDIF 716 ! ec_pbe(rhodw,0,grhodw,0) 717 ! zeta = 1.0_DP 718 IF (rhodw > small) THEN 719 v2_tmp = 0.0_DP 720 ! 721 rs = (pi34/rhodw)**third 722 CALL pw_spin( rs, -1.0_DP, ec_u, vc_u(1), vc_u(2) ) !<GPU:pw_spin=>pw_spin_d> 723 ! 724 IF (SQRT(grhodw2) > small) THEN 725 ! 726 CALL pbec_spin( rhodw, -1.0_DP+small, grhodw2, 1, & !<GPU:pbec_spin=>pbec_spin_d> 727 ecdw_0, v1_0v(1), v1_0v(2), v2_tmp ) 728 ! 729 v1up_0 = v1_0v(1) 730 v1dw_0 = v1_0v(2) 731 ! 732 ELSE 733 ecdw_0 = 0.0_DP 734 v1dw_0 = 0.0_DP 735 v2dw_0 = 0.0_DP 736 ENDIF 737 ! 738 ecdw_0 = ecdw_0/rhodw + ec_u 739 v1dw_0 = (v1dw_0 + vc_u(2)-ecdw_0)/rhodw 740 v2dw_0 = v2_tmp/rhodw*grhodw 741 ELSE 742 ecdw_0 = 0.0_DP 743 v1dw_0 = 0.0_DP 744 v2dw_0 = 0.0_DP 745 ENDIF 746 ! 747 IF (ecdw_0 > ec_pbe) THEN 748 ecdw_til = ecdw_0 749 v1dw_dwtil = v1dw_0 750 v2dw_dwtil = v2dw_0 751 v1dw_uptil = 0.0_DP 752 v2dw_uptil = 0.0_DP 753 ELSE 754 ecdw_til = ec_pbe 755 v1dw_dwtil = v1dw_pbe 756 v2dw_dwtil = v2dw_pbe 757 v1dw_uptil = v1up_pbe 758 v2dw_uptil = v2dw_pbe 759 ENDIF 760 !cccccccccccccccccccccccccccccccccccccccccc-------checked 761 sum = (rhoup*ecup_til+rhodw*ecdw_til)/rho 762 dsum1up = (ecup_til-ecdw_til)*rhodw/rho**2 & 763 + (rhoup*v1up_uptil + rhodw*v1dw_uptil)/rho 764 dsum1dw = (ecdw_til-ecup_til)*rhoup/rho**2 & 765 + (rhodw*v1dw_dwtil + rhoup*v1up_dwtil)/rho 766 ! vector 767 dsum2up = (rhoup*v2up_uptil + rhodw*v2dw_uptil)/rho 768 dsum2dw = (rhodw*v2dw_dwtil + rhoup*v2up_dwtil)/rho 769 !ccccccccccccccccccccccccccccccccccccccccc---------checked 770 aa = zeta 771 ! bb = (rho*(grhoup-grhodw) - (rhoup-rhodw)*grho)**2 & 772 ! /(4.0_DP*fac*rho**(14.0_DP/3.0_DP)) 773 bb = 0.0_DP 774 DO ipol = 1, 3 775 term(ipol) = rhodw*grhoup(ipol)-rhoup*grhodw(ipol) 776 bb = bb + term(ipol)**2 777 ENDDO 778 ! vector 779 term = term/(fac*rho**(14.0_DP/3.0_DP)) 780 bb = bb/(fac*rho**(14.0_DP/3.0_DP)) 781 ! bb = (rhodw*grhoup-rhoup*grhodw)**2/fac*rho**(-14.0_DP/3.0_DP) 782 aa2 = aa*aa 783 ca0 = 0.53_DP+aa2*(0.87_DP+aa2*(0.50_DP+aa2*2.26_DP)) 784 dca0da = aa*(1.74_DP+aa2*(2.0_DP+aa2*13.56_DP)) 785 ! 786 IF (ABS(aa) <= 1.0_DP-small) THEN 787 term3 = (1.0_DP+aa)**(-p43) + (1.0_DP-aa)**(-p43) 788 term1 = (1.0_DP+bb*0.50_DP*term3) 789 term2 = (1.0_DP+aa)**(-7.0_DP/3.0_DP) + (1.0_DP-aa)**(-7.0_DP/3.0_DP) 790 cab = ca0/term1**4 791 dcabda = (dca0da/ca0 + 8.0_DP/3.0_DP*bb*term2/term1)*cab 792 dcabdb = -2.0_DP*cab*term3/term1 793 ELSE 794 cab = 0.0_DP 795 dcabda = 0.0_DP 796 dcabdb = 0.0_DP 797 ENDIF 798 ! 799 da1up = 2.0_DP*rhodw/rho**2 800 da1dw = -2.0_DP*rhoup/rho**2 801 db1up = -2.0_DP*(grhodw(1)*term(1)+grhodw(2)*term(2)+grhodw(3)*term(3)) & 802 -14.0_DP/3.0_DP*bb/rho 803 db1dw = 2.0_DP*(grhoup(1)*term(1)+grhoup(2)*term(2)+grhoup(3)*term(3)) & 804 -14.0_DP/3.0_DP*bb/rho 805 !vector, not scalar 806 db2up = term*rhodw*2.0_DP 807 db2dw = -term*rhoup*2.0_DP 808 ! 809 dcab1up = dcabda*da1up + dcabdb*db1up 810 dcab1dw = dcabda*da1dw + dcabdb*db1dw 811 !vector, not scalar 812 dcab2up = dcabdb*db2up 813 dcab2dw = dcabdb*db2dw 814 !cccccccccccccccccccccccccccccccccccccccccccccccccccccc------checked 815 tauw = 0.1250_DP*grho2/rho 816 z = tauw/tau 817 z2 = z*z 818 ! 819 term1 = 1.0_DP+cab*z2 820 term2 = (1.0_DP+cab)*z2 821 ec_rev = ec_pbe*term1-term2*sum 822 ! 823 drev1up = v1up_pbe*term1 & 824 + ec_pbe*(z2*dcab1up - 2.0_DP*cab*z2/rho) & 825 + (2.0_DP*term2/rho - z2*dcab1up)*sum & 826 - term2*dsum1up 827 ! 828 drev1dw = v1dw_pbe*term1 & 829 + ec_pbe*(z2*dcab1dw - 2.0_DP*cab*z2/rho) & 830 + (2.0_DP*term2/rho - z2*dcab1dw)*sum & 831 - term2*dsum1dw 832 ! 833 ! vector, not scalar 834 drev2up = v2up_pbe*term1 & 835 + ec_pbe*(z2*dcab2up+0.5_DP*cab*z/(rho*tau)*grhovec) & 836 - (term2*4.0_DP/grho2*grhovec + z2*dcab2up)*sum & 837 - term2*dsum2up 838 ! 839 drev2dw = v2dw_pbe*term1 & 840 + ec_pbe*(z2*dcab2dw+0.5_DP*cab*z/(rho*tau)*grhovec) & 841 - (term2*4.0_DP/grho2*grhovec + z2*dcab2dw)*sum & 842 - term2*dsum2dw 843 ! 844 drev3 = ((1.0_DP+cab)*sum-ec_pbe*cab)*2.0_DP*z2/tau 845 !ccccccccccccccccccccccccccccccccccccccccccccccccccc----checked 846 term1 = ec_rev*(1.0_DP+2.8_DP*ec_rev*z2*z) 847 term2 = (1.0_DP+5.6_DP*ec_rev*z2*z)*rho 848 term3 = -8.4_DP*ec_rev*ec_rev*z2*z 849 ! 850 v1up = term1 + term2*drev1up + term3 851 v1dw = term1 + term2*drev1dw + term3 852 ! 853 term3 = term3*rho 854 v3 = term2*drev3 + term3/tau 855 ! 856 term3 = -2.0_DP*term3/grho2 !grho/|grho|^2 = 1/grho 857 v2up = term2*drev2up + term3*grhovec 858 v2dw = term2*drev2dw + term3*grhovec 859 ! 860 ! call pw_spin((pi34/rho)**third,zeta,ec_u,vcup_u,vcdw_u) 861 sc = rho*ec_rev*(1.0_DP+2.8_DP*ec_rev*z2*z) !-rho*ec_u 862 ! v1up=v1up-vcup_u 863 ! v1dw=v1dw-vcdw_u 864 ! 865 RETURN 866 ! 867END SUBROUTINE metac_spin 868! 869! END TPSSS 870!========================================================================= 871! 872! --- M06L --- 873! 874! input: - rho 875! - grho2=|\nabla rho|^2 876! - tau = the kinetic energy density 877! It is defined as summ_i( |nabla phi_i|**2 ) 878! 879! definition: E_x = \int ex dr 880! 881! output: ex (rho, grho, tau) 882! v1x= D(E_x)/D(rho) 883! v2x= D(E_x)/D( D rho/D r_alpha ) / |\nabla rho| 884! ( v2x = 1/|grho| * dsx / d|grho| = 2 * dsx / dgrho2 ) 885! v3x= D(E_x)/D(tau) 886! 887! ec, v1c, v2c, v3c as above for correlation 888! 889!------------------------------------------------------------------------- 890SUBROUTINE m06lxc( rho, grho2, tau, ex, ec, v1x, v2x, v3x, v1c, v2c, v3c ) !<GPU:DEVICE> 891 !----------------------------------------------------------------------- 892 ! 893 USE kinds, ONLY : DP 894 ! 895 IMPLICIT NONE 896 ! 897 REAL(DP), INTENT(IN) :: rho 898 !! the total charge 899 REAL(DP), INTENT(IN) :: grho2 900 !! square of gradient of rho 901 REAL(DP), INTENT(IN) :: tau 902 !! the total kinetic energy density 903 REAL(DP), INTENT(OUT) :: ex 904 !! exchange energy 905 REAL(DP), INTENT(OUT) :: ec 906 !! correlation energy 907 REAL(DP), INTENT(OUT) :: v1x 908 !! derivatives of exchange wr. rho 909 REAL(DP), INTENT(OUT) :: v2x 910 !! derivatives of exchange wr. grho 911 REAL(DP), INTENT(OUT) :: v3x 912 !! derivatives of exchange wr. tau 913 REAL(DP), INTENT(OUT) :: v1c 914 !! derivatives of correlation wr. rho 915 REAL(DP), INTENT(OUT) :: v2c 916 !! derivatives of correlation wr. grho 917 REAL(DP), INTENT(OUT) :: v3c 918 !! derivatives of correlation wr. tau 919 ! 920 ! ... local variables 921 ! 922 REAL(DP) :: rhoa, rhob, grho2a, grho2b, taua, taub, v1cb, v2cb, v3cb 923 REAL(DP), PARAMETER :: zero = 0.0_dp, two = 2.0_dp, four = 4.0_dp 924 ! 925 ! 926 rhoa = rho / two ! one component only 927 rhob = rhoa 928 ! 929 grho2a = grho2 / four 930 grho2b = grho2a 931 ! 932 taua = tau * two * 0.5_dp ! Taua, which is Tau_sigma is half Tau 933 taub = taua ! Tau is defined as summ_i( |nabla phi_i|**2 ) 934 ! in the M06L routine 935 ! 936 CALL m06lx( rhoa, grho2a, taua, ex, v1x, v2x, v3x ) !<GPU:m06lx=>m06lx_d> 937 ! 938 ex = two * ex ! Add the two components up + dw 939 ! 940 v2x = 0.5_dp * v2x 941 ! 942 CALL m06lc( rhoa, rhob, grho2a, grho2b, taua, taub, ec, v1c, v2c, v3c, & !<GPU:m06lc=>m06lc_d> 943 v1cb, v2cb, v3cb ) 944 ! 945 v2c = 0.5_dp * v2c 946 ! 947END SUBROUTINE m06lxc 948! 949! 950!------------------------------------------------------------------------- 951SUBROUTINE m06lxc_spin( rhoup, rhodw, grhoup2, grhodw2, tauup, taudw, & !<GPU:DEVICE> 952 ex, ec, v1xup, v1xdw, v2xup, v2xdw, v3xup, v3xdw, & 953 v1cup, v1cdw, v2cup, v2cdw, v3cup, v3cdw ) 954 !----------------------------------------------------------------------- 955 ! 956 USE kinds, ONLY : DP 957 ! 958 IMPLICIT NONE 959 ! 960 REAL(DP), INTENT(IN) :: rhoup, rhodw, grhoup2, grhodw2, tauup, taudw 961 REAL(DP), INTENT(OUT) :: ex, ec, v1xup, v1xdw, v2xup, v2xdw, v3xup, v3xdw, & 962 v1cup, v1cdw, v2cup, v2cdw, v3cup, v3cdw 963 ! 964 REAL(DP) :: exup, exdw, taua, taub 965 REAL(DP), PARAMETER :: zero = 0.0_dp, two = 2.0_dp 966 ! 967 taua = tauup * two ! Tau is defined as summ_i( |nabla phi_i|**2 ) 968 taub = taudw * two ! in the rest of the routine 969 ! 970 CALL m06lx( rhoup, grhoup2, taua, exup, v1xup, v2xup, v3xup ) !<GPU:m06lx=>m06lx_d> 971 CALL m06lx( rhodw, grhodw2, taub, exdw, v1xdw, v2xdw, v3xdw ) !<GPU:m06lx=>m06lx_d> 972 ! 973 ex = exup + exdw 974 ! 975 CALL m06lc( rhoup, rhodw, grhoup2, grhodw2, taua, taub, & !<GPU:m06lc=>m06lc_d> 976 ec, v1cup, v2cup, v3cup, v1cdw, v2cdw, v3cdw ) 977 ! 978END SUBROUTINE m06lxc_spin 979! ! 980! 981!=============================== M06L exchange ========================== 982! 983!------------------------------------------------------------------------------- 984SUBROUTINE m06lx( rho, grho2, tau, ex, v1x, v2x, v3x ) !<GPU:DEVICE> 985 !--------------------------------------------------------------------------- 986 !! M06L exchange. 987 ! 988 USE kinds, ONLY : DP 989 USE constants, ONLY : pi 990 ! 991 IMPLICIT NONE 992 ! 993 REAL(DP), INTENT(IN) :: rho 994 !! the total charge 995 REAL(DP), INTENT(IN) :: grho2 996 !! square of gradient of rho 997 REAL(DP), INTENT(IN) :: tau 998 !! the total kinetic energy density 999 REAL(DP), INTENT(OUT) :: ex 1000 !! exchange energy 1001 REAL(DP), INTENT(OUT) :: v1x 1002 !! derivatives of exchange wr. rho 1003 REAL(DP), INTENT(OUT) :: v2x 1004 !! derivatives of exchange wr. grho 1005 REAL(DP), INTENT(OUT) :: v3x 1006 !! derivatives of exchange wr. tau 1007 ! 1008 ! ... local variables 1009 ! 1010 REAL(DP) :: v1x_unif, ex_unif, ex_pbe, sx_pbe, v1x_pbe, v2x_pbe 1011 ! ex_unif: lda \epsilon_x(rho) 1012 ! v2x = 1/|grho| * dsx / d|grho| = 2 * dsx / dgrho2 1013 ! 1014 REAL(DP), PARAMETER :: zero = 0._dp, one = 1.0_dp, two = 2.0_dp, three = 3.0_dp, & 1015 four = 4.0_dp, five = 5.0_dp, six = 6.0_dp, eight = 8.0_dp, & 1016 f12 = one/two, f13 = one/three, f23 = two/three, & 1017 f53 = five/three, f83 = eight/three, f43 = four/three, & 1018 pi34 = pi*three/four, pi2 = pi*pi, small = 1.d-10 1019 ! 1020 REAL(DP) :: d0, d1, d2, d3, d4, d5, CF, CT, CX, alpha 1021 REAL(DP), DIMENSION(0:11) :: at 1022 INTEGER :: i 1023 ! 1024 ! VSXC98 variables (LDA part) 1025 ! 1026 REAL(DP) :: xs, xs2, grho, rhom83, rho13, rho43, zs, gh 1027 REAL(DP) :: hg, dhg_dxs2, dhg_dzs 1028 REAL(DP) :: dxs2_drho, dxs2_dgrho2, dzs_drho, dzs_dtau 1029 REAL(DP) :: ex_vs98, v1x_vs98, v2x_vs98, v3x_vs98, v2x_vs98_g 1030 ! 1031 ! GGA and MGGA variables 1032 ! 1033 REAL(DP) :: tau_unif, ts, ws, fws, dfws, dfws_drho, dfws_dtau, & 1034 dws_dts, dts_dtau, dts_drho 1035 ! 1036 ! set parameters 1037 ! 1038 at(0) = 3.987756d-01 1039 at(1) = 2.548219d-01 1040 at(2) = 3.923994d-01 1041 at(3) = -2.103655d+00 1042 at(4) = -6.302147d+00 1043 at(5) = 1.097615d+01 1044 at(6) = 3.097273d+01 1045 at(7) = -2.318489d+01 1046 at(8) = -5.673480d+01 1047 at(9) = 2.160364d+01 1048 at(10) = 3.421814d+01 1049 at(11) = -9.049762d+00 1050 ! 1051 d0 = 6.012244d-01 1052 d1 = 4.748822d-03 1053 d2 = -8.635108d-03 1054 d3 = -9.308062d-06 1055 d4 = 4.482811d-05 1056 d5 = zero 1057 ! 1058 alpha = 1.86726d-03 1059 ! 1060 IF (rho < small .AND. tau < small) THEN 1061 ex = zero 1062 v1x = zero 1063 v2x = zero 1064 v3x = zero 1065 RETURN 1066 ENDIF 1067 ! 1068 ! ... VSXC98 functional (LDA part) 1069 ! 1070 ! set variables 1071 ! 1072 CF = (three/five) * (six*pi2)**f23 1073 CT = CF / two 1074 CX = -(three/two) * (three/(four*pi))**f13 ! Cx LSDA 1075 ! 1076 ! IF (rho >= small .AND. grho>=small) THEN 1077 ! 1078 grho = SQRT(grho2) 1079 rho43 = rho**f43 1080 rho13 = rho**f13 1081 rhom83 = one / rho**f83 1082 ! 1083 xs = grho / rho43 1084 xs2 = xs * xs 1085 zs = tau / rho**f53 - CF 1086 gh = one + alpha * (xs2 + zs) 1087 ! 1088 IF (gh >= small) THEN 1089 CALL gvt4( xs2, zs, d0, d1, d2, d3, d4, d5, alpha, hg, dhg_dxs2, dhg_dzs ) !<GPU:gvt4=>gvt4_d> 1090 ELSE 1091 hg = zero 1092 dhg_dxs2 = zero 1093 dhg_dzs = zero 1094 ENDIF 1095 ! 1096 dxs2_drho = -f83*xs2/rho 1097 dxs2_dgrho2 = rhom83 1098 dzs_drho = -f53*tau*rhom83 1099 dzs_dtau = one/rho**f53 1100 ! 1101 ex_unif = CX * rho43 1102 ex_vs98 = ex_unif * hg 1103 v1x_vs98 = CX * ( f43 * hg * rho**f13 ) + & 1104 ex_unif * ( dhg_dxs2*dxs2_drho + dhg_dzs*dzs_drho ) 1105 v2x_vs98 = two * ex_unif * dhg_dxs2 * dxs2_dgrho2 1106 v3x_vs98 = ex_unif * dhg_dzs * dzs_dtau 1107 ! 1108 ! ... mo6lx functional 1109 ! 1110 tau_unif = CF * rho**f53 ! Tau is define as summ_i( |nabla phi_i|**2 ) 1111 ts = tau_unif / tau 1112 ws = (ts - one)/(ts + one) 1113 ! 1114 fws = zero 1115 dfws = zero 1116 ! 1117 DO i = 0, 11 1118 fws = fws + at(i)*ws**i 1119 dfws = dfws + i*at(i)*ws**(i-1) 1120 ENDDO 1121 ! 1122 dws_dts = two/((ts+1)**2) 1123 dts_drho = ( (six*pi*pi*rho)**f23 )/tau 1124 dts_dtau = -ts/tau 1125 dfws_drho = dfws*dws_dts*dts_drho 1126 dfws_dtau = dfws*dws_dts*dts_dtau 1127 ! 1128 CALL pbex_m06l( two*rho, four*grho2, sx_pbe, v1x_pbe, v2x_pbe ) !<GPU:pbex_m06l=>pbex_m06l_d> 1129 ! 1130 v1x_unif = f43 * CX * rho13 1131 ! 1132 sx_pbe = f12 * sx_pbe 1133 v1x_pbe = v1x_pbe + v1x_unif 1134 v2x_pbe = two * v2x_pbe 1135 ! 1136 ex_pbe = sx_pbe + ex_unif 1137 ! 1138 ! ... energy and potential 1139 ! 1140 ex = ex_vs98 + ex_pbe*fws 1141 ! 1142 v1x = v1x_vs98 + v1x_pbe*fws + ex_pbe*dfws_drho 1143 v2x = v2x_vs98 + v2x_pbe*fws 1144 v3X = v3x_vs98 + ex_pbe*dfws_dtau 1145 ! 1146END SUBROUTINE m06lx 1147! 1148! 1149!------------------------------------------------------------------- 1150SUBROUTINE pbex_m06l( rho, grho2, sx, v1x, v2x ) !<GPU:DEVICE> 1151 !--------------------------------------------------------------- 1152 !! PBE exchange (without Slater exchange): 1153 !! J.P.Perdew, K.Burke, M.Ernzerhof, PRL 77, 3865 (1996). 1154 ! 1155 !! v2x = 1/|grho| * dsx / d|grho| = 2 * dsx / dgrho2 1156 ! 1157 USE kinds 1158 USE constants, ONLY : pi 1159 ! 1160 IMPLICIT NONE 1161 ! 1162 REAL(DP), INTENT(IN) :: rho 1163 !! charge density 1164 REAL(DP), INTENT(IN) :: grho2 1165 !! squared gradient 1166 REAL(DP), INTENT(OUT) :: sx 1167 !! energy 1168 REAL(DP), INTENT(OUT) :: v1x 1169 !! potential w.r. rho 1170 REAL(DP), INTENT(OUT) :: v2x 1171 !! potential w.r. grho 1172 ! 1173 ! ... local variables 1174 ! 1175 INTEGER :: iflag 1176 REAL(DP) :: grho, rho43, xs, xs2, dxs2_drho, dxs2_dgrho2 1177 REAL(DP) :: CX, denom, C1, C2, ex, Fx, dFx_dxs2, dex_drho 1178 ! 1179 REAL(DP), PARAMETER :: mu = 0.21951_dp, ka = 0.804_dp, & 1180 one = 1.0_dp, two = 2.0_dp, three = 3.0_dp, & 1181 four = 4.0_dp, six = 6.0_dp, eight = 8.0_dp, & 1182 f13 = one/three, f23 = two/three, f43 = four/three, & 1183 f34 = three/four, f83 = eight/three 1184 ! 1185 CX = f34 * (three/pi)**f13 ! Cx LDA 1186 denom = four * (three*pi**two)**f23 1187 C1 = mu / denom 1188 C2 = mu / (ka * denom) 1189 ! 1190 grho = SQRT(grho2) 1191 rho43 = rho**f43 1192 xs = grho / rho43 1193 xs2 = xs * xs 1194 ! 1195 dxs2_drho = -f83 * xs2 / rho 1196 dxs2_dgrho2 = one /rho**f83 1197 ! 1198 ex = - CX * rho43 1199 dex_drho = - f43 * CX * rho**f13 1200 ! 1201 Fx = C1*xs2 / (one + C2*xs2) 1202 dFx_dxs2 = C1 / (one + C2*xs2)**2 1203 ! 1204 ! Energy 1205 ! 1206 sx = Fx * ex 1207 ! 1208 ! Potential 1209 ! 1210 v1x = dex_drho * Fx + ex * dFx_dxs2 * dxs2_drho 1211 v2x = two * ex * dFx_dxs2* dxs2_dgrho2 1212 ! 1213END SUBROUTINE pbex_m06l 1214! 1215! 1216!=============================== M06L correlation ========================== 1217! 1218!--------------------------------------------------------------------------------- 1219SUBROUTINE m06lc( rhoa, rhob, grho2a, grho2b, taua, taub, ec, v1c_up, v2c_up, & !<GPU:DEVICE> 1220 v3c_up, v1c_dw, v2c_dw, v3c_dw ) 1221 !------------------------------------------------------------------------------- 1222 !! M06L correlation. 1223 ! 1224 USE kinds, ONLY : DP 1225 USE constants, ONLY : pi 1226 ! 1227 IMPLICIT NONE 1228 ! 1229 REAL(DP), INTENT(IN) :: rhoa 1230 !! charge density up 1231 REAL(DP), INTENT(IN) :: rhob 1232 !! charge density down 1233 REAL(DP), INTENT(IN) :: grho2a 1234 !! squared gradient up 1235 REAL(DP), INTENT(IN) :: grho2b 1236 !! squared gradient down 1237 REAL(DP), INTENT(IN) :: taua 1238 !! kinetic energy density up 1239 REAL(DP), INTENT(IN) :: taub 1240 !! kinetic energyt density down 1241 REAL(DP), INTENT(OUT) :: ec 1242 !! correlation energy 1243 REAL(DP), INTENT(OUT) :: v1c_up 1244 !! corr. potential w.r. rho - up 1245 REAL(DP), INTENT(OUT) :: v2c_up 1246 !! corr. potential w.r. rho - up 1247 REAL(DP), INTENT(OUT) :: v3c_up 1248 !! corr. potential w.r. rho - up 1249 REAL(DP), INTENT(OUT) :: v1c_dw 1250 !! corr. potential w.r. rho - down 1251 REAL(DP), INTENT(OUT) :: v2c_dw 1252 !! corr. potential w.r. grho - down 1253 REAL(DP), INTENT(OUT) :: v3c_dw 1254 !! corr. potential w.r. tau - down 1255 ! 1256 ! ... local variables 1257 ! 1258 REAL(DP) :: rs, zeta 1259 REAL(DP) :: vc_v(2) 1260 ! 1261 REAL(DP), PARAMETER :: zero = 0._dp, one = 1.0_dp, two = 2.0_dp, three = 3.0_dp, & 1262 four = 4.0_dp, five = 5.0_dp, six = 6.0_dp, eight = 8.0_dp, & 1263 f12 = one/two, f13 = one/three, f23 = two/three, & 1264 f53 = five/three, f83 = eight/three, f43 = four/three, & 1265 pi34 = three/(four*pi), pi2 = pi*pi, f35 = three/five, & 1266 small = 1.d-10 1267 ! 1268 ! parameters of the MO6Lc functional 1269 ! 1270 REAL(DP), DIMENSION(0:4):: cs, cab 1271 ! 1272 REAL(DP) :: ds0, ds1, ds2, ds3, ds4, ds5, CF, alpha, Ds, & 1273 dab0, dab1, dab2, dab3, dab4, dab5, gama_ab, gama_s, & 1274 alpha_s, alpha_ab 1275 ! 1276 ! functions and variables 1277 ! 1278 REAL(DP) :: ec_pw_a, ec_pw_b, ec_pw_ab 1279 ! 1280 REAL(DP) :: vv, vc_pw_a, vc_pw_b, vc_pw_ab, vc_pw_up, vc_pw_dw, Ecaa, Ecbb, Ecab, & 1281 Ec_UEG_ab, Ec_UEG_aa, Ec_UEG_bb, decab_drhoa, decab_drhob, & 1282 v1_ab_up, v1_ab_dw, v2_ab_up, v2_ab_dw, v3_ab_up, v3_ab_dw, & 1283 v1_aa_up, v2_aa_up, v3_aa_up, v1_bb_dw, v2_bb_dw, v3_bb_dw 1284 ! 1285 REAL(DP) :: xsa, xs2a, rsa, grhoa, xsb, xs2b, grhob, rsb, zsa, zsb, & 1286 xs2ab, zsab, rho, & 1287 dxs2a_drhoa, dxs2b_drhob, dxs2a_dgrhoa2, dxs2b_dgrhob2, & 1288 dzsa_drhoa, dzsb_drhob, dzsa_dtaua, dzsb_dtaub 1289 ! 1290 REAL(DP) :: hga, dhga_dxs2a, dhga_dzsa, hgb, dhgb_dxs2b, dhgb_dzsb, & 1291 hgab, dhgab_dxs2ab, dhgab_dzsab, & 1292 Dsa, Dsb, dDsa_dxs2a, dDsa_dzsa, dDsb_dxs2b, dDsb_dzsb, & 1293 gsa, gsb, gsab, dgsa_dxs2a, dgsb_dxs2b, dgsab_dxs2ab, num 1294 ! 1295 INTEGER :: ifunc 1296 ! 1297 ! 1298 dab0 = 3.957626d-01 1299 dab1 = -5.614546d-01 1300 dab2 = 1.403963d-02 1301 dab3 = 9.831442d-04 1302 dab4 = -3.577176d-03 1303 dab5 = zero 1304 ! 1305 cab(0) = 6.042374d-01 1306 cab(1) = 1.776783d+02 1307 cab(2) = -2.513252d+02 1308 cab(3) = 7.635173d+01 1309 cab(4) = -1.255699d+01 1310 ! 1311 gama_ab = 0.0031_dp 1312 alpha_ab = 0.00304966_dp 1313 ! 1314 ds0 = 4.650534d-01 1315 ds1 = 1.617589d-01 1316 ds2 = 1.833657d-01 1317 ds3 = 4.692100d-04 1318 ds4 = -4.990573d-03 1319 ds5 = zero 1320 ! 1321 cs(0) = 5.349466d-01 1322 cs(1) = 5.396620d-01 1323 cs(2) = -3.161217d+01 1324 cs(3) = 5.149592d+01 1325 cs(4) = -2.919613d+01 1326 ! 1327 gama_s = 0.06_dp 1328 alpha_s = 0.00515088_dp 1329 ! 1330 CF = f35 * (six*pi2)**f23 1331 ! 1332 ifunc = 1 ! iflag=1 J.P. Perdew and Y. Wang, PRB 45, 13244 (1992) 1333 ! 1334 ! ... Ecaa 1335 ! 1336 IF (rhoa < small .AND. taua < small) THEN 1337 ! 1338 Ecaa = zero 1339 v1_aa_up = zero 1340 v2_aa_up = zero 1341 v3_aa_up = zero 1342 ! 1343 ELSE 1344 ! 1345 rsa = (pi34/rhoa)**f13 1346 grhoa = SQRT(grho2a) 1347 xsa = grhoa / rhoa**f43 1348 xs2a = xsa * xsa 1349 zsa = taua/rhoa**f53 - CF 1350 ! 1351 dxs2a_drhoa = -f83*xs2a/rhoa 1352 dxs2a_dgrhoa2 = one/(rhoa**f83) 1353 ! 1354 dzsa_drhoa = -f53*taua/(rhoa**f83) 1355 dzsa_dtaua = one/rhoa**f53 1356 ! 1357 Dsa = one - xs2a/(four * (zsa + CF)) 1358 dDsa_dxs2a = - one/(four * (zsa + CF)) 1359 dDsa_dzsa = xs2a/(four * (zsa + CF)**2) 1360 ! 1361 ec_pw_a = zero 1362 vc_pw_a = zero 1363 ! 1364 rs = rsa 1365 zeta = one 1366 CALL pw_spin( rs, zeta, ec_pw_a, vc_v(1), vc_v(2) ) !<GPU:pw_spin=>pw_spin_d> 1367 vc_pw_a = vc_v(1) 1368 vv = vc_v(2) 1369 ! 1370 CALL gvt4( xs2a, zsa, ds0, ds1, ds2, ds3, ds4, ds5, alpha_s, hga, dhga_dxs2a, dhga_dzsa ) !<GPU:gvt4=>gvt4_d> 1371 CALL gfunc( cs, gama_s, xs2a, gsa, dgsa_dxs2a ) !<GPU:gfunc=>gfunc_d> 1372 ! 1373 Ec_UEG_aa = rhoa*ec_pw_a 1374 num = (dgsa_dxs2a + dhga_dxs2a)*Dsa + (gsa + hga)*dDsa_dxs2a 1375 ! 1376 ! 1377 Ecaa = Ec_UEG_aa * (gsa + hga) * Dsa 1378 ! 1379 v1_aa_up = vc_pw_a * (gsa + hga) * Dsa + & 1380 Ec_UEG_aa * num * dxs2a_drhoa + & 1381 Ec_UEG_aa * (dhga_dzsa*Dsa + (gsa + hga)*dDsa_dzsa) * dzsa_drhoa 1382 ! 1383 v2_aa_up = two * Ec_UEG_aa * num * dxs2a_dgrhoa2 1384 ! 1385 v3_aa_up = Ec_UEG_aa * (dhga_dzsa*Dsa + (gsa + hga)*dDsa_dzsa) * dzsa_dtaua 1386 ! 1387 ENDIF 1388 ! 1389 ! ... Ecbb 1390 ! 1391 IF (rhob < small .AND. taub < small) THEN 1392 ! 1393 Ecbb = zero 1394 v1_bb_dw = zero 1395 v2_bb_dw = zero 1396 v3_bb_dw = zero 1397 ! 1398 ELSE 1399 ! 1400 rsb = (pi34/rhob)**f13 1401 grhob = SQRT(grho2b) 1402 xsb = grhob / rhob**f43 1403 xs2b = xsb * xsb 1404 zsb = taub/rhob**f53 - CF 1405 1406 dxs2b_drhob = -f83*xs2b/rhob 1407 dxs2b_dgrhob2 = one /rhob**f83 1408 1409 dzsb_drhob = -f53*taub/(rhob**f83) 1410 dzsb_dtaub = one/rhob**f53 1411 1412 Dsb = one - xs2b/(four * (zsb + CF)) 1413 dDsb_dxs2b = - one/(four * (zsb + CF)) 1414 dDsb_dzsb = xs2b/(four * (zsb + CF)**2) 1415 ! 1416 zeta = one 1417 rs = rsb 1418 CALL pw_spin( rs, zeta, ec_pw_b, vc_v(1), vc_v(2) ) !<GPU:pw_spin=>pw_spin_d> 1419 vc_pw_b = vc_v(1) 1420 vv = vc_v(2) 1421 ! 1422 CALL gvt4( xs2b, zsb, ds0, ds1, ds2, ds3, ds4, ds5, alpha_s, hgb, dhgb_dxs2b, dhgb_dzsb ) !<GPU:gvt4=>gvt4_d> 1423 CALL gfunc( cs, gama_s, xs2b, gsb, dgsb_dxs2b ) !<GPU:gfunc=>gfunc_d> 1424 ! 1425 Ec_UEG_bb = rhob*ec_pw_b 1426 num = (dgsb_dxs2b + dhgb_dxs2b)*Dsb + (gsb + hgb)*dDsb_dxs2b 1427 ! 1428 Ecbb = Ec_UEG_bb * (gsb + hgb) * Dsb 1429 ! 1430 v1_bb_dw = vc_pw_b * (gsb + hgb) * Dsb + & 1431 Ec_UEG_bb * num * dxs2b_drhob + & 1432 Ec_UEG_bb * (dhgb_dzsb*Dsb + (gsb + hgb)*dDsb_dzsb)*dzsb_drhob 1433 ! 1434 v2_bb_dw = two * Ec_UEG_bb * num * dxs2b_dgrhob2 1435 ! 1436 v3_bb_dw = Ec_UEG_bb * (dhgb_dzsb*Dsb + (gsb + hgb)*dDsb_dzsb)*dzsb_dtaub 1437 ! 1438 ENDIF 1439 ! 1440 ! ... Ecab 1441 ! 1442 IF (rhoa < small .AND. rhob < small) THEN 1443 ! 1444 Ecab = zero 1445 v1_ab_up = zero 1446 v1_ab_dw = zero 1447 v2_ab_up = zero 1448 v2_ab_dw = zero 1449 v3_ab_up = zero 1450 v3_ab_dw = zero 1451 ! 1452 ELSE 1453 ! 1454 xs2ab = xs2a + xs2b 1455 zsab = zsa + zsb 1456 rho = rhoa + rhob 1457 zeta = (rhoa - rhob)/rho 1458 rs = (pi34/rho)**f13 1459 ! 1460 CALL gvt4( xs2ab, zsab, dab0, dab1, dab2, dab3, dab4, dab5, alpha_ab, hgab, & !<GPU:gvt4=>gvt4_d> 1461 dhgab_dxs2ab, dhgab_dzsab ) 1462 ! 1463 CALL pw_spin( rs, zeta, ec_pw_ab, vc_v(1), vc_v(2) ) !<GPU:pw_spin=>pw_spin_d> 1464 vc_pw_up=vc_v(1) ; vc_pw_dw=vc_v(2) 1465 ! 1466 CALL gfunc( cab, gama_ab, xs2ab, gsab, dgsab_dxs2ab ) !<GPU:gfunc=>gfunc_d> 1467 ! 1468 decab_drhoa = vc_pw_up - vc_pw_a 1469 decab_drhob = vc_pw_dw - vc_pw_b 1470 ! 1471 Ec_UEG_ab = ec_pw_ab*rho - ec_pw_a*rhoa - ec_pw_b*rhob 1472 ! 1473 Ecab = Ec_UEG_ab * (gsab + hgab) 1474 ! 1475 v1_ab_up = decab_drhoa * (gsab + hgab) + & 1476 Ec_UEG_ab * (dgsab_dxs2ab + dhgab_dxs2ab) * dxs2a_drhoa + & 1477 Ec_UEG_ab * dhgab_dzsab * dzsa_drhoa 1478 ! 1479 v1_ab_dw = decab_drhob * (gsab + hgab) + & 1480 Ec_UEG_ab * (dgsab_dxs2ab + dhgab_dxs2ab) * dxs2b_drhob + & 1481 Ec_UEG_ab * dhgab_dzsab * dzsb_drhob 1482 ! 1483 v2_ab_up = two * Ec_UEG_ab * (dgsab_dxs2ab + dhgab_dxs2ab) * dxs2a_dgrhoa2 1484 v2_ab_dw = two * Ec_UEG_ab * (dgsab_dxs2ab + dhgab_dxs2ab) * dxs2b_dgrhob2 1485 1486 v3_ab_up = Ec_UEG_ab * dhgab_dzsab * dzsa_dtaua 1487 v3_ab_dw = Ec_UEG_ab * dhgab_dzsab * dzsb_dtaub 1488 ! 1489 ENDIF 1490 ! 1491 ! ... ec and vc 1492 ! 1493 ec = Ecaa + Ecbb + Ecab 1494 ! 1495 v1c_up = v1_aa_up + v1_ab_up 1496 v2c_up = v2_aa_up + v2_ab_up 1497 v3c_up = v3_aa_up + v3_ab_up 1498 ! 1499 v1c_dw = v1_bb_dw + v1_ab_dw 1500 v2c_dw = v2_bb_dw + v2_ab_dw 1501 v3c_dw = v3_bb_dw + v3_ab_dw 1502 ! 1503END SUBROUTINE m06lc 1504 ! 1505 ! 1506 !------------------------------------------------------------------- 1507SUBROUTINE gfunc( cspin, gama, xspin, gs, dgs_dx ) !<GPU:DEVICE> 1508 !----------------------------------------------------------------- 1509 ! 1510 USE kinds, ONLY : DP 1511 ! 1512 IMPLICIT NONE 1513 ! 1514 REAL(DP), INTENT(IN) :: cspin(0:4) 1515 REAL(DP), INTENT(IN) :: xspin, gama 1516 REAL(DP), INTENT(OUT) :: gs, dgs_dx 1517 ! 1518 REAL(DP) :: de, d2, x1, x2, x3, x4 1519 REAL(DP), PARAMETER :: one=1.0d0, two=2.0d0, & 1520 three=3.0d0, four=4.0d0 1521 ! 1522 de = one/(one + gama*xspin) 1523 d2 = de**2 1524 x1 = gama * xspin * de 1525 x2 = x1**2 1526 x3 = x1**3 1527 x4 = x1**4 1528 ! 1529 gs = cspin(0) + cspin(1)*x1 + cspin(2)*x2 + cspin(3)*x3 + cspin(4)*x4 1530 dgs_dx = gama*d2*(cspin(1) + two*cspin(2)*x1 + three*cspin(3)*x2 + four*cspin(4)*x3) 1531 ! 1532 END SUBROUTINE gfunc 1533 ! 1534! 1535! 1536!------------------------------------------------------------------------- 1537SUBROUTINE gvt4( x, z, a, b, c, d, e, f, alpha, hg, dh_dx, dh_dz ) !<GPU:DEVICE> 1538 !---------------------------------------------------------------------- 1539 ! 1540 USE kinds, ONLY : DP 1541 ! 1542 IMPLICIT NONE 1543 ! 1544 REAL(DP), INTENT(IN) :: X, z, a, b, c, d, e, f, alpha 1545 REAL(DP), INTENT(OUT) :: hg, dh_dx, dh_dz 1546 ! 1547 REAL(DP) :: gamma, gamma2, gamma3 1548 REAL(DP), PARAMETER :: one=1.0_dp, two=2.0_dp, three=3.0_dp 1549 ! 1550 gamma = one + alpha*(x+z) 1551 gamma2 = gamma*gamma 1552 gamma3 = gamma2*gamma 1553 ! 1554 hg = a/gamma + (b*x + c*z)/gamma2 + (d*x*x + e*x*z + f*z*z)/gamma3 1555 ! 1556 dh_dx = ( -alpha*a + b + (two*x*(d - alpha*b) + z*(e - two*alpha*c))/ gamma & 1557 - three*alpha*(d*x*x + e*x*z + f*z*z)/gamma2 )/gamma2 1558 ! 1559 dh_dz = ( -alpha*a + c + (two*z*(f - alpha*c) + x*(e -two*alpha*b))/ gamma & 1560 - three*alpha*(d*x*x + e*x*z + f*z*z)/gamma2 )/gamma2 1561 ! 1562 RETURN 1563 ! 1564END SUBROUTINE gvt4 1565! 1566! END M06L 1567!========================================================================= 1568END MODULE 1569