1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Calculate the Thomas-Fermi kinetic energy functional 8!> plus the von Weizsaecker term 9!> \par History 10!> JGH (26.02.2003) : OpenMP enabled 11!> fawzi (04.2004) : adapted to the new xc interface 12!> \author JGH (18.02.2002) 13! ************************************************************************************************** 14MODULE xc_tfw 15 USE cp_array_utils, ONLY: cp_3d_r_p_type 16 USE kinds, ONLY: dp 17 USE xc_derivative_set_types, ONLY: xc_derivative_set_type,& 18 xc_dset_get_derivative 19 USE xc_derivative_types, ONLY: xc_derivative_get,& 20 xc_derivative_type 21 USE xc_functionals_utilities, ONLY: set_util 22 USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type 23 USE xc_rho_set_types, ONLY: xc_rho_set_get,& 24 xc_rho_set_type 25#include "../base/base_uses.f90" 26 27 IMPLICIT NONE 28 29 PRIVATE 30 31! *** Global parameters *** 32 33 REAL(KIND=dp), PARAMETER :: pi = 3.14159265358979323846264338_dp 34 REAL(KIND=dp), PARAMETER :: f13 = 1.0_dp/3.0_dp, & 35 f23 = 2.0_dp*f13, & 36 f43 = 4.0_dp*f13, & 37 f53 = 5.0_dp*f13 38 39 PUBLIC :: tfw_lda_info, tfw_lda_eval, tfw_lsd_info, tfw_lsd_eval 40 41 REAL(KIND=dp) :: cf, flda, flsd, fvw 42 REAL(KIND=dp) :: eps_rho 43 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_tfw' 44 45CONTAINS 46 47! ************************************************************************************************** 48!> \brief ... 49!> \param cutoff ... 50! ************************************************************************************************** 51 SUBROUTINE tfw_init(cutoff) 52 53 REAL(KIND=dp), INTENT(IN) :: cutoff 54 55 eps_rho = cutoff 56 CALL set_util(cutoff) 57 58 cf = 0.3_dp*(3.0_dp*pi*pi)**f23 59 flda = cf 60 flsd = flda*2.0_dp**f23 61 fvw = 1.0_dp/72.0_dp 62 63 END SUBROUTINE tfw_init 64 65! ************************************************************************************************** 66!> \brief ... 67!> \param reference ... 68!> \param shortform ... 69!> \param needs ... 70!> \param max_deriv ... 71! ************************************************************************************************** 72 SUBROUTINE tfw_lda_info(reference, shortform, needs, max_deriv) 73 CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform 74 TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs 75 INTEGER, INTENT(out), OPTIONAL :: max_deriv 76 77 IF (PRESENT(reference)) THEN 78 reference = "Thomas-Fermi-Weizsaecker kinetic energy functional {LDA version}" 79 END IF 80 IF (PRESENT(shortform)) THEN 81 shortform = "TF+vW kinetic energy functional {LDA}" 82 END IF 83 IF (PRESENT(needs)) THEN 84 needs%rho = .TRUE. 85 needs%rho_1_3 = .TRUE. 86 needs%norm_drho = .TRUE. 87 END IF 88 IF (PRESENT(max_deriv)) max_deriv = 3 89 90 END SUBROUTINE tfw_lda_info 91 92! ************************************************************************************************** 93!> \brief ... 94!> \param reference ... 95!> \param shortform ... 96!> \param needs ... 97!> \param max_deriv ... 98! ************************************************************************************************** 99 SUBROUTINE tfw_lsd_info(reference, shortform, needs, max_deriv) 100 CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform 101 TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs 102 INTEGER, INTENT(out), OPTIONAL :: max_deriv 103 104 IF (PRESENT(reference)) THEN 105 reference = "Thomas-Fermi-Weizsaecker kinetic energy functional" 106 END IF 107 IF (PRESENT(shortform)) THEN 108 shortform = "TF+vW kinetic energy functional" 109 END IF 110 IF (PRESENT(needs)) THEN 111 needs%rho_spin = .TRUE. 112 needs%rho_spin_1_3 = .TRUE. 113 needs%norm_drho = .TRUE. 114 END IF 115 IF (PRESENT(max_deriv)) max_deriv = 3 116 117 END SUBROUTINE tfw_lsd_info 118 119! ************************************************************************************************** 120!> \brief ... 121!> \param rho_set ... 122!> \param deriv_set ... 123!> \param order ... 124! ************************************************************************************************** 125 SUBROUTINE tfw_lda_eval(rho_set, deriv_set, order) 126 TYPE(xc_rho_set_type), POINTER :: rho_set 127 TYPE(xc_derivative_set_type), POINTER :: deriv_set 128 INTEGER, INTENT(in) :: order 129 130 CHARACTER(len=*), PARAMETER :: routineN = 'tfw_lda_eval', routineP = moduleN//':'//routineN 131 132 INTEGER :: handle, npoints 133 INTEGER, DIMENSION(:, :), POINTER :: bo 134 REAL(KIND=dp) :: epsilon_rho 135 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: s 136 REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: e_0, e_ndrho, e_ndrho_ndrho, e_rho, & 137 e_rho_ndrho, e_rho_ndrho_ndrho, e_rho_rho, e_rho_rho_ndrho, e_rho_rho_rho, grho, r13, rho 138 TYPE(xc_derivative_type), POINTER :: deriv 139 140 CALL timeset(routineN, handle) 141 NULLIFY (bo) 142 143 CPASSERT(ASSOCIATED(rho_set)) 144 CPASSERT(rho_set%ref_count > 0) 145 CPASSERT(ASSOCIATED(deriv_set)) 146 CPASSERT(deriv_set%ref_count > 0) 147 CALL xc_rho_set_get(rho_set, rho_1_3=r13, rho=rho, & 148 norm_drho=grho, local_bounds=bo, rho_cutoff=epsilon_rho) 149 npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1) 150 CALL tfw_init(epsilon_rho) 151 152 ALLOCATE (s(npoints)) 153 CALL calc_s(rho, grho, s, npoints) 154 155 IF (order >= 0) THEN 156 deriv => xc_dset_get_derivative(deriv_set, "", & 157 allocate_deriv=.TRUE.) 158 CALL xc_derivative_get(deriv, deriv_data=e_0) 159 160 CALL tfw_u_0(rho, r13, s, e_0, npoints) 161 END IF 162 IF (order >= 1 .OR. order == -1) THEN 163 deriv => xc_dset_get_derivative(deriv_set, "(rho)", & 164 allocate_deriv=.TRUE.) 165 CALL xc_derivative_get(deriv, deriv_data=e_rho) 166 deriv => xc_dset_get_derivative(deriv_set, "(norm_drho)", & 167 allocate_deriv=.TRUE.) 168 CALL xc_derivative_get(deriv, deriv_data=e_ndrho) 169 170 CALL tfw_u_1(rho, grho, r13, s, e_rho, e_ndrho, npoints) 171 END IF 172 IF (order >= 2 .OR. order == -2) THEN 173 deriv => xc_dset_get_derivative(deriv_set, "(rho)(rho)", & 174 allocate_deriv=.TRUE.) 175 CALL xc_derivative_get(deriv, deriv_data=e_rho_rho) 176 deriv => xc_dset_get_derivative(deriv_set, "(rho)(norm_drho)", & 177 allocate_deriv=.TRUE.) 178 CALL xc_derivative_get(deriv, deriv_data=e_rho_ndrho) 179 deriv => xc_dset_get_derivative(deriv_set, & 180 "(norm_drho)(norm_drho)", allocate_deriv=.TRUE.) 181 CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho) 182 183 CALL tfw_u_2(rho, grho, r13, s, e_rho_rho, e_rho_ndrho, & 184 e_ndrho_ndrho, npoints) 185 END IF 186 IF (order >= 3 .OR. order == -3) THEN 187 deriv => xc_dset_get_derivative(deriv_set, "(rho)(rho)(rho)", & 188 allocate_deriv=.TRUE.) 189 CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho) 190 deriv => xc_dset_get_derivative(deriv_set, & 191 "(rho)(rho)(norm_drho)", allocate_deriv=.TRUE.) 192 CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_ndrho) 193 deriv => xc_dset_get_derivative(deriv_set, & 194 "(rho)(norm_drho)(norm_drho)", allocate_deriv=.TRUE.) 195 CALL xc_derivative_get(deriv, deriv_data=e_rho_ndrho_ndrho) 196 197 CALL tfw_u_3(rho, grho, r13, s, e_rho_rho_rho, e_rho_rho_ndrho, & 198 e_rho_ndrho_ndrho, npoints) 199 END IF 200 IF (order > 3 .OR. order < -3) THEN 201 CPABORT("derivatives bigger than 3 not implemented") 202 END IF 203 204 DEALLOCATE (s) 205 CALL timestop(handle) 206 END SUBROUTINE tfw_lda_eval 207 208! ************************************************************************************************** 209!> \brief ... 210!> \param rho ... 211!> \param grho ... 212!> \param s ... 213!> \param npoints ... 214! ************************************************************************************************** 215 SUBROUTINE calc_s(rho, grho, s, npoints) 216 REAL(KIND=dp), DIMENSION(*), INTENT(in) :: rho, grho 217 REAL(KIND=dp), DIMENSION(*), INTENT(out) :: s 218 INTEGER, INTENT(in) :: npoints 219 220 INTEGER :: ip 221 222!$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)& 223!$OMP SHARED(npoints,rho,eps_rho,s,grho) 224 DO ip = 1, npoints 225 IF (rho(ip) < eps_rho) THEN 226 s(ip) = 0.0_dp 227 ELSE 228 s(ip) = grho(ip)*grho(ip)/rho(ip) 229 END IF 230 END DO 231 END SUBROUTINE calc_s 232 233! ************************************************************************************************** 234!> \brief ... 235!> \param rho_set ... 236!> \param deriv_set ... 237!> \param order ... 238! ************************************************************************************************** 239 SUBROUTINE tfw_lsd_eval(rho_set, deriv_set, order) 240 TYPE(xc_rho_set_type), POINTER :: rho_set 241 TYPE(xc_derivative_set_type), POINTER :: deriv_set 242 INTEGER, INTENT(in) :: order 243 244 CHARACTER(len=*), PARAMETER :: routineN = 'tfw_lsd_eval', routineP = moduleN//':'//routineN 245 CHARACTER(len=12), DIMENSION(2), PARAMETER :: & 246 norm_drho_spin_name = (/"(norm_drhoa)", "(norm_drhob)"/) 247 CHARACTER(len=6), DIMENSION(2), PARAMETER :: rho_spin_name = (/"(rhoa)", "(rhob)"/) 248 249 INTEGER :: handle, i, ispin, npoints 250 INTEGER, DIMENSION(:, :), POINTER :: bo 251 REAL(KIND=dp) :: epsilon_rho 252 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: s 253 REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: e_0, e_ndrho, e_ndrho_ndrho, e_rho, & 254 e_rho_ndrho, e_rho_ndrho_ndrho, & 255 e_rho_rho, e_rho_rho_ndrho, & 256 e_rho_rho_rho 257 TYPE(cp_3d_r_p_type), DIMENSION(2) :: norm_drho, rho, rho_1_3 258 TYPE(xc_derivative_type), POINTER :: deriv 259 260 CALL timeset(routineN, handle) 261 NULLIFY (deriv, bo) 262 DO i = 1, 2 263 NULLIFY (norm_drho(i)%array, rho(i)%array, rho_1_3(i)%array) 264 END DO 265 266 CPASSERT(ASSOCIATED(rho_set)) 267 CPASSERT(rho_set%ref_count > 0) 268 CPASSERT(ASSOCIATED(deriv_set)) 269 CPASSERT(deriv_set%ref_count > 0) 270 CALL xc_rho_set_get(rho_set, rhoa_1_3=rho_1_3(1)%array, & 271 rhob_1_3=rho_1_3(2)%array, rhoa=rho(1)%array, & 272 rhob=rho(2)%array, norm_drhoa=norm_drho(1)%array, & 273 norm_drhob=norm_drho(2)%array, rho_cutoff=epsilon_rho, & 274 local_bounds=bo) 275 npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1) 276 CALL tfw_init(epsilon_rho) 277 278 ALLOCATE (s(npoints)) 279 280 DO ispin = 1, 2 281 CALL calc_s(rho(ispin)%array, norm_drho(ispin)%array, s, npoints) 282 283 IF (order >= 0) THEN 284 deriv => xc_dset_get_derivative(deriv_set, "", & 285 allocate_deriv=.TRUE.) 286 CALL xc_derivative_get(deriv, deriv_data=e_0) 287 288 CALL tfw_p_0(rho(ispin)%array, & 289 rho_1_3(ispin)%array, s, e_0, npoints) 290 END IF 291 IF (order >= 1 .OR. order == -1) THEN 292 deriv => xc_dset_get_derivative(deriv_set, rho_spin_name(ispin), & 293 allocate_deriv=.TRUE.) 294 CALL xc_derivative_get(deriv, deriv_data=e_rho) 295 deriv => xc_dset_get_derivative(deriv_set, norm_drho_spin_name(ispin), & 296 allocate_deriv=.TRUE.) 297 CALL xc_derivative_get(deriv, deriv_data=e_ndrho) 298 299 CALL tfw_p_1(rho(ispin)%array, norm_drho(ispin)%array, & 300 rho_1_3(ispin)%array, s, e_rho, e_ndrho, npoints) 301 END IF 302 IF (order >= 2 .OR. order == -2) THEN 303 deriv => xc_dset_get_derivative(deriv_set, rho_spin_name(ispin)// & 304 rho_spin_name(ispin), allocate_deriv=.TRUE.) 305 CALL xc_derivative_get(deriv, deriv_data=e_rho_rho) 306 deriv => xc_dset_get_derivative(deriv_set, rho_spin_name(ispin)// & 307 norm_drho_spin_name(ispin), allocate_deriv=.TRUE.) 308 CALL xc_derivative_get(deriv, deriv_data=e_rho_ndrho) 309 deriv => xc_dset_get_derivative(deriv_set, norm_drho_spin_name(ispin)// & 310 norm_drho_spin_name(ispin), allocate_deriv=.TRUE.) 311 CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho) 312 313 CALL tfw_p_2(rho(ispin)%array, norm_drho(ispin)%array, & 314 rho_1_3(ispin)%array, s, e_rho_rho, e_rho_ndrho, & 315 e_ndrho_ndrho, npoints) 316 END IF 317 IF (order >= 3 .OR. order == -3) THEN 318 deriv => xc_dset_get_derivative(deriv_set, rho_spin_name(ispin)// & 319 rho_spin_name(ispin)//rho_spin_name(ispin), & 320 allocate_deriv=.TRUE.) 321 CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho) 322 deriv => xc_dset_get_derivative(deriv_set, rho_spin_name(ispin)// & 323 rho_spin_name(ispin)//norm_drho_spin_name(ispin), & 324 allocate_deriv=.TRUE.) 325 CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_ndrho) 326 deriv => xc_dset_get_derivative(deriv_set, rho_spin_name(ispin)// & 327 norm_drho_spin_name(ispin)//norm_drho_spin_name(ispin), & 328 allocate_deriv=.TRUE.) 329 CALL xc_derivative_get(deriv, deriv_data=e_rho_ndrho_ndrho) 330 331 CALL tfw_p_3(rho(ispin)%array, norm_drho(ispin)%array, & 332 rho_1_3(ispin)%array, s, e_rho_rho_rho, e_rho_rho_ndrho, & 333 e_rho_ndrho_ndrho, npoints) 334 END IF 335 IF (order > 3 .OR. order < -3) THEN 336 CPABORT("derivatives bigger than 3 not implemented") 337 END IF 338 END DO 339 340 DEALLOCATE (s) 341 CALL timestop(handle) 342 END SUBROUTINE tfw_lsd_eval 343 344! ************************************************************************************************** 345!> \brief ... 346!> \param rho ... 347!> \param r13 ... 348!> \param s ... 349!> \param e_0 ... 350!> \param npoints ... 351! ************************************************************************************************** 352 SUBROUTINE tfw_u_0(rho, r13, s, e_0, npoints) 353 354 REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, r13, s 355 REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_0 356 INTEGER, INTENT(in) :: npoints 357 358 INTEGER :: ip 359 360!$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)& 361!$OMP SHARED(npoints,rho,eps_rho,e_0,flda,r13,s,fvw) 362 DO ip = 1, npoints 363 364 IF (rho(ip) > eps_rho) THEN 365 366 e_0(ip) = e_0(ip) + flda*r13(ip)*r13(ip)*rho(ip) + fvw*s(ip) 367 368 END IF 369 370 END DO 371 372 END SUBROUTINE tfw_u_0 373 374! ************************************************************************************************** 375!> \brief ... 376!> \param rho ... 377!> \param grho ... 378!> \param r13 ... 379!> \param s ... 380!> \param e_rho ... 381!> \param e_ndrho ... 382!> \param npoints ... 383! ************************************************************************************************** 384 SUBROUTINE tfw_u_1(rho, grho, r13, s, e_rho, e_ndrho, npoints) 385 386 REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, grho, r13, s 387 REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_rho, e_ndrho 388 INTEGER, INTENT(in) :: npoints 389 390 INTEGER :: ip 391 REAL(KIND=dp) :: f 392 393 f = f53*flda 394 395!$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)& 396!$OMP SHARED(npoints,rho,eps_rho,e_rho,e_ndrho,grho,s,r13,f,fvw) 397 DO ip = 1, npoints 398 399 IF (rho(ip) > eps_rho) THEN 400 401 e_rho(ip) = e_rho(ip) + f*r13(ip)*r13(ip) - fvw*s(ip)/rho(ip) 402 e_ndrho(ip) = e_ndrho(ip) + 2.0_dp*fvw*grho(ip)/rho(ip) 403 404 END IF 405 406 END DO 407 408 END SUBROUTINE tfw_u_1 409 410! ************************************************************************************************** 411!> \brief ... 412!> \param rho ... 413!> \param grho ... 414!> \param r13 ... 415!> \param s ... 416!> \param e_rho_rho ... 417!> \param e_rho_ndrho ... 418!> \param e_ndrho_ndrho ... 419!> \param npoints ... 420! ************************************************************************************************** 421 SUBROUTINE tfw_u_2(rho, grho, r13, s, e_rho_rho, e_rho_ndrho, e_ndrho_ndrho, & 422 npoints) 423 424 REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, grho, r13, s 425 REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_rho_rho, e_rho_ndrho, e_ndrho_ndrho 426 INTEGER, INTENT(in) :: npoints 427 428 INTEGER :: ip 429 REAL(KIND=dp) :: f 430 431 f = f23*f53*flda 432 433!$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)& 434!$OMP SHARED(npoints,rho,eps_rho,e_rho_rho,e_rho_ndrho,e_ndrho_ndrho,grho,f,fvw) 435 DO ip = 1, npoints 436 437 IF (rho(ip) > eps_rho) THEN 438 439 e_rho_rho(ip) = e_rho_rho(ip) + f/r13(ip) + 2.0_dp*fvw*s(ip)/(rho(ip)*rho(ip)) 440 e_rho_ndrho(ip) = e_rho_ndrho(ip) - 2.0_dp*fvw*grho(ip)/(rho(ip)*rho(ip)) 441 e_ndrho_ndrho(ip) = e_ndrho_ndrho(ip) + 2.0_dp*fvw/rho(ip) 442 443 END IF 444 445 END DO 446 447 END SUBROUTINE tfw_u_2 448 449! ************************************************************************************************** 450!> \brief ... 451!> \param rho ... 452!> \param grho ... 453!> \param r13 ... 454!> \param s ... 455!> \param e_rho_rho_rho ... 456!> \param e_rho_rho_ndrho ... 457!> \param e_rho_ndrho_ndrho ... 458!> \param npoints ... 459! ************************************************************************************************** 460 SUBROUTINE tfw_u_3(rho, grho, r13, s, e_rho_rho_rho, e_rho_rho_ndrho, & 461 e_rho_ndrho_ndrho, npoints) 462 463 REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, grho, r13, s 464 REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_rho_rho_rho, e_rho_rho_ndrho, & 465 e_rho_ndrho_ndrho 466 INTEGER, INTENT(in) :: npoints 467 468 INTEGER :: ip 469 REAL(KIND=dp) :: f 470 471 f = -f13*f23*f53*flda 472 473!$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)& 474!$OMP SHARED(npoints,rho,eps_rho,e_rho_rho_rho,r13,s,e_rho_rho_ndrho,e_rho_ndrho_ndrho,f,fvw) 475 DO ip = 1, npoints 476 477 IF (rho(ip) > eps_rho) THEN 478 479 e_rho_rho_rho(ip) = e_rho_rho_rho(ip) + f/(r13(ip)*rho(ip)) & 480 - 6.0_dp*fvw*s(ip)/(rho(ip)*rho(ip)*rho(ip)) 481 e_rho_rho_ndrho(ip) = e_rho_rho_ndrho(ip) & 482 + 4.0_dp*fvw*grho(ip)/(rho(ip)*rho(ip)*rho(ip)) 483 e_rho_ndrho_ndrho(ip) = e_rho_ndrho_ndrho(ip) & 484 - 2.0_dp*fvw/(rho(ip)*rho(ip)) 485 END IF 486 487 END DO 488 489 END SUBROUTINE tfw_u_3 490 491! ************************************************************************************************** 492!> \brief ... 493!> \param rhoa ... 494!> \param r13a ... 495!> \param sa ... 496!> \param e_0 ... 497!> \param npoints ... 498! ************************************************************************************************** 499 SUBROUTINE tfw_p_0(rhoa, r13a, sa, e_0, npoints) 500 501 REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rhoa, r13a, sa 502 REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_0 503 INTEGER, INTENT(in) :: npoints 504 505 INTEGER :: ip 506 507!$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)& 508!$OMP SHARED(npoints, rhoa,eps_rho,e_0,r13a,sa,flsd,fvw) 509 DO ip = 1, npoints 510 511 IF (rhoa(ip) > eps_rho) THEN 512 e_0(ip) = e_0(ip) + flsd*r13a(ip)*r13a(ip)*rhoa(ip) + fvw*sa(ip) 513 END IF 514 515 END DO 516 517 END SUBROUTINE tfw_p_0 518 519! ************************************************************************************************** 520!> \brief ... 521!> \param rhoa ... 522!> \param grhoa ... 523!> \param r13a ... 524!> \param sa ... 525!> \param e_rho ... 526!> \param e_ndrho ... 527!> \param npoints ... 528! ************************************************************************************************** 529 SUBROUTINE tfw_p_1(rhoa, grhoa, r13a, sa, e_rho, e_ndrho, npoints) 530 531 REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rhoa, grhoa, r13a, sa 532 REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_rho, e_ndrho 533 INTEGER, INTENT(in) :: npoints 534 535 INTEGER :: ip 536 REAL(KIND=dp) :: f 537 538 f = f53*flsd 539 540!$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)& 541!$OMP SHARED(npoints,rhoa,eps_rho,r13a,sa,fvw,grhoa,e_rho,e_ndrho,f) 542 DO ip = 1, npoints 543 544 IF (rhoa(ip) > eps_rho) THEN 545 e_rho(ip) = e_rho(ip) + f*r13a(ip)*r13a(ip) - fvw*sa(ip)/rhoa(ip) 546 e_ndrho(ip) = e_ndrho(ip) + 2.0_dp*fvw*grhoa(ip)/rhoa(ip) 547 END IF 548 549 END DO 550 551 END SUBROUTINE tfw_p_1 552 553! ************************************************************************************************** 554!> \brief ... 555!> \param rhoa ... 556!> \param grhoa ... 557!> \param r13a ... 558!> \param sa ... 559!> \param e_rho_rho ... 560!> \param e_rho_ndrho ... 561!> \param e_ndrho_ndrho ... 562!> \param npoints ... 563! ************************************************************************************************** 564 SUBROUTINE tfw_p_2(rhoa, grhoa, r13a, sa, e_rho_rho, e_rho_ndrho, & 565 e_ndrho_ndrho, npoints) 566 567 REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rhoa, grhoa, r13a, sa 568 REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_rho_rho, e_rho_ndrho, e_ndrho_ndrho 569 INTEGER, INTENT(in) :: npoints 570 571 INTEGER :: ip 572 REAL(KIND=dp) :: f 573 574 f = f23*f53*flsd 575 576!$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)& 577!$OMP SHARED(npoints,rhoa,eps_rho,e_rho_rho,f,fvw,r13a,sa,e_rho_ndrho,e_ndrho_ndrho) 578 DO ip = 1, npoints 579 580 IF (rhoa(ip) > eps_rho) THEN 581 e_rho_rho(ip) = e_rho_rho(ip) & 582 + f/r13a(ip) + 2.0_dp*fvw*sa(ip)/(rhoa(ip)*rhoa(ip)) 583 e_rho_ndrho(ip) = e_rho_ndrho(ip) & 584 - 2.0_dp*fvw*grhoa(ip)/(rhoa(ip)*rhoa(ip)) 585 e_ndrho_ndrho(ip) = e_ndrho_ndrho(ip) + 2.0_dp*fvw/rhoa(ip) 586 END IF 587 588 END DO 589 590 END SUBROUTINE tfw_p_2 591 592! ************************************************************************************************** 593!> \brief ... 594!> \param rhoa ... 595!> \param grhoa ... 596!> \param r13a ... 597!> \param sa ... 598!> \param e_rho_rho_rho ... 599!> \param e_rho_rho_ndrho ... 600!> \param e_rho_ndrho_ndrho ... 601!> \param npoints ... 602! ************************************************************************************************** 603 SUBROUTINE tfw_p_3(rhoa, grhoa, r13a, sa, e_rho_rho_rho, e_rho_rho_ndrho, & 604 e_rho_ndrho_ndrho, npoints) 605 606 REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rhoa, grhoa, r13a, sa 607 REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_rho_rho_rho, e_rho_rho_ndrho, & 608 e_rho_ndrho_ndrho 609 INTEGER, INTENT(in) :: npoints 610 611 INTEGER :: ip 612 REAL(KIND=dp) :: f 613 614 f = -f13*f23*f53*flsd 615 616!$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)& 617!$OMP SHARED(npoints,rhoa,eps_rho,e_rho_rho_rho,e_rho_rho_ndrho,e_rho_ndrho_ndrho,f,fvw,sa,grhoa) 618 DO ip = 1, npoints 619 620 IF (rhoa(ip) > eps_rho) THEN 621 e_rho_rho_rho(ip) = e_rho_rho_rho(ip) & 622 + f/(r13a(ip)*rhoa(ip)) & 623 - 6.0_dp*fvw*sa(ip)/(rhoa(ip)*rhoa(ip)*rhoa(ip)) 624 e_rho_rho_ndrho(ip) = e_rho_rho_ndrho(ip) & 625 + 4.0_dp*fvw*grhoa(ip)/(rhoa(ip)*rhoa(ip)*rhoa(ip)) 626 e_rho_ndrho_ndrho(ip) = e_rho_ndrho_ndrho(ip) & 627 - 2.0_dp*fvw/(rhoa(ip)*rhoa(ip)) 628 END IF 629 630 END DO 631 632 END SUBROUTINE tfw_p_3 633 634END MODULE xc_tfw 635 636