1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Calculates the exchange energy for the pbe hole model in a truncated 8!> coulomb potential, considering the long range part only. Can be used 9!> as longrange correction to a truncated Hartree Fock calculation 10!> \par History 11!> Manuel Guidon (01.2009) : initial version 12!> \author Manuel Guidon (01.2009) 13! ************************************************************************************************** 14 15MODULE xc_xpbe_hole_t_c_lr 16 USE input_section_types, ONLY: section_vals_type,& 17 section_vals_val_get 18 USE kinds, ONLY: dp 19 USE mathconstants, ONLY: euler,& 20 pi 21 USE mathlib, ONLY: expint 22 USE xc_derivative_set_types, ONLY: xc_derivative_set_type,& 23 xc_dset_get_derivative 24 USE xc_derivative_types, ONLY: xc_derivative_get,& 25 xc_derivative_type 26 USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type 27 USE xc_rho_set_types, ONLY: xc_rho_set_get,& 28 xc_rho_set_type 29#include "../base/base_uses.f90" 30 31 IMPLICIT NONE 32 33 PRIVATE 34 35! *** Global parameters *** 36 37 PUBLIC :: xpbe_hole_t_c_lr_lda_eval, xpbe_hole_t_c_lr_lda_info, & 38 xpbe_hole_t_c_lr_lda_calc_1, xpbe_hole_t_c_lr_lda_calc_2, & 39 xpbe_hole_t_c_lr_lsd_eval, xpbe_hole_t_c_lr_lsd_info 40 41 REAL(KIND=dp), PARAMETER :: a1 = 0.00979681_dp, & 42 a2 = 0.04108340_dp, & 43 a3 = 0.18744000_dp, & 44 a4 = 0.00120824_dp, & 45 a5 = 0.0347188_dp 46 REAL(KIND=dp), PARAMETER :: A = 1.0161144_dp, & 47 B = -0.37170836_dp, & 48 C = -0.077215461_dp, & 49 D = 0.57786348_dp, & 50 E = -0.051955731_dp, & 51 F2 = 0.47965830_dp, & 52 F1 = 6.4753871_dp, & 53 clda = -0.73855876638202240588423_dp 54 REAL(KIND=dp), PARAMETER :: expcutoff = 700.0_dp 55 REAL(KIND=dp), PARAMETER :: smax = 8.572844_dp, & 56 sconst = 18.79622316_dp, & 57 scutoff = 8.3_dp 58 REAL(KIND=dp), PARAMETER :: gcutoff = 0.08_dp, & 59 g1 = -0.02628417880_dp/E, & 60 g2 = -0.07117647788_dp/E, & 61 g3 = 0.08534541323_dp/E, & 62 g4 = 0.0_dp 63 REAL(KIND=dp), PARAMETER :: wcutoff = 14.0_dp 64 65 REAL(KIND=dp), PARAMETER :: EPS_RHO = 0.00000001_dp 66 67 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_xpbe_hole_t_c_lr' 68 69CONTAINS 70 71! ************************************************************************************************** 72!> \brief returns various information on the functional 73!> \param reference string with the reference of the actual functional 74!> \param shortform string with the shortform of the functional name 75!> \param needs the components needed by this functional are set to 76!> true (does not set the unneeded components to false) 77!> \param max_deriv controls the number of derivatives 78!> \par History 79!> 01.2009 created [mguidon] 80!> \author mguidon 81! ************************************************************************************************** 82 SUBROUTINE xpbe_hole_t_c_lr_lda_info(reference, shortform, needs, max_deriv) 83 CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform 84 TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs 85 INTEGER, INTENT(out), OPTIONAL :: max_deriv 86 87 IF (PRESENT(reference)) THEN 88 reference = "{LDA version}" 89 END IF 90 IF (PRESENT(shortform)) THEN 91 shortform = "{LDA}" 92 END IF 93 IF (PRESENT(needs)) THEN 94 needs%rho = .TRUE. 95 needs%norm_drho = .TRUE. 96 END IF 97 IF (PRESENT(max_deriv)) max_deriv = 1 98 99 END SUBROUTINE xpbe_hole_t_c_lr_lda_info 100 101! ************************************************************************************************** 102!> \brief returns various information on the functional 103!> \param reference string with the reference of the actual functional 104!> \param shortform string with the shortform of the functional name 105!> \param needs the components needed by this functional are set to 106!> true (does not set the unneeded components to false) 107!> \param max_deriv controls the number of derivatives 108!> \par History 109!> 01.2009 created [mguidon] 110!> \author mguidon 111! ************************************************************************************************** 112 SUBROUTINE xpbe_hole_t_c_lr_lsd_info(reference, shortform, needs, max_deriv) 113 CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform 114 TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs 115 INTEGER, INTENT(out), OPTIONAL :: max_deriv 116 117 IF (PRESENT(reference)) THEN 118 reference = "{LSD version}" 119 END IF 120 IF (PRESENT(shortform)) THEN 121 shortform = "{LSD}" 122 END IF 123 IF (PRESENT(needs)) THEN 124 needs%rho_spin = .TRUE. 125 needs%norm_drho_spin = .TRUE. 126 END IF 127 IF (PRESENT(max_deriv)) max_deriv = 1 128 129 END SUBROUTINE xpbe_hole_t_c_lr_lsd_info 130 131! ************************************************************************************************** 132!> \brief evaluates the pbe-hole exchange in a truncated coulomb potential 133!> \param rho_set the density where you want to evaluate the functional 134!> \param deriv_set place where to store the functional derivatives (they are 135!> added to the derivatives) 136!> \param order degree of the derivative that should be evaluated, 137!> if positive all the derivatives up to the given degree are evaluated, 138!> if negative only the given degree is calculated 139!> \param params input parameters (scaling,cutoff) 140!> \par History 141!> 01.2009 created [Manuel Guidon] 142!> \author Manuel Guidon 143!> \note 144! ************************************************************************************************** 145 SUBROUTINE xpbe_hole_t_c_lr_lda_eval(rho_set, deriv_set, order, params) 146 147 TYPE(xc_rho_set_type), POINTER :: rho_set 148 TYPE(xc_derivative_set_type), POINTER :: deriv_set 149 INTEGER, INTENT(IN) :: order 150 TYPE(section_vals_type), POINTER :: params 151 152 CHARACTER(len=*), PARAMETER :: routineN = 'xpbe_hole_t_c_lr_lda_eval', & 153 routineP = moduleN//':'//routineN 154 155 INTEGER :: handle, npoints 156 INTEGER, DIMENSION(:, :), POINTER :: bo 157 REAL(kind=dp) :: epsilon_rho, R, sx 158 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: dummy, e_0, e_ndrho, e_rho, norm_drho, & 159 rho 160 TYPE(xc_derivative_type), POINTER :: deriv 161 162 CALL timeset(routineN, handle) 163 164 NULLIFY (bo) 165 166 CALL section_vals_val_get(params, "SCALE_X", r_val=sx) 167 CALL section_vals_val_get(params, "CUTOFF_RADIUS", r_val=R) 168 169 CPASSERT(ASSOCIATED(rho_set)) 170 CPASSERT(rho_set%ref_count > 0) 171 CPASSERT(ASSOCIATED(deriv_set)) 172 CPASSERT(deriv_set%ref_count > 0) 173 174 CALL xc_rho_set_get(rho_set, rho=rho, & 175 norm_drho=norm_drho, local_bounds=bo, rho_cutoff=epsilon_rho) 176 npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1) 177 178 dummy => rho 179 180 e_0 => dummy 181 e_rho => dummy 182 e_ndrho => dummy 183 184 IF (order >= 0) THEN 185 deriv => xc_dset_get_derivative(deriv_set, "", & 186 allocate_deriv=.TRUE.) 187 CALL xc_derivative_get(deriv, deriv_data=e_0) 188 END IF 189 IF (order >= 1 .OR. order == -1) THEN 190 deriv => xc_dset_get_derivative(deriv_set, "(rho)", & 191 allocate_deriv=.TRUE.) 192 CALL xc_derivative_get(deriv, deriv_data=e_rho) 193 deriv => xc_dset_get_derivative(deriv_set, "(norm_drho)", & 194 allocate_deriv=.TRUE.) 195 CALL xc_derivative_get(deriv, deriv_data=e_ndrho) 196 END IF 197 IF (order > 1 .OR. order < -1) THEN 198 CPABORT("derivatives bigger than 2 not implemented") 199 END IF 200 201 IF (R == 0.0_dp) THEN 202 CPABORT("Cutoff_Radius 0.0 not implemented") 203 END IF 204 205!$OMP PARALLEL DEFAULT(NONE) & 206!$OMP SHARED(npoints, order, rho, norm_drho, e_0, e_rho) & 207!$OMP SHARED(e_ndrho, epsilon_rho) & 208!$OMP SHARED(sx, r) 209 210 CALL xpbe_hole_t_c_lr_lda_calc(npoints, order, rho=rho, norm_drho=norm_drho, & 211 e_0=e_0, e_rho=e_rho, e_ndrho=e_ndrho, & 212 epsilon_rho=epsilon_rho, & 213 sx=sx, R=R) 214 215!$OMP END PARALLEL 216 217 CALL timestop(handle) 218 219 END SUBROUTINE xpbe_hole_t_c_lr_lda_eval 220 221! ************************************************************************************************** 222!> \brief intermediate level routine in order to decide which branch has to be 223!> taken 224!> \param npoints number of points on the grid 225!> \param order order of the derivatives 226!> \param rho value of density on the grid 227!> \param norm_drho value of gradient on the grid 228!> \param e_0 derivatives of the energy on the grid 229!> \param e_rho derivatives of the energy on the grid 230!> \param e_ndrho derivatives of the energy on the grid 231!> \param epsilon_rho cutoffs 232!> \param sx functional parameters 233!> \param R functional parameters 234!> \par History 235!> 01.2009 created [Manuel Guidon] 236!> \author Manuel Guidon 237!> \note For numerical reasons there are two different branches 238!> 239! ************************************************************************************************** 240 SUBROUTINE xpbe_hole_t_c_lr_lda_calc(npoints, order, rho, norm_drho, e_0, e_rho, e_ndrho, & 241 epsilon_rho, sx, R) 242 243 INTEGER, INTENT(in) :: npoints, order 244 REAL(kind=dp), DIMENSION(1:npoints), INTENT(inout) :: rho, norm_drho, e_0, e_rho, e_ndrho 245 REAL(kind=dp), INTENT(in) :: epsilon_rho, sx, R 246 247 CHARACTER(len=*), PARAMETER :: routineN = 'xpbe_hole_t_c_lr_lda_calc', & 248 routineP = moduleN//':'//routineN 249 250 INTEGER :: ip 251 REAL(dp) :: my_ndrho, my_rho 252 REAL(KIND=dp) :: ss, ss2, sscale, t1, t2, t3, t4, t6, t7, & 253 t8 254 255!$OMP DO 256 257 DO ip = 1, npoints 258 my_rho = MAX(rho(ip), 0.0_dp) 259 IF (my_rho > epsilon_rho) THEN 260 my_ndrho = MAX(norm_drho(ip), EPSILON(0.0_dp)*1.e4_dp) 261 ! *** Do some precalculation in order to catch the correct branch afterwards 262 sscale = 1.0_dp 263 t1 = pi**2 264 t2 = t1*my_rho 265 t3 = t2**(0.1e1_dp/0.3e1_dp) 266 t4 = 0.1e1_dp/t3 267 t6 = my_ndrho*t4 268 t7 = 0.1e1_dp/my_rho 269 t8 = t7*sscale 270 ss = 0.3466806371753173524216762e0_dp*t6*t8 271 IF (ss > scutoff) THEN 272 ss2 = ss*ss 273 sscale = (smax*ss2 - sconst)/(ss2*ss) 274 END IF 275 IF (ss*sscale > gcutoff) THEN 276 CALL xpbe_hole_t_c_lr_lda_calc_1(e_0(ip), e_rho(ip), e_ndrho(ip), & 277 my_rho, & 278 my_ndrho, sscale, sx, R, order) 279 ELSE 280 CALL xpbe_hole_t_c_lr_lda_calc_2(e_0(ip), e_rho(ip), e_ndrho(ip), & 281 my_rho, & 282 my_ndrho, sscale, sx, R, order) 283 END IF 284 END IF 285 END DO 286 287!$OMP END DO 288 289 END SUBROUTINE xpbe_hole_t_c_lr_lda_calc 290 291! ************************************************************************************************** 292!> \brief evaluates the pbe-hole exchange in a truncated coulomb potential 293!> \param rho_set the density where you want to evaluate the functional 294!> \param deriv_set place where to store the functional derivatives (they are 295!> added to the derivatives) 296!> \param order degree of the derivative that should be evaluated, 297!> if positive all the derivatives up to the given degree are evaluated, 298!> if negative only the given degree is calculated 299!> \param params input parameters (scaling,cutoff) 300!> \par History 301!> 01.2009 created [Manuel Guidon] 302!> \author Manuel Guidon 303!> \note 304! ************************************************************************************************** 305 SUBROUTINE xpbe_hole_t_c_lr_lsd_eval(rho_set, deriv_set, order, params) 306 307 TYPE(xc_rho_set_type), POINTER :: rho_set 308 TYPE(xc_derivative_set_type), POINTER :: deriv_set 309 INTEGER, INTENT(IN) :: order 310 TYPE(section_vals_type), POINTER :: params 311 312 CHARACTER(len=*), PARAMETER :: routineN = 'xpbe_hole_t_c_lr_lsd_eval', & 313 routineP = moduleN//':'//routineN 314 315 INTEGER :: handle, npoints 316 INTEGER, DIMENSION(:, :), POINTER :: bo 317 REAL(kind=dp) :: epsilon_rho, R, sx 318 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: dummy, e_0, e_ndrhoa, e_ndrhob, e_rhoa, & 319 e_rhob, norm_drhoa, norm_drhob, rhoa, & 320 rhob 321 TYPE(xc_derivative_type), POINTER :: deriv 322 323 CALL timeset(routineN, handle) 324 325 NULLIFY (bo) 326 327 CALL section_vals_val_get(params, "SCALE_X", r_val=sx) 328 CALL section_vals_val_get(params, "CUTOFF_RADIUS", r_val=R) 329 330 CPASSERT(ASSOCIATED(rho_set)) 331 CPASSERT(rho_set%ref_count > 0) 332 CPASSERT(ASSOCIATED(deriv_set)) 333 CPASSERT(deriv_set%ref_count > 0) 334 335 CALL xc_rho_set_get(rho_set, rhoa=rhoa, rhob=rhob, & 336 norm_drhoa=norm_drhoa, norm_drhob=norm_drhob, local_bounds=bo, rho_cutoff=epsilon_rho) 337 npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1) 338 339 dummy => rhoa 340 341 e_0 => dummy 342 e_rhoa => dummy 343 e_rhob => dummy 344 e_ndrhoa => dummy 345 e_ndrhob => dummy 346 347 IF (order >= 0) THEN 348 deriv => xc_dset_get_derivative(deriv_set, "", & 349 allocate_deriv=.TRUE.) 350 CALL xc_derivative_get(deriv, deriv_data=e_0) 351 END IF 352 IF (order >= 1 .OR. order == -1) THEN 353 deriv => xc_dset_get_derivative(deriv_set, "(rhoa)", & 354 allocate_deriv=.TRUE.) 355 CALL xc_derivative_get(deriv, deriv_data=e_rhoa) 356 deriv => xc_dset_get_derivative(deriv_set, "(rhob)", & 357 allocate_deriv=.TRUE.) 358 CALL xc_derivative_get(deriv, deriv_data=e_rhob) 359 deriv => xc_dset_get_derivative(deriv_set, "(norm_drhoa)", & 360 allocate_deriv=.TRUE.) 361 CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa) 362 deriv => xc_dset_get_derivative(deriv_set, "(norm_drhob)", & 363 allocate_deriv=.TRUE.) 364 CALL xc_derivative_get(deriv, deriv_data=e_ndrhob) 365 END IF 366 IF (order > 1 .OR. order < -1) THEN 367 CPABORT("derivatives bigger than 2 not implemented") 368 END IF 369 370 IF (R == 0.0_dp) THEN 371 CPABORT("Cutoff_Radius 0.0 not implemented") 372 END IF 373 374!$OMP PARALLEL DEFAULT(NONE) & 375!$OMP SHARED(npoints, order, rhoa, norm_drhoa, e_0, e_rhoa) & 376!$OMP SHARED(epsilon_rho, sx, r) & 377!$OMP SHARED(rhob, norm_drhob, e_rhob, e_ndrhoa, e_ndrhob) 378 379 CALL xpbe_hole_t_c_lr_lsd_calc(npoints, order, rho=rhoa, norm_drho=norm_drhoa, & 380 e_0=e_0, e_rho=e_rhoa, e_ndrho=e_ndrhoa, & 381 epsilon_rho=epsilon_rho, sx=sx, R=R) 382 CALL xpbe_hole_t_c_lr_lsd_calc(npoints, order, rho=rhob, norm_drho=norm_drhob, & 383 e_0=e_0, e_rho=e_rhob, e_ndrho=e_ndrhob, & 384 epsilon_rho=epsilon_rho, sx=sx, R=R) 385 386!$OMP END PARALLEL 387 388 CALL timestop(handle) 389 390 END SUBROUTINE xpbe_hole_t_c_lr_lsd_eval 391 392! ************************************************************************************************** 393!> \brief intermediate level routine in order to decide which branch has to be 394!> taken 395!> \param npoints number of points on the grid 396!> \param order order of the derivatives 397!> \param rho density on the grid 398!> \param norm_drho gradient on the grid 399!> \param e_0 derivatives of the energy on the grid 400!> \param e_rho derivatives of the energy on the grid 401!> \param e_ndrho derivatives of the energy on the grid 402!> \param epsilon_rho cutoffs 403!> \param sx functional parameters 404!> \param R functional parameters 405!> \par History 406!> 01.2009 created [Manuel Guidon] 407!> \author Manuel Guidon 408!> \note For numerical reasons there are two different branches. This code calls 409!> the lda version applying spin scaling relations 410! ************************************************************************************************** 411 SUBROUTINE xpbe_hole_t_c_lr_lsd_calc(npoints, order, rho, norm_drho, e_0, e_rho, e_ndrho, & 412 epsilon_rho, sx, R) 413 414 INTEGER, INTENT(in) :: npoints, order 415 REAL(kind=dp), DIMENSION(1:npoints), INTENT(inout) :: rho, norm_drho, e_0, e_rho, e_ndrho 416 REAL(kind=dp), INTENT(in) :: epsilon_rho, sx, R 417 418 CHARACTER(len=*), PARAMETER :: routineN = 'xpbe_hole_t_c_lr_lsd_calc', & 419 routineP = moduleN//':'//routineN 420 421 INTEGER :: ip 422 REAL(dp) :: my_ndrho, my_rho 423 REAL(KIND=dp) :: e_tmp, ss, ss2, sscale, t1, t2, t3, t4, & 424 t6, t7, t8 425 426!$OMP DO 427 428 DO ip = 1, npoints 429 my_rho = 2.0_dp*MAX(rho(ip), 0.0_dp) 430 IF (my_rho > epsilon_rho) THEN 431 my_ndrho = 2.0_dp*MAX(norm_drho(ip), 0.0_dp) 432 433 ! *** Do some precalculation in order to catch the correct branch afterwards 434 sscale = 1.0_dp 435 t1 = pi**2 436 t2 = t1*my_rho 437 t3 = t2**(0.1e1_dp/0.3e1_dp) 438 t4 = 0.1e1_dp/t3 439 t6 = my_ndrho*t4 440 t7 = 0.1e1_dp/my_rho 441 t8 = t7*sscale 442 ss = 0.3466806371753173524216762e0_dp*t6*t8 443 IF (ss > scutoff) THEN 444 ss2 = ss*ss 445 sscale = (smax*ss2 - sconst)/(ss2*ss) 446 END IF 447 e_tmp = 0.0_dp 448 IF (ss*sscale > gcutoff) THEN 449 CALL xpbe_hole_t_c_lr_lda_calc_1(e_tmp, e_rho(ip), e_ndrho(ip), & 450 my_rho, & 451 my_ndrho, sscale, sx, R, order) 452 ELSE 453 CALL xpbe_hole_t_c_lr_lda_calc_2(e_tmp, e_rho(ip), e_ndrho(ip), & 454 my_rho, & 455 my_ndrho, sscale, sx, R, order) 456 END IF 457 e_0(ip) = e_0(ip) + 0.5_dp*e_tmp 458 459 END IF 460 END DO 461 462!$OMP END DO 463 464 END SUBROUTINE xpbe_hole_t_c_lr_lsd_calc 465 466! ************************************************************************************************** 467!> \brief low level routine that calculates the energy derivatives in one point 468!> \param e_0 derivatives of the energy on the grid 469!> \param e_rho derivatives of the energy on the grid 470!> \param e_ndrho derivatives of the energy on the grid 471!> \param rho value of density on the grid 472!> \param ndrho value of gradient on the grid 473!> \param sscale functional parameters 474!> \param sx functional parameters 475!> \param R functional parameters 476!> \param order order of the derivatives 477!> \par History 478!> 01.2009 created [Manuel Guidon] 479!> \author Manuel Guidon 480! ************************************************************************************************** 481 SUBROUTINE xpbe_hole_t_c_lr_lda_calc_1(e_0, e_rho, e_ndrho, & 482 rho, ndrho, sscale, sx, R, order) 483 REAL(KIND=dp), INTENT(INOUT) :: e_0, e_rho, e_ndrho 484 REAL(KIND=dp), INTENT(IN) :: rho, ndrho, sscale, sx, R 485 INTEGER, INTENT(IN) :: order 486 487 REAL(KIND=dp) :: p, pndrho, prho, q, qndrho, qrho, sndrho, srho, ssval, t1, t103, t104, & 488 t105, t106, t107, t108, t109, t11, t113, t115, t116, t117, t118, t119, t12, t121, t123, & 489 t125, t126, t129, t13, t133, t136, t14, t140, t142, t143, t144, t147, t150, t152, t155, & 490 t156, t159, t162, t163, t164, t167, t17, t170, t173, t175, t176, t177, t178, t181, t183, & 491 t187, t19, t190, t194, t195, t196, t199, t2, t202, t209, t21, t214, t216, t217, t218, & 492 t22, t222, t223, t224, t227, t228, t229, t232, t234, t235, t236, t237, t240, t244, t248, & 493 t25, t251, t255, t256, t257, t261, t262, t265, t269, t27, t273, t276 494 REAL(KIND=dp) :: t279, t280, t281, t29, t292, t3, t308, t31, t312, t314, t318, t320, t322, & 495 t33, t331, t334, t339, t34, t341, t347, t349, t35, t359, t362, t368, t37, t370, t373, & 496 t38, t382, t385, t39, t393, t4, t40, t401, t402, t404, t405, t409, t41, t42, t420, t421, & 497 t424, t425, t426, t429, t430, t431, t437, t44, t445, t446, t45, t451, t46, t473, t475, & 498 t479, t48, t484, t497, t498, t5, t50, t503, t504, t51, t517, t52, t523, t524, t527, t53, & 499 t533, t534, t54, t541, t56, t568, t57, t58, t581, t590, t6, t60, t603, t604, t61, t611, & 500 t612, t64, t640, t65, t653, t667, t675, t677, t69, t7, t71, t716 501 REAL(KIND=dp) :: t717, t720, t723, t73, t739, t758, t762, t77, t8, t800, t803, t808, t85, & 502 t86, t87, t88, t91, t92, t93, t94, t95, t98, t99 503 504 IF (order >= 0) THEN 505 t1 = 3**(0.1e1_dp/0.3e1_dp) 506 t2 = t1**2 507 t3 = ndrho*t2 508 t4 = pi**2 509 t5 = t4*rho 510 t6 = t5**(0.1e1_dp/0.3e1_dp) 511 t7 = 0.1e1_dp/t6 512 t8 = 0.1e1_dp/rho 513 ssval = t3*t7*t8/0.6e1_dp 514 t11 = red(rho, ndrho) 515 t12 = t11**2 516 t13 = sscale**2 517 t14 = t12*t13 518 t17 = t12**2 519 t19 = t13**2 520 t21 = a1*t12*t13 + a2*t17*t19 521 t22 = t14*t21 522 t25 = t17*t11 523 t27 = t19*sscale 524 t29 = t17*t12 525 t31 = t19*t13 526 t33 = 1 + a3*t17*t19 + a4*t25*t27 + a5*t29*t31 527 t34 = 0.1e1_dp/t33 528 t35 = R**2 529 t37 = t6**2 530 t38 = t2*t37 531 t39 = t34*t35*t38 532 q = t22*t39 533 t40 = t21*t34 534 t41 = 0.1e1_dp/A 535 t42 = t40*t41 536 p = 0.9e1_dp/0.4e1_dp*t14*t42 537 t44 = rho**(0.1e1_dp/0.3e1_dp) 538 t45 = t44*rho 539 t46 = exei(P, Q) 540 t48 = expint(1, q) 541 t50 = t14*t40 542 t51 = D + t50 543 t52 = t51*t35 544 t53 = t52*t38 545 t54 = expint(1, t53) 546 t56 = t51**2 547 t57 = t56*t51 548 t58 = 0.1e1_dp/t57 549 t60 = D*C 550 t61 = D**2 551 t64 = D*t21 552 t65 = t34*B 553 t69 = SQRT(pi) 554 t71 = F1*t21 555 t73 = t71*t34 + F2 556 t77 = C*(1 + t73*t12*t13) 557 t85 = t69*(15*E + 6*t77*t51 + 4*B*t56 + 8*A*t57) 558 t86 = SQRT(t51) 559 t87 = t86*t57 560 t88 = 0.1e1_dp/t87 561 t91 = SQRT(A) 562 t92 = pi*t91 563 t93 = EXP(p) 564 t94 = t11*sscale 565 t95 = SQRT(t42) 566 t98 = erf(0.3e1_dp/0.2e1_dp*t94*t95) 567 t99 = 1 - t98 568 t103 = 0.3e1_dp/0.4e1_dp*pi + t85*t88/0.16e2_dp - 0.3e1_dp/0.4e1_dp*t92 & 569 *t93*t99 570 t104 = 0.1e1_dp/t69 571 t105 = t103*t104 572 t106 = 0.1e1_dp/t12 573 t107 = 0.1e1_dp/t13 574 t108 = t106*t107 575 t109 = t108*t87 576 t113 = t40*C + REAL(2*t64*t65, KIND=dp) - 0.32e2_dp/0.15e2_dp*t105*t109 & 577 + t60*t73 578 t115 = t17*t19 579 t116 = t21**2 580 t117 = t33**2 581 t118 = 0.1e1_dp/t117 582 t119 = t116*t118 583 t121 = C*t73 584 t123 = t119*B + t40*t121 585 t125 = t35*t2 586 t126 = t61*C 587 t129 = E*t21 588 t133 = D*t103*t104 589 t136 = t34*C 590 t140 = REAL(2*t129*t34, KIND=dp) - 0.32e2_dp/0.15e2_dp*t133*t109 + REAL(2 & 591 *t64*t136, KIND=dp) + t126*t73 592 t142 = t105*t106 593 t143 = t107*t87 594 t144 = t143*t40 595 t147 = t136*t73 596 t150 = C*t116 597 t152 = -0.32e2_dp/0.15e2_dp*t142*t144 + REAL(2*t64*t147, KIND=dp) + t150 & 598 *t118 599 t155 = t29*t31*C 600 t156 = t73*t116 601 t159 = t126 + 2*D*E + t14*t140 + t115*t152 + t155*t156* & 602 t118 603 t162 = t35**2 604 t163 = t162*t1 605 t164 = t6*t5 606 t167 = t61*t103*t104 607 t170 = t34*E 608 t173 = -0.16e2_dp/0.15e2_dp*t167*t109 + REAL(2*t64*t170, KIND=dp) 609 t175 = t34*t103 610 t176 = t64*t175 611 t177 = t104*t106 612 t178 = t177*t143 613 t181 = E*t116 614 t183 = -0.32e2_dp/0.15e2_dp*t176*t178 + t181*t118 615 t187 = t104*t87*t119 616 t190 = t61*E + t14*t173 + t115*t183 - 0.16e2_dp/0.15e2_dp*t115 & 617 *t103*t187 618 t194 = 2*E + t60 + t61*B + t14*t113 + t115*t123 + t125*t37 & 619 *t159 + 3*t163*t164*t190 620 t195 = t58*t194 621 t196 = D*t35 622 t199 = EXP(-t196*t38 - q) 623 t202 = -0.4e1_dp/0.9e1_dp*A*t46 + 0.4e1_dp/0.9e1_dp*A*t48 - 0.4e1_dp/ & 624 0.9e1_dp*A*t54 - 0.4e1_dp/0.9e1_dp*t195*t199 625 e_0 = e_0 + (t45*t202*Clda)*sx 626 END IF 627 IF (order >= 1 .OR. order >= -1) THEN 628 t209 = rho**2 629 srho = -t3/t164*t8*t4/0.18e2_dp - t3*t7/t209/0.6e1_dp 630 631 t214 = t2*t7 632 sndrho = t214*t8/0.6e1_dp 633 t216 = t11*t13 634 t217 = t216*t40 635 t218 = dsdrho(rho, ndrho) 636 t222 = 2*t217*t125*t37*t218 637 t223 = a1*t11 638 t224 = t13*t218 639 t227 = t12*t11 640 t228 = a2*t227 641 t229 = t19*t218 642 t232 = 2*t223*t224 + 4*t228*t229 643 t234 = t14*t232*t39 644 t235 = t21*t118 645 t236 = t14*t235 646 t237 = a3*t227 647 t240 = a4*t17 648 t244 = a5*t25 649 t248 = 4*t237*t229 + 5*t240*t27*t218 + 6*t244*t31*t218 650 t251 = t236*t125*t37*t248 651 t255 = 0.2e1_dp/0.3e1_dp*t50*t125*t7*t4 652 qrho = t222 + t234 - t251 + t255 653 t256 = t216*t21 654 t257 = t34*t41 655 t261 = t232*t34 656 t262 = t261*t41 657 t265 = t118*t41 658 prho = 0.9e1_dp/0.2e1_dp*t256*t257*t218 + 0.9e1_dp/0.4e1_dp*t14*t262 & 659 - 0.9e1_dp/0.4e1_dp*t22*t265*t248 660 t269 = dsdndrho(rho) 661 t273 = t13*t269 662 t276 = t19*t269 663 t279 = 2*t223*t273 + 4*t228*t276 664 t280 = t279*t34 665 t281 = t280*t41 666 t292 = 4*t237*t276 + 5*t240*t27*t269 + 6*t244*t31*t269 667 pndrho = 0.9e1_dp/0.2e1_dp*t256*t257*t269 + 0.9e1_dp/0.4e1_dp*t14* & 668 t281 - 0.9e1_dp/0.4e1_dp*t22*t265*t292 669 qndrho = 2*t217*t125*t37*t269 + t14*t279*t39 - t236*t125 & 670 *t37*t292 671 t308 = dexeirho(P, Q, Prho, Qrho) 672 t312 = EXP(-q) 673 t314 = t312*t106*t107 674 t318 = 0.1e1_dp/t35 675 t320 = 0.1e1_dp/t37 676 t322 = 0.1e1_dp/t21*t33*t318*t1*t320 677 t331 = 2*t216*t40*t218 + t14*t261 - t14*t235*t248 678 t334 = t214*t4 679 t339 = EXP(-t53) 680 t341 = 0.1e1_dp/t51 681 t347 = t56**2 682 t349 = 0.1e1_dp/t347*t194 683 t359 = D*t232 684 t362 = t118*B 685 t368 = t118*t248 686 t370 = F1*t232*t34 - t71*t368 687 t373 = t73*t11 688 t382 = B*t51 689 t385 = A*t56 690 t393 = 0.1e1_dp/t86/t347 691 t401 = t92*t93 692 t402 = SQRT(0.31415926535897932385e1_dp) 693 t404 = EXP(-p) 694 t405 = 0.1e1_dp/t402*t404 695 t409 = 0.1e1_dp/t95 696 t420 = REAL(t69*(6*C*(t370*t12*t13 + 2*t373*t224)*t51 & 697 + 6*t77*t331 + 8*t382*t331 + 24*t385*t331)*t88, KIND=dp)/ & 698 0.16e2_dp - 0.7e1_dp/0.32e2_dp*REAL(t85, KIND=dp)*REAL(t393, KIND=dp)* & 699 REAL(t331, KIND=dp) - 0.3e1_dp & 700 /0.4e1_dp*t92*prho*t93*t99 + 0.3e1_dp/0.2e1_dp*t401*t405 & 701 *(0.3e1_dp/0.2e1_dp*t218*sscale*t95 + 0.3e1_dp/0.4e1_dp*t94*t409 & 702 *(t262 - t235*t41*t248)) 703 t421 = t420*t104 704 t424 = 0.1e1_dp/t227 705 t425 = t105*t424 706 t426 = t143*t218 707 t429 = t86*t56 708 t430 = t107*t429 709 t431 = t430*t331 710 t437 = t227*t19 711 t445 = 0.1e1_dp/t117/t33 712 t446 = t116*t445 713 t451 = t121*t248 714 t473 = t424*t107 715 t475 = t473*t87*t218 716 t479 = t108*t429*t331 717 t484 = t118*C 718 t497 = t105*t473 719 t498 = t87*t21 720 t503 = t105*t108 721 t504 = t429*t21 722 t517 = t64*t118 723 t523 = C*t21 724 t524 = t118*t232 725 t527 = t445*t248 726 t533 = t25*t31*C 727 t534 = t118*t218 728 t541 = t73*t21 729 t568 = t118*E 730 t581 = t64*t118*t103 731 t590 = t104*t424 732 t603 = t437*t105 733 t604 = t87*t116 734 t611 = t115*t105 735 t612 = t429*t116 736 e_rho = e_rho + (0.4e1_dp/0.3e1_dp*t44*t202*Clda + t45*(-0.4e1_dp/0.9e1_dp* & 737 A*t308 - 0.4e1_dp/0.27e2_dp*A*qrho*t314*t322 + 0.4e1_dp/0.27e2_dp & 738 *A*(t331*t35*t38 + 0.2e1_dp/0.3e1_dp*t52*t334)*t339*t341 & 739 *t318*t1*t320 + 0.4e1_dp/0.3e1_dp*t349*t199*t331 - 0.4e1_dp & 740 /0.9e1_dp*t58*(REAL(2*t216*t113*t218, KIND=dp) + t14*(t261*C & 741 - t235*C*REAL(t248, KIND=dp) + REAL(2*t359*t65, KIND=dp) - REAL(2*t64*t362 & 742 *t248, KIND=dp) - 0.32e2_dp/0.15e2_dp*t421*t109 + 0.64e2_dp/0.15e2_dp*t425 & 743 *t426 - 0.112e3_dp/0.15e2_dp*t142*t431 + t60*t370) + REAL(4* & 744 t437*t123*t218, KIND=dp) + t115*(0.2e1_dp*t235*B*t232 - 0.2e1_dp*t446 & 745 *B*REAL(t248, KIND=dp) + t261*t121 - t235*t451 + t40*C*t370) & 746 + 0.2e1_dp/0.3e1_dp*t125*t7*t159*t4 + t125*t37*(REAL(2* & 747 t216*t140*t218, KIND=dp) + t14*(0.2e1_dp*E*t232*t34 - REAL(2*t129 & 748 *t368, KIND=dp) - 0.32e2_dp/0.15e2_dp*D*t420*t104*t109 + 0.64e2_dp/0.15e2_dp & 749 *t133*t475 - 0.112e3_dp/0.15e2_dp*t133*t479 + REAL(2*t359 & 750 *t136, KIND=dp) - REAL(2*t64*t484*t248, KIND=dp) + t126*t370) + REAL(4* & 751 t437*t152*t218, KIND=dp) + t115*(-0.32e2_dp/0.15e2_dp*t421*t106*t144 & 752 + 0.64e2_dp/0.15e2_dp*t497*t498*t34*REAL(t218, KIND=dp) - 0.112e3_dp/0.15e2_dp & 753 *t503*t504*t34*t331 - 0.32e2_dp/0.15e2_dp*t142*t143* & 754 t261 + 0.32e2_dp/0.15e2_dp*t503*t498*REAL(t368, KIND=dp) + REAL(2*t359 & 755 *t147, KIND=dp) - 0.2e1_dp*t517*t451 + 0.2e1_dp*REAL(t64, KIND=dp)*REAL(t136, KIND=dp)* & 756 t370 + REAL(2*t523*t524, KIND=dp) - REAL(2*t150*t527, KIND=dp)) + REAL(6*t533 & 757 *t156*t534, KIND=dp) + t155*t370*t116*t118 + 0.2e1_dp*t155*t541 & 758 *REAL(t524, KIND=dp) - 0.2e1_dp*t155*REAL(t156, KIND=dp)*REAL(t527, KIND=dp)) + 0.4e1_dp & 759 *t163*t6*t190*t4 + 0.3e1_dp*t163*t164*(REAL(2*t216*t173 & 760 *t218, KIND=dp) + t14*(-0.16e2_dp/0.15e2_dp*t61*t420*t104*t109 + & 761 0.32e2_dp/0.15e2_dp*t167*t475 - 0.56e2_dp/0.15e2_dp*t167*t479 + REAL(2 & 762 *t359*t170, KIND=dp) - REAL(2*t64*t568*t248, KIND=dp)) + REAL(4*t437 & 763 *t183*t218, KIND=dp) + t115*(-0.32e2_dp/0.15e2_dp*REAL(t359, KIND=dp)*REAL(t175, KIND=dp) & 764 *REAL(t178, KIND=dp) + 0.32e2_dp/0.15e2_dp*t581*t177*t143*REAL(t248, KIND=dp) & 765 - 0.32e2_dp/0.15e2_dp*REAL(t64, KIND=dp)*t34*t420*REAL(t178, KIND=dp) + 0.64e2_dp & 766 /0.15e2_dp*t176*t590*t426 - 0.112e3_dp/0.15e2_dp*t176*t177* & 767 t431 + REAL(2*t129*t524, KIND=dp) - REAL(2*t181*t527, KIND=dp)) - 0.64e2_dp/ & 768 0.15e2_dp*REAL(t603, KIND=dp)*REAL(t604, KIND=dp)*REAL(t534, KIND=dp) - 0.16e2_dp/0.15e2_dp* & 769 t115*t420*t187 - 0.56e2_dp/0.15e2_dp*t611*t612*t118*t331 - & 770 0.32e2_dp/0.15e2_dp*t611*t498*REAL(t524, KIND=dp) + 0.32e2_dp/0.15e2_dp*t611 & 771 *REAL(t604, KIND=dp)*REAL(t527, KIND=dp)))*t199 - 0.4e1_dp/0.9e1_dp*t195*(-0.2e1_dp & 772 /0.3e1_dp*t196*t334 - t222 - t234 + t251 - t255)*t199)* & 773 Clda)*sx 774 t640 = dexeindrho(P, Q, Pndrho, Qndrho) 775 t653 = 2*t216*t40*t269 + t14*t280 - t14*t235*t292 776 t667 = D*t279 777 t675 = t118*t292 778 t677 = F1*t279*t34 - t71*t675 779 t716 = REAL(t69*(6*C*(t677*t12*t13 + 2*t373*t273)*t51 & 780 + 6*t77*t653 + 8*t382*t653 + 24*t385*t653)*t88, KIND=dp)/ & 781 0.16e2_dp - 0.7e1_dp/0.32e2_dp*REAL(t85, KIND=dp)*REAL(t393, KIND=dp)* & 782 REAL(t653, KIND=dp) - 0.3e1_dp & 783 /0.4e1_dp*t92*pndrho*t93*t99 + 0.3e1_dp/0.2e1_dp*t401*t405 & 784 *(0.3e1_dp/0.2e1_dp*t269*sscale*t95 + 0.3e1_dp/0.4e1_dp*t94* & 785 t409*(t281 - t235*t41*t292)) 786 t717 = t716*t104 787 t720 = t143*t269 788 t723 = t430*t653 789 t739 = t121*t292 790 t758 = t473*t87*t269 791 t762 = t108*t429*t653 792 t800 = t118*t279 793 t803 = t445*t292 794 t808 = t118*t269 795 e_ndrho = e_ndrho + (t45*(-0.4e1_dp/0.9e1_dp*A*t640 - 0.4e1_dp/0.27e2_dp*A*qndrho & 796 *t314*t322 + 0.4e1_dp/0.9e1_dp*A*t653*t339*t341 + 0.4e1_dp & 797 /0.3e1_dp*t349*t199*t653 - 0.4e1_dp/0.9e1_dp*t58*(REAL(2*t216 & 798 *t113*t269, KIND=dp) + t14*(t280*C - t235*C*REAL(t292, KIND=dp) + REAL(2 & 799 *t667*t65, KIND=dp) - REAL(2*t64*t362*t292, KIND=dp) - 0.32e2_dp/0.15e2_dp & 800 *t717*t109 + 0.64e2_dp/0.15e2_dp*t425*t720 - 0.112e3_dp/0.15e2_dp* & 801 t142*t723 + t60*t677) + REAL(4*t437*t123*t269, KIND=dp) + t115* & 802 (0.2e1_dp*t235*B*t279 - 0.2e1_dp*t446*B*REAL(t292, KIND=dp) + t280* & 803 t121 - t235*t739 + t40*C*t677) + t125*t37*(REAL(2*t216 & 804 *t140*t269, KIND=dp) + t14*(0.2e1_dp*E*t279*t34 - REAL(2*t129* & 805 t675, KIND=dp) - 0.32e2_dp/0.15e2_dp*D*t716*t104*t109 + 0.64e2_dp/0.15e2_dp & 806 *t133*t758 - 0.112e3_dp/0.15e2_dp*t133*t762 + REAL(2*t667* & 807 t136, KIND=dp) - REAL(2*t64*t484*t292, KIND=dp) + t126*t677) + REAL(4*t437 & 808 *t152*t269, KIND=dp) + t115*(-0.32e2_dp/0.15e2_dp*t717*t106*t144 + & 809 0.64e2_dp/0.15e2_dp*t497*t498*t34*REAL(t269, KIND=dp) - 0.112e3_dp/0.15e2_dp & 810 *t503*t504*t34*t653 - 0.32e2_dp/0.15e2_dp*t142*t143*t280 & 811 + 0.32e2_dp/0.15e2_dp*t503*t498*REAL(t675, KIND=dp) + REAL(2*t667* & 812 t147, KIND=dp) - 0.2e1_dp*t517*t739 + 0.2e1_dp*REAL(t64, KIND=dp)*REAL(t136, KIND=dp)*t677 & 813 + REAL(2*t523*t800, KIND=dp) - REAL(2*t150*t803, KIND=dp)) + REAL(6*t533 & 814 *t156*t808, KIND=dp) + t155*t677*t116*t118 + 0.2e1_dp*t155*t541 & 815 *REAL(t800, KIND=dp) - 0.2e1_dp*t155*REAL(t156, KIND=dp)*REAL(t803, KIND=dp)) + 0.3e1_dp*t163 & 816 *t164*(REAL(2*t216*t173*t269, KIND=dp) + t14*(-0.16e2_dp/0.15e2_dp & 817 *t61*t716*t104*t109 + 0.32e2_dp/0.15e2_dp*t167*t758 - 0.56e2_dp & 818 /0.15e2_dp*t167*t762 + REAL(2*t667*t170, KIND=dp) - REAL(2*t64 & 819 *t568*t292, KIND=dp)) + REAL(4*t437*t183*t269, KIND=dp) + t115*(-0.32e2_dp & 820 /0.15e2_dp*REAL(t667, KIND=dp)*REAL(t175, KIND=dp)*REAL(t178, KIND=dp) + 0.32e2_dp/0.15e2_dp & 821 *t581*t177*t143*REAL(t292, KIND=dp) - 0.32e2_dp/0.15e2_dp*REAL(t64, KIND=dp)* & 822 t34*t716*REAL(t178, KIND=dp) + 0.64e2_dp/0.15e2_dp*t176*t590*t720 - 0.112e3_dp & 823 /0.15e2_dp*t176*t177*t723 + REAL(2*t129*t800, KIND=dp) - REAL(2 & 824 *t181*t803, KIND=dp)) - 0.64e2_dp/0.15e2_dp*REAL(t603, KIND=dp)*REAL(t604, KIND=dp)* & 825 REAL(t808, KIND=dp) - 0.16e2_dp/0.15e2_dp*t115*t716*t187 - 0.56e2_dp/0.15e2_dp & 826 *t611*t612*t118*t653 - 0.32e2_dp/0.15e2_dp*t611*t498*REAL(t800, KIND=dp) & 827 + 0.32e2_dp/0.15e2_dp*t611*REAL(t604, KIND=dp)*REAL(t803, KIND=dp)))*t199 & 828 + 0.4e1_dp/0.9e1_dp*t195*qndrho*t199)*Clda)*sx 829 END IF 830 831 END SUBROUTINE xpbe_hole_t_c_lr_lda_calc_1 832 833! ************************************************************************************************** 834!> \brief low level routine that calculates the energy derivatives in one point 835!> \param e_0 derivatives of the energy on the grid 836!> \param e_rho derivatives of the energy on the grid 837!> \param e_ndrho derivatives of the energy on the grid 838!> \param rho value of density on the grid 839!> \param ndrho value of gradient on the grid 840!> \param sscale functional parameters 841!> \param sx functional parameters 842!> \param R functional parameters 843!> \param order order of the derivatives 844!> \par History 845!> 01.2009 created [Manuel Guidon] 846!> \author Manuel Guidon 847! ************************************************************************************************** 848 SUBROUTINE xpbe_hole_t_c_lr_lda_calc_2(e_0, e_rho, e_ndrho, & 849 rho, ndrho, sscale, sx, R, order) 850 REAL(KIND=dp), INTENT(INOUT) :: e_0, e_rho, e_ndrho 851 REAL(KIND=dp), INTENT(IN) :: rho, ndrho, sscale, sx, R 852 INTEGER, INTENT(IN) :: order 853 854 REAL(KIND=dp) :: p, pndrho, prho, q, qndrho, qrho, sndrho, srho, ssval, t1, t102, t106, t11, & 855 t110, t113, t115, t117, t118, t119, t12, t122, t125, t126, t127, t128, t13, t130, t133, & 856 t135, t138, t14, t140, t142, t143, t146, t150, t151, t152, t155, t158, t165, t17, t170, & 857 t172, t173, t174, t178, t179, t180, t183, t184, t185, t188, t19, t190, t191, t192, t193, & 858 t196, t197, t2, t200, t204, t207, t21, t211, t212, t213, t217, t22, t221, t225, t229, & 859 t232, t235, t236, t242, t248, t25, t264, t268, t27, t272, t276, t278, t280, t287, t289, & 860 t29, t292, t297, t299, t3, t305, t307, t31, t317, t320, t324 861 REAL(KIND=dp) :: t327, t33, t330, t333, t334, t338, t34, t340, t344, t35, t352, t353, t358, & 862 t37, t38, t380, t39, t394, t398, t399, t4, t40, t401, t406, t407, t408, t41, t415, t435, & 863 t44, t45, t454, t46, t461, t48, t485, t496, t498, t5, t50, t51, t512, t52, t524, t525, & 864 t529, t53, t531, t54, t545, t56, t579, t58, t581, t586, t6, t60, t61, t64, t65, t7, t74, & 865 t75, t77, t79, t8, t81, t83, t84, t85, t86, t89, t91, t93, t94, t95, t97 866 867 IF (order >= 0) THEN 868 t1 = 3**(0.1e1_dp/0.3e1_dp) 869 t2 = t1**2 870 t3 = ndrho*t2 871 t4 = pi**2 872 t5 = t4*rho 873 t6 = t5**(0.1e1_dp/0.3e1_dp) 874 t7 = 0.1e1_dp/t6 875 t8 = 0.1e1_dp/rho 876 ssval = t3*t7*t8/0.6e1_dp 877 878 t11 = red(rho, ndrho) 879 t12 = t11**2 880 t13 = sscale**2 881 t14 = t12*t13 882 t17 = t12**2 883 t19 = t13**2 884 t21 = a1*t12*t13 + a2*t17*t19 885 t22 = t14*t21 886 t25 = t17*t11 887 t27 = t19*sscale 888 t29 = t17*t12 889 t31 = t19*t13 890 t33 = 1 + a3*t17*t19 + a4*t25*t27 + a5*t29*t31 891 t34 = 0.1e1_dp/t33 892 t35 = R**2 893 t37 = t6**2 894 t38 = t2*t37 895 t39 = t34*t35*t38 896 q = t22*t39 897 t40 = t21*t34 898 t41 = 0.1e1_dp/A 899 p = 0.9e1_dp/0.4e1_dp*t14*t40*t41 900 t44 = rho**(0.1e1_dp/0.3e1_dp) 901 t45 = t44*rho 902 t46 = exei(P, Q) 903 t48 = expint(1, q) 904 t50 = t14*t40 905 t51 = D + t50 906 t52 = t51*t35 907 t53 = t52*t38 908 t54 = expint(1, t53) 909 t56 = t51**2 910 t58 = 0.1e1_dp/t56/t51 911 t60 = D*C 912 t61 = D**2 913 t64 = D*t21 914 t65 = t34*B 915 t74 = g1 + g2*t12*t13 + g3*t17*t19 + g4*t25*t27 916 t75 = E*t74 917 t77 = F1*t21 918 t79 = t77*t34 + F2 919 t81 = t40*C + 2*t64*t65 + 2*t75 + t60*t79 920 t83 = t17*t19 921 t84 = t21**2 922 t85 = t33**2 923 t86 = 0.1e1_dp/t85 924 t89 = C*t79 925 t91 = t84*t86*B + t40*t89 926 t93 = t35*t2 927 t94 = t61*C 928 t95 = D*E 929 t97 = E*t21 930 t102 = t34*C 931 t106 = 2*t97*t34 + 2*t95*t74 + 2*t64*t102 + t94*t79 932 t110 = t102*t79 933 t113 = C*t84 934 t115 = 2*t75*t40 + 2*t64*t110 + t113*t86 935 t117 = t29*t31 936 t118 = t117*C 937 t119 = t79*t84 938 t122 = t94 + 2*t95 + t14*t106 + t83*t115 + t118*t119*t86 939 t125 = t35**2 940 t126 = t125*t1 941 t127 = t6*t5 942 t128 = t61*E 943 t130 = t34*E 944 t133 = t128*t74 + 2*t64*t130 945 t135 = t130*t74 946 t138 = E*t84 947 t140 = 2*t64*t135 + t138*t86 948 t142 = t117*E 949 t143 = t74*t84 950 t146 = t128 + t14*t133 + t83*t140 + t142*t143*t86 951 t150 = 2*E + t60 + t61*B + t14*t81 + t83*t91 + t93*t37* & 952 t122 + 3*t126*t127*t146 953 t151 = t58*t150 954 t152 = D*t35 955 t155 = EXP(-t152*t38 - q) 956 t158 = -0.4e1_dp/0.9e1_dp*A*t46 + 0.4e1_dp/0.9e1_dp*A*t48 - 0.4e1_dp/ & 957 0.9e1_dp*A*t54 - 0.4e1_dp/0.9e1_dp*t151*t155 958 e_0 = e_0 + (t45*t158*Clda)*sx 959 END IF 960 IF (order >= 1 .OR. order == -1) THEN 961 t165 = rho**2 962 srho = -t3/t127*t8*t4/0.18e2_dp - t3*t7/t165/0.6e1_dp 963 t170 = t2*t7 964 sndrho = t170*t8/0.6e1_dp 965 t172 = t11*t13 966 t173 = t172*t40 967 t174 = dsdrho(rho, ndrho) 968 t178 = 2*t173*t93*t37*t174 969 t179 = a1*t11 970 t180 = t13*t174 971 t183 = t12*t11 972 t184 = a2*t183 973 t185 = t19*t174 974 t188 = 2*t179*t180 + 4*t184*t185 975 t190 = t14*t188*t39 976 t191 = t21*t86 977 t192 = t14*t191 978 t193 = a3*t183 979 t196 = a4*t17 980 t197 = t27*t174 981 t200 = a5*t25 982 t204 = 4*t193*t185 + 5*t196*t197 + 6*t200*t31*t174 983 t207 = t192*t93*t37*t204 984 t211 = 0.2e1_dp/0.3e1_dp*t50*t93*t7*t4 985 qrho = t178 + t190 - t207 + t211 986 t212 = t172*t21 987 t213 = t34*t41 988 t217 = t188*t34 989 t221 = t86*t41 990 prho = 0.9e1_dp/0.2e1_dp*t212*t213*t174 + 0.9e1_dp/0.4e1_dp*t14*t217 & 991 *t41 - 0.9e1_dp/0.4e1_dp*t22*t221*t204 992 t225 = dsdndrho(rho) 993 t229 = t13*t225 994 t232 = t19*t225 995 t235 = 2*t179*t229 + 4*t184*t232 996 t236 = t235*t34 997 t242 = t27*t225 998 t248 = 4*t193*t232 + 5*t196*t242 + 6*t200*t31*t225 999 pndrho = 0.9e1_dp/0.2e1_dp*t212*t213*t225 + 0.9e1_dp/0.4e1_dp*t14* & 1000 t236*t41 - 0.9e1_dp/0.4e1_dp*t22*t221*t248 1001 qndrho = 2*t173*t93*t37*t225 + t14*t235*t39 - t192*t93 & 1002 *t37*t248 1003 t264 = dexeirho(P, Q, Prho, Qrho) 1004 t268 = EXP(-q) 1005 t272 = t268/t12/t13 1006 t276 = 0.1e1_dp/t35 1007 t278 = 0.1e1_dp/t37 1008 t280 = 0.1e1_dp/t21*t33*t276*t1*t278 1009 t287 = t191*t204 1010 t289 = 2*t172*t40*t174 + t14*t217 - t14*t287 1011 t292 = t170*t4 1012 t297 = EXP(-t53) 1013 t299 = 0.1e1_dp/t51 1014 t305 = t56**2 1015 t307 = 0.1e1_dp/t305*t150 1016 t317 = D*t188 1017 t320 = t86*B 1018 t324 = g2*t11 1019 t327 = g3*t183 1020 t330 = g4*t17 1021 t333 = 2*t324*t180 + 4*t327*t185 + 5*t330*t197 1022 t334 = E*t333 1023 t338 = t86*t204 1024 t340 = F1*t188*t34 - t77*t338 1025 t344 = t183*t19 1026 t352 = 0.1e1_dp/t85/t33 1027 t353 = t84*t352 1028 t358 = t89*t204 1029 t380 = t86*C 1030 t394 = t64*t86 1031 t398 = C*t21 1032 t399 = t86*t188 1033 t401 = t352*t204 1034 t406 = t25*t31 1035 t407 = t406*C 1036 t408 = t86*t174 1037 t415 = t79*t21 1038 t435 = t86*E 1039 t454 = t406*E 1040 t461 = t74*t21 1041 e_rho = e_rho + (0.4e1_dp/0.3e1_dp*t44*t158*Clda + t45*(-0.4e1_dp/0.9e1_dp* & 1042 A*t264 - 0.4e1_dp/0.27e2_dp*A*qrho*t272*t280 + 0.4e1_dp/0.27e2_dp & 1043 *A*(t289*t35*t38 + 0.2e1_dp/0.3e1_dp*t52*t292)*t297*t299 & 1044 *t276*t1*t278 + 0.4e1_dp/0.3e1_dp*t307*t155*t289 - 0.4e1_dp & 1045 /0.9e1_dp*t58*(REAL(2*t172*t81*t174, KIND=dp) + REAL(t14*(t217 & 1046 *C - t191*C*t204 + 2*t317*t65 - 2*t64*t320*t204 + 2 & 1047 *t334 + t60*t340), KIND=dp) + REAL(4*t344*t91*t174, KIND=dp) + REAL(t83* & 1048 (2*t191*B*t188 - 2*t353*B*t204 + t217*t89 - t191*t358 & 1049 + t40*C*t340), KIND=dp) + 0.2e1_dp/0.3e1_dp*t93*t7*t122*t4 + t93 & 1050 *t37*REAL(2*t172*t106*t174 + t14*(2*E*t188*t34 & 1051 - 2*t97*t338 + 2*t95*t333 + 2*t317*t102 - 2*t64*t380 & 1052 *t204 + t94*t340) + 4*t344*t115*t174 + t83*(2*t334 & 1053 *t40 + 2*t75*t217 - 2*t75*t287 + 2*t317*t110 - 2*t394 & 1054 *t358 + 2*t64*t102*t340 + 2*t398*t399 - 2*t113* & 1055 t401) + 6*t407*t119*t408 + t118*t340*t84*t86 + 2*t118 & 1056 *t415*t399 - 2*t118*t119*t401, KIND=dp) + 0.4e1_dp*t126*t6*t146 & 1057 *t4 + 0.3e1_dp*t126*t127*REAL(2*t172*t133*t174 + t14 & 1058 *(t128*t333 + 2*t317*t130 - 2*t64*t435*t204) + 4*t344 & 1059 *t140*t174 + t83*(2*t317*t135 - 2*t394*t75*t204 & 1060 + 2*t64*t130*t333 + 2*t97*t399 - 2*t138*t401) + 6* & 1061 t454*t143*t408 + t142*t333*t84*t86 + 2*t142*t461*t399 & 1062 - 2*t142*t143*t401, KIND=dp))*t155 - 0.4e1_dp/0.9e1_dp*t151*(-0.2e1_dp & 1063 /0.3e1_dp*t152*t292 - t178 - t190 + t207 - t211)*t155)* & 1064 Clda)*sx 1065 t485 = dexeindrho(P, Q, Pndrho, Qndrho) 1066 t496 = t191*t248 1067 t498 = 2*t172*t40*t225 + t14*t236 - t14*t496 1068 t512 = D*t235 1069 t524 = 2*t324*t229 + 4*t327*t232 + 5*t330*t242 1070 t525 = E*t524 1071 t529 = t86*t248 1072 t531 = F1*t235*t34 - t77*t529 1073 t545 = t89*t248 1074 t579 = t86*t235 1075 t581 = t352*t248 1076 t586 = t86*t225 1077 e_ndrho = e_ndrho + (t45*(-0.4e1_dp/0.9e1_dp*A*t485 - 0.4e1_dp/0.27e2_dp*A*qndrho & 1078 *t272*t280 + 0.4e1_dp/0.9e1_dp*A*t498*t297*t299 + 0.4e1_dp & 1079 /0.3e1_dp*t307*t155*t498 - 0.4e1_dp/0.9e1_dp*t58*REAL(2*t172 & 1080 *t81*t225 + t14*(t236*C - t191*C*t248 + 2*t512*t65 & 1081 - 2*t64*t320*t248 + 2*t525 + t60*t531) + 4*t344*t91 & 1082 *t225 + t83*(2*t191*B*t235 - 2*t353*B*t248 + t236 & 1083 *t89 - t191*t545 + t40*C*t531) + t93*t37*(2*t172* & 1084 t106*t225 + t14*(2*E*t235*t34 - 2*t97*t529 + 2*t95 & 1085 *t524 + 2*t512*t102 - 2*t64*t380*t248 + t94*t531) + & 1086 4*t344*t115*t225 + t83*(2*t525*t40 + 2*t75*t236 - & 1087 2*t75*t496 + 2*t512*t110 - 2*t394*t545 + 2*t64*t102 & 1088 *t531 + 2*t398*t579 - 2*t113*t581) + 6*t407*t119* & 1089 t586 + t118*t531*t84*t86 + 2*t118*t415*t579 - 2*t118 & 1090 *t119*t581) + 3*t126*t127*(2*t172*t133*t225 + t14 & 1091 *(t128*t524 + 2*t512*t130 - 2*t64*t435*t248) + 4*t344 & 1092 *t140*t225 + t83*(2*t512*t135 - 2*t394*t75*t248 & 1093 + 2*t64*t130*t524 + 2*t97*t579 - 2*t138*t581) + 6* & 1094 t454*t143*t586 + t142*t524*t84*t86 + 2*t142*t461*t579 & 1095 - 2*t142*t143*t581), KIND=dp)*t155 + 0.4e1_dp/0.9e1_dp*t151*qndrho & 1096 *t155)*Clda)*sx 1097 END IF 1098 1099 END SUBROUTINE xpbe_hole_t_c_lr_lda_calc_2 1100 1101! ************************************************************************************************** 1102!> \brief These functions evaluate products exp(x)*Ei(x) and pi*exp(x)*erfc(sqrt(x)), 1103!> as well as their derivatives with respect to various combinations of 1104!> rho and norm_drho. 1105!> \param P = 9/4*s^2*H/A 1106!> \param Q = s^2*H*R^2*kf^2 1107!> \return ... 1108!> \par History 1109!> 01.2009 created [Manuel Guidon] 1110!> \author Manuel Guidon 1111!> \note 1112!> - In order to avoid numerical instabilities, these routines use Taylor- 1113!> expansions for the above core-products for large arguments. 1114!> - red calculates the reduced gradient in a save way (i.e. using a fixed 1115!> cutoff EPS_RHO) 1116! ************************************************************************************************** 1117 FUNCTION exei(P, Q) 1118 REAL(dp), INTENT(IN) :: P, Q 1119 REAL(dp) :: exei 1120 1121 REAL(dp) :: Q2, Q3, Q4, tmp 1122 1123 exei = 0.0_dp 1124 IF (P < expcutoff) THEN 1125 !Use exact product 1126 IF (P + Q < 0.5_dp) THEN 1127 tmp = -euler - LOG(P + Q) + P + Q 1128 tmp = tmp - 0.25_dp*(P + Q)**2 + 1.0_dp/18.0_dp*(P + Q)**3 - 1.0_dp/96.0_dp*(P + Q)**4 1129 tmp = tmp + 1.0_dp/600.0_dp*(P + Q)**5 1130 exei = EXP(P)*tmp 1131 ELSE 1132 exei = EXP(P)*expint(1, Q + P) 1133 END IF 1134 ELSE 1135 !Use approximation 1136 tmp = 1.0_dp/P 1137 ! *** 1st order 1138 exei = tmp 1139 ! *** 2nd order 1140 tmp = tmp/P 1141 exei = exei - (Q + 1.0_dp)*tmp 1142 ! *** 3rd order 1143 tmp = tmp/P 1144 Q2 = Q*Q 1145 exei = exei + (2.0_dp*Q + Q2 + 2.0_dp)*tmp 1146 ! *** 4th order 1147 tmp = tmp/P 1148 Q3 = Q2*Q 1149 exei = exei - (3.0_dp*Q2 + 6.0_dp*Q + Q3 + 6.0_dp)*tmp 1150 ! *** 5th order 1151 tmp = tmp/P 1152 Q4 = Q3*Q 1153 exei = exei + (24.0_dp + 4.0_dp*Q3 + Q4 + 12.0_dp*Q2 + 24.0_dp*Q)*tmp 1154 1155 ! *** scaling factor 1156 exei = EXP(-Q)*exei 1157 END IF 1158 END FUNCTION exei 1159 1160! ************************************************************************************************** 1161!> \brief ... 1162!> \param P ... 1163!> \param Q ... 1164!> \param dPrho ... 1165!> \param dQrho ... 1166!> \return ... 1167! ************************************************************************************************** 1168 FUNCTION dexeirho(P, Q, dPrho, dQrho) 1169 REAL(dp), INTENT(IN) :: P, Q, dPrho, dQrho 1170 REAL(dp) :: dexeirho 1171 1172 dexeirho = dPrho*exei(P, Q) - (dPrho + dQrho)/(P + Q)*EXP(-Q) 1173 END FUNCTION dexeirho 1174 1175! ************************************************************************************************** 1176!> \brief ... 1177!> \param P ... 1178!> \param Q ... 1179!> \param dPndrho ... 1180!> \param dQndrho ... 1181!> \return ... 1182! ************************************************************************************************** 1183 FUNCTION dexeindrho(P, Q, dPndrho, dQndrho) 1184 REAL(dp), INTENT(IN) :: P, Q, dPndrho, dQndrho 1185 REAL(dp) :: dexeindrho 1186 1187 dexeindrho = dPndrho*exei(P, Q) - (dPndrho + dQndrho)/(P + Q)*EXP(-Q) 1188 END FUNCTION dexeindrho 1189 1190! ************************************************************************************************** 1191!> \brief ... 1192!> \param rho ... 1193!> \param ndrho ... 1194!> \return ... 1195! ************************************************************************************************** 1196 FUNCTION red(rho, ndrho) 1197 REAL(dp), INTENT(IN) :: rho, ndrho 1198 REAL(dp) :: red 1199 1200 red = 1.0_dp/6.0_dp*ndrho*3.0_dp**(2.0_dp/3.0_dp) 1201 red = red/(pi**(2.0_dp/3.0_dp)) 1202 red = red*MAX(1.0_dp/rho**(4.0_dp/3.0_dp), EPS_RHO) 1203 END FUNCTION red 1204 1205! ************************************************************************************************** 1206!> \brief ... 1207!> \param rho ... 1208!> \param ndrho ... 1209!> \return ... 1210! ************************************************************************************************** 1211 FUNCTION dsdrho(rho, ndrho) 1212 REAL(dp), INTENT(IN) :: rho, ndrho 1213 REAL(dp) :: dsdrho 1214 1215 dsdrho = -2.0_dp/9.0_dp*ndrho*3.0**(2.0_dp/3.0_dp) 1216 dsdrho = dsdrho/(pi**(2.0_dp/3.0_dp)) 1217 dsdrho = dsdrho*MAX(1.0_dp/rho**(7.0_dp/3.0_dp), EPS_RHO) 1218 END FUNCTION dsdrho 1219 1220! ************************************************************************************************** 1221!> \brief ... 1222!> \param rho ... 1223!> \return ... 1224! ************************************************************************************************** 1225 FUNCTION dsdndrho(rho) 1226 REAL(dp), INTENT(IN) :: rho 1227 REAL(dp) :: dsdndrho 1228 1229 dsdndrho = 1.0_dp/6.0_dp*3.0_dp**(2.0_dp/3.0_dp) 1230 dsdndrho = dsdndrho/(pi**(2.0_dp/3.0_dp)) 1231 dsdndrho = dsdndrho*MAX(1.0_dp/rho**(4.0_dp/3.0_dp), EPS_RHO) 1232 END FUNCTION dsdndrho 1233 1234END MODULE xc_xpbe_hole_t_c_lr 1235 1236