1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Calculation of non local dispersion functionals 8!> Some routines adapted from: 9!> Copyright (C) 2001-2009 Quantum ESPRESSO group 10!> Copyright (C) 2009 Brian Kolb, Timo Thonhauser - Wake Forest University 11!> This file is distributed under the terms of the 12!> GNU General Public License. See the file `License' 13!> in the root directory of the present distribution, 14!> or http://www.gnu.org/copyleft/gpl.txt . 15!> \author JGH 16! ************************************************************************************************** 17MODULE qs_dispersion_nonloc 18 USE bibliography, ONLY: Dion2004,& 19 Romanperez2009,& 20 Sabatini2013,& 21 cite_reference 22 USE cp_files, ONLY: close_file,& 23 open_file 24 USE cp_para_types, ONLY: cp_para_env_type 25 USE input_constants, ONLY: vdw_nl_DRSLL,& 26 vdw_nl_LMKLL,& 27 vdw_nl_RVV10,& 28 xc_vdw_fun_nonloc 29 USE kinds, ONLY: default_string_length,& 30 dp 31 USE mathconstants, ONLY: pi 32 USE message_passing, ONLY: mp_bcast,& 33 mp_sum 34 USE pw_grid_types, ONLY: HALFSPACE,& 35 pw_grid_type 36 USE pw_methods, ONLY: pw_axpy,& 37 pw_derive,& 38 pw_transfer 39 USE pw_pool_types, ONLY: pw_pool_create_pw,& 40 pw_pool_give_back_pw,& 41 pw_pool_type 42 USE pw_types, ONLY: COMPLEXDATA1D,& 43 REALDATA3D,& 44 REALSPACE,& 45 RECIPROCALSPACE,& 46 pw_p_type,& 47 pw_type 48 USE qs_dispersion_types, ONLY: qs_dispersion_type 49 USE virial_types, ONLY: virial_type 50#include "./base/base_uses.f90" 51 52 IMPLICIT NONE 53 54 PRIVATE 55 56 REAL(KIND=dp), PARAMETER :: epsr = 1.e-12_dp 57 58 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_dispersion_nonloc' 59 60 PUBLIC :: qs_dispersion_nonloc_init, calculate_dispersion_nonloc 61 62! ************************************************************************************************** 63 64CONTAINS 65 66! ************************************************************************************************** 67!> \brief ... 68!> \param dispersion_env ... 69!> \param para_env ... 70! ************************************************************************************************** 71 SUBROUTINE qs_dispersion_nonloc_init(dispersion_env, para_env) 72 TYPE(qs_dispersion_type), POINTER :: dispersion_env 73 TYPE(cp_para_env_type), POINTER :: para_env 74 75 CHARACTER(len=*), PARAMETER :: routineN = 'qs_dispersion_nonloc_init', & 76 routineP = moduleN//':'//routineN 77 78 CHARACTER(LEN=default_string_length) :: filename 79 INTEGER :: funit, handle, nqs, nr_points, q1_i, & 80 q2_i, vdw_type 81 82 CALL timeset(routineN, handle) 83 84 SELECT CASE (dispersion_env%nl_type) 85 CASE DEFAULT 86 CPABORT("Unknown vdW-DF functional") 87 CASE (vdw_nl_DRSLL, vdw_nl_LMKLL) 88 CALL cite_reference(Dion2004) 89 CASE (vdw_nl_RVV10) 90 CALL cite_reference(Sabatini2013) 91 END SELECT 92 CALL cite_reference(RomanPerez2009) 93 94 vdw_type = dispersion_env%type 95 SELECT CASE (vdw_type) 96 CASE DEFAULT 97 ! do nothing 98 CASE (xc_vdw_fun_nonloc) 99 ! setup information on non local functionals 100 filename = dispersion_env%kernel_file_name 101 IF (para_env%ionode) THEN 102 ! Read the kernel information from file "filename" 103 CALL open_file(file_name=filename, unit_number=funit, file_form="FORMATTED") 104 READ (funit, *) nqs, nr_points 105 READ (funit, *) dispersion_env%r_max 106 END IF 107 CALL mp_bcast(nqs, para_env%source, para_env%group) 108 CALL mp_bcast(nr_points, para_env%source, para_env%group) 109 CALL mp_bcast(dispersion_env%r_max, para_env%source, para_env%group) 110 ALLOCATE (dispersion_env%q_mesh(nqs), dispersion_env%kernel(0:nr_points, nqs, nqs), & 111 dispersion_env%d2phi_dk2(0:nr_points, nqs, nqs)) 112 dispersion_env%nqs = nqs 113 dispersion_env%nr_points = nr_points 114 IF (para_env%ionode) THEN 115 !! Read in the values of the q points used to generate this kernel 116 READ (funit, "(1p, 4e23.14)") dispersion_env%q_mesh 117 !! For each pair of q values, read in the function phi_q1_q2(k). 118 !! That is, the fourier transformed kernel function assuming q1 and q2 119 !! for all the values of r used. 120 DO q1_i = 1, nqs 121 DO q2_i = 1, q1_i 122 READ (funit, "(1p, 4e23.14)") dispersion_env%kernel(0:nr_points, q1_i, q2_i) 123 dispersion_env%kernel(0:nr_points, q2_i, q1_i) = dispersion_env%kernel(0:nr_points, q1_i, q2_i) 124 END DO 125 END DO 126 !! Again, for each pair of q values (q1 and q2), read in the value 127 !! of the second derivative of the above mentiond Fourier transformed 128 !! kernel function phi_alpha_beta(k). These are used for spline 129 !! interpolation of the Fourier transformed kernel. 130 DO q1_i = 1, nqs 131 DO q2_i = 1, q1_i 132 READ (funit, "(1p, 4e23.14)") dispersion_env%d2phi_dk2(0:nr_points, q1_i, q2_i) 133 dispersion_env%d2phi_dk2(0:nr_points, q2_i, q1_i) = dispersion_env%d2phi_dk2(0:nr_points, q1_i, q2_i) 134 END DO 135 END DO 136 CALL close_file(unit_number=funit) 137 END IF 138 CALL mp_bcast(dispersion_env%q_mesh, para_env%source, para_env%group) 139 CALL mp_bcast(dispersion_env%kernel, para_env%source, para_env%group) 140 CALL mp_bcast(dispersion_env%d2phi_dk2, para_env%source, para_env%group) 141 ! 2nd derivates for interpolation 142 ALLOCATE (dispersion_env%d2y_dx2(nqs, nqs)) 143 CALL initialize_spline_interpolation(dispersion_env%q_mesh, dispersion_env%d2y_dx2) 144 ! 145 dispersion_env%q_cut = dispersion_env%q_mesh(nqs) 146 dispersion_env%q_min = dispersion_env%q_mesh(1) 147 dispersion_env%dk = 2.0_dp*pi/dispersion_env%r_max 148 149 END SELECT 150 151 CALL timestop(handle) 152 153 END SUBROUTINE qs_dispersion_nonloc_init 154 155! ************************************************************************************************** 156!> \brief Calculates the non-local vdW functional using the method of Soler 157!> For spin polarized cases we use E(a,b) = E(a+b), i.e. total density 158!> \param vxc_rho ... 159!> \param rho_r ... 160!> \param rho_g ... 161!> \param edispersion ... 162!> \param dispersion_env ... 163!> \param energy_only ... 164!> \param pw_pool ... 165!> \param xc_pw_pool ... 166!> \param para_env ... 167!> \param virial ... 168! ************************************************************************************************** 169 SUBROUTINE calculate_dispersion_nonloc(vxc_rho, rho_r, rho_g, edispersion, & 170 dispersion_env, energy_only, pw_pool, xc_pw_pool, para_env, virial) 171 TYPE(pw_p_type), DIMENSION(:), POINTER :: vxc_rho, rho_r, rho_g 172 REAL(KIND=dp), INTENT(OUT) :: edispersion 173 TYPE(qs_dispersion_type), POINTER :: dispersion_env 174 LOGICAL, INTENT(IN) :: energy_only 175 TYPE(pw_pool_type), POINTER :: pw_pool, xc_pw_pool 176 TYPE(cp_para_env_type), POINTER :: para_env 177 TYPE(virial_type), OPTIONAL, POINTER :: virial 178 179 CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_dispersion_nonloc', & 180 routineP = moduleN//':'//routineN 181 INTEGER, DIMENSION(3, 3), PARAMETER :: nd = RESHAPE((/1, 0, 0, 0, 1, 0, 0, 0, 1/), (/3, 3/)) 182 183 INTEGER :: handle, i, i_grid, idir, ispin, nl_type, & 184 np, nspin, p, q, r, s 185 INTEGER, DIMENSION(1:3) :: hi, lo, n 186 LOGICAL :: use_virial 187 REAL(KIND=dp) :: b_value, beta, const, Ec_nl, sumnp 188 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: dq0_dgradrho, dq0_drho, hpot, potential, & 189 q0, rho 190 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: drho, thetas 191 TYPE(pw_grid_type), POINTER :: grid 192 TYPE(pw_p_type), ALLOCATABLE, DIMENSION(:) :: thetas_g 193 TYPE(pw_p_type), ALLOCATABLE, DIMENSION(:, :) :: drho_r 194 TYPE(pw_type), POINTER :: tmp_g, tmp_r, vxc_g, vxc_r 195 196 CALL timeset(routineN, handle) 197 198 CPASSERT(ASSOCIATED(rho_r)) 199 CPASSERT(ASSOCIATED(rho_g)) 200 CPASSERT(ASSOCIATED(pw_pool)) 201 202 IF (PRESENT(virial)) THEN 203 use_virial = virial%pv_calculate .AND. (.NOT. virial%pv_numer) 204 ELSE 205 use_virial = .FALSE. 206 ENDIF 207 IF (use_virial) THEN 208 CPASSERT(.NOT. energy_only) 209 END IF 210 IF (.NOT. energy_only) THEN 211 CPASSERT(ASSOCIATED(vxc_rho)) 212 END IF 213 214 nl_type = dispersion_env%nl_type 215 216 b_value = dispersion_env%b_value 217 beta = 0.03125_dp*(3.0_dp/(b_value**2.0_dp))**0.75_dp 218 nspin = SIZE(rho_r) 219 220 const = 1.0_dp/(3.0_dp*SQRT(pi)*b_value**1.5_dp)/(pi**0.75_dp) 221 222 ! temporary arrays for FFT 223 CALL pw_pool_create_pw(pw_pool, tmp_g, use_data=COMPLEXDATA1D, in_space=RECIPROCALSPACE) 224 CALL pw_pool_create_pw(pw_pool, tmp_r, use_data=REALDATA3D, in_space=REALSPACE) 225 226 ! get density derivatives 227 ALLOCATE (drho_r(3, nspin)) 228 DO ispin = 1, nspin 229 DO idir = 1, 3 230 NULLIFY (drho_r(idir, ispin)%pw) 231 CALL pw_pool_create_pw(pw_pool, drho_r(idir, ispin)%pw, & 232 use_data=REALDATA3D, in_space=REALSPACE) 233 CALL pw_transfer(rho_g(ispin)%pw, tmp_g) 234 CALL pw_derive(tmp_g, nd(:, idir)) 235 CALL pw_transfer(tmp_g, drho_r(idir, ispin)%pw) 236 END DO 237 END DO 238 239 np = SIZE(tmp_r%cr3d) 240 ALLOCATE (rho(np), drho(np, 3)) !in the following loops, rho and drho _will_ have the same bounds 241 DO i = 1, 3 242 lo(i) = LBOUND(tmp_r%cr3d, i) 243 hi(i) = UBOUND(tmp_r%cr3d, i) 244 n(i) = hi(i) - lo(i) + 1 245 END DO 246!$OMP PARALLEL DO DEFAULT(NONE) & 247!$OMP SHARED(n,rho) & 248!$OMP COLLAPSE(3) 249 DO r = 0, n(3) - 1 250 DO q = 0, n(2) - 1 251 DO p = 0, n(1) - 1 252 rho(r*n(2)*n(1) + q*n(1) + p + 1) = 0.0_dp 253 END DO 254 END DO 255 END DO 256!$OMP END PARALLEL DO 257 DO i = 1, 3 258!$OMP PARALLEL DO DEFAULT(NONE) & 259!$OMP SHARED(n,i,drho) & 260!$OMP COLLAPSE(3) 261 DO r = 0, n(3) - 1 262 DO q = 0, n(2) - 1 263 DO p = 0, n(1) - 1 264 drho(r*n(2)*n(1) + q*n(1) + p + 1, i) = 0.0_dp 265 END DO 266 END DO 267 END DO 268!$OMP END PARALLEL DO 269 END DO 270 DO ispin = 1, nspin 271 CPASSERT(rho_r(ispin)%pw%in_use == REALDATA3D) 272 CALL pw_transfer(rho_g(ispin)%pw, tmp_g) 273 CALL pw_transfer(tmp_g, tmp_r) 274!$OMP PARALLEL DO DEFAULT(NONE) & 275!$OMP SHARED(n,lo,rho,tmp_r) & 276!$OMP PRIVATE(s) & 277!$OMP COLLAPSE(3) 278 DO r = 0, n(3) - 1 279 DO q = 0, n(2) - 1 280 DO p = 0, n(1) - 1 281 s = r*n(2)*n(1) + q*n(1) + p + 1 282 rho(s) = rho(s) + tmp_r%cr3d(p + lo(1), q + lo(2), r + lo(3)) 283 END DO 284 END DO 285 END DO 286!$OMP END PARALLEL DO 287 DO i = 1, 3 288!$OMP PARALLEL DO DEFAULT(NONE) & 289!$OMP SHARED(ispin,i,n,lo,drho,drho_r) & 290!$OMP PRIVATE(s) & 291!$OMP COLLAPSE(3) 292 DO r = 0, n(3) - 1 293 DO q = 0, n(2) - 1 294 DO p = 0, n(1) - 1 295 s = r*n(2)*n(1) + q*n(1) + p + 1 296 drho(s, i) = drho(s, i) + drho_r(i, ispin)%pw%cr3d(p + lo(1), q + lo(2), r + lo(3)) 297 END DO 298 END DO 299 END DO 300!$OMP END PARALLEL DO 301 END DO 302 END DO 303 !! --------------------------------------------------------------------------------- 304 !! Find the value of q0 for all assigned grid points. q is defined in equations 305 !! 11 and 12 of DION and q0 is the saturated version of q defined in equation 306 !! 5 of SOLER. This routine also returns the derivatives of the q0s with respect 307 !! to the charge-density and the gradient of the charge-density. These are needed 308 !! for the potential calculated below. 309 !! --------------------------------------------------------------------------------- 310 311 IF (energy_only) THEN 312 ALLOCATE (q0(np)) 313 SELECT CASE (nl_type) 314 CASE DEFAULT 315 CPABORT("Unknown vdW-DF functional") 316 CASE (vdw_nl_DRSLL, vdw_nl_LMKLL) 317 CALL get_q0_on_grid_eo_vdw(rho, drho, q0, dispersion_env) 318 CASE (vdw_nl_RVV10) 319 CALL get_q0_on_grid_eo_rvv10(rho, drho, q0, dispersion_env) 320 END SELECT 321 ELSE 322 ALLOCATE (q0(np), dq0_drho(np), dq0_dgradrho(np)) 323 SELECT CASE (nl_type) 324 CASE DEFAULT 325 CPABORT("Unknown vdW-DF functional") 326 CASE (vdw_nl_DRSLL, vdw_nl_LMKLL) 327 CALL get_q0_on_grid_vdw(rho, drho, q0, dq0_drho, dq0_dgradrho, dispersion_env) 328 CASE (vdw_nl_RVV10) 329 CALL get_q0_on_grid_rvv10(rho, drho, q0, dq0_drho, dq0_dgradrho, dispersion_env) 330 END SELECT 331 END IF 332 333 !! --------------------------------------------------------------------------------------------- 334 !! Here we allocate and calculate the theta functions appearing in equations 8-12 of SOLER. 335 !! They are defined as rho*P_i(q0(rho, gradient_rho)) for vdW and vdW2 or 336 !! constant*rho^(3/4)*P_i(q0(rho, gradient_rho)) for rVV10 where P_i is a polynomial that 337 !! interpolates a Kronecker delta function at the point q_i (taken from the q_mesh) and q0 is 338 !! the saturated version of q. 339 !! q is defined in equations 11 and 12 of DION and the saturation proceedure is defined in equation 5 340 !! of soler. This is the biggest memory consumer in the method since the thetas array is 341 !! (total # of FFT points)*Nqs complex numbers. In a parallel run, each processor will hold the 342 !! values of all the theta functions on just the points assigned to it. 343 !! -------------------------------------------------------------------------------------------------- 344 !! thetas are stored in reciprocal space as theta_i(k) because this is the way they are used later 345 !! for the convolution (equation 11 of SOLER). 346 !! -------------------------------------------------------------------------------------------------- 347 ALLOCATE (thetas(np, dispersion_env%nqs)) 348 !! Interpolate the P_i polynomials defined in equation 3 in SOLER for the particular 349 !! q0 values we have. 350 CALL spline_interpolation(dispersion_env%q_mesh, dispersion_env%d2y_dx2, q0, thetas) 351 352 !! Form the thetas where theta is defined as rho*p_i(q0) for vdW and vdW2 or 353 !! constant*rho^(3/4)*p_i(q0) for rVV10 354 !! ------------------------------------------------------------------------------------ 355 SELECT CASE (nl_type) 356 CASE DEFAULT 357 CPABORT("Unknown vdW-DF functional") 358 CASE (vdw_nl_DRSLL, vdw_nl_LMKLL) 359!$OMP PARALLEL DO DEFAULT( NONE ) & 360!$OMP SHARED( dispersion_env, thetas, rho) 361 DO i = 1, dispersion_env%nqs 362 thetas(:, i) = thetas(:, i)*rho(:) 363 END DO 364!$OMP END PARALLEL DO 365 CASE (vdw_nl_RVV10) 366!$OMP PARALLEL DO DEFAULT( NONE ) & 367!$OMP SHARED( np, rho, dispersion_env, thetas, const ) & 368!$OMP PRIVATE( i ) & 369!$OMP SCHEDULE(DYNAMIC) ! use dynamic to allow for possibility of cases having (rho(i_grid) .LE. epsr) 370 DO i_grid = 1, np 371 IF (rho(i_grid) > epsr) THEN 372 DO i = 1, dispersion_env%nqs 373 thetas(i_grid, i) = thetas(i_grid, i)*const*rho(i_grid)**(0.75_dp) 374 END DO 375 ELSE 376 thetas(i_grid, :) = 0.0_dp 377 END IF 378 END DO 379!$OMP END PARALLEL DO 380 END SELECT 381 !! ------------------------------------------------------------------------------------ 382 !! Get thetas in reciprocal space. 383 DO i = 1, 3 384 lo(i) = LBOUND(tmp_r%cr3d, i) 385 hi(i) = UBOUND(tmp_r%cr3d, i) 386 n(i) = hi(i) - lo(i) + 1 387 END DO 388 ALLOCATE (thetas_g(dispersion_env%nqs)) 389 DO i = 1, dispersion_env%nqs 390!$OMP PARALLEL DO DEFAULT(NONE) & 391!$OMP SHARED(i,n,lo,thetas,tmp_r) & 392!$OMP COLLAPSE(3) 393 DO r = 0, n(3) - 1 394 DO q = 0, n(2) - 1 395 DO p = 0, n(1) - 1 396 tmp_r%cr3d(p + lo(1), q + lo(2), r + lo(3)) = thetas(r*n(2)*n(1) + q*n(1) + p + 1, i) 397 END DO 398 END DO 399 END DO 400!$OMP END PARALLEL DO 401 NULLIFY (thetas_g(i)%pw) 402 CALL pw_pool_create_pw(pw_pool, thetas_g(i)%pw, use_data=COMPLEXDATA1D, in_space=RECIPROCALSPACE) 403 CALL pw_transfer(tmp_r, thetas_g(i)%pw) 404 END DO 405 grid => thetas_g(1)%pw%pw_grid 406 !! --------------------------------------------------------------------------------------------- 407 !! Carry out the integration in equation 8 of SOLER. This also turns the thetas array into the 408 !! precursor to the u_i(k) array which is inverse fourier transformed to get the u_i(r) functions 409 !! of SOLER equation 11. Add the energy we find to the output variable etxc. 410 !! -------------------------------------------------------------------------------------------------- 411 sumnp = np 412 CALL mp_sum(sumnp, para_env%group) 413 IF (use_virial) THEN 414 ! calculates kernel contribution to stress 415 CALL vdW_energy(thetas_g, dispersion_env, Ec_nl, energy_only, virial) 416 SELECT CASE (nl_type) 417 CASE (vdw_nl_RVV10) 418 Ec_nl = 0.5_dp*Ec_nl + beta*SUM(rho(:))*grid%vol/sumnp 419 END SELECT 420 ! calculates energy contribution to stress 421 ! potential contribution to stress is calculated together with other potentials (Hxc) 422 DO idir = 1, 3 423 virial%pv_xc(idir, idir) = virial%pv_xc(idir, idir) + Ec_nl 424 END DO 425 ELSE 426 CALL vdW_energy(thetas_g, dispersion_env, Ec_nl, energy_only) 427 SELECT CASE (nl_type) 428 CASE (vdw_nl_RVV10) 429 Ec_nl = 0.5_dp*Ec_nl + beta*SUM(rho(:))*grid%vol/sumnp 430 END SELECT 431 END IF 432 CALL mp_sum(Ec_nl, para_env%group) 433 edispersion = Ec_nl 434 435 IF (energy_only) THEN 436 DEALLOCATE (q0) 437 ELSE 438 !! ---------------------------------------------------------------------------- 439 !! Inverse Fourier transform the u_i(k) to get the u_i(r) of SOLER equation 11. 440 !!----------------------------------------------------------------------------- 441 DO i = 1, 3 442 lo(i) = LBOUND(tmp_r%cr3d, i) 443 hi(i) = UBOUND(tmp_r%cr3d, i) 444 n(i) = hi(i) - lo(i) + 1 445 END DO 446 DO i = 1, dispersion_env%nqs 447 CALL pw_transfer(thetas_g(i)%pw, tmp_r) 448!$OMP PARALLEL DO DEFAULT(NONE) & 449!$OMP SHARED(i,n,lo,thetas,tmp_r) & 450!$OMP COLLAPSE(3) 451 DO r = 0, n(3) - 1 452 DO q = 0, n(2) - 1 453 DO p = 0, n(1) - 1 454 thetas(r*n(2)*n(1) + q*n(1) + p + 1, i) = tmp_r%cr3d(p + lo(1), q + lo(2), r + lo(3)) 455 END DO 456 END DO 457 END DO 458!$OMP END PARALLEL DO 459 END DO 460 !! ------------------------------------------------------------------------- 461 !! Here we allocate the array to hold the potential. This is calculated via 462 !! equation 10 of SOLER, using the u_i(r) calculated from equations 11 and 463 !! 12 of SOLER. Each processor allocates the array to be the size of the 464 !! full grid because, as can be seen in SOLER equation 9, processors need 465 !! to access grid points outside their allocated regions. 466 !! ------------------------------------------------------------------------- 467 ALLOCATE (potential(np), hpot(np)) 468 IF (use_virial) THEN 469 ! calculates gradient contribution to stress 470 grid => tmp_g%pw_grid 471 CALL get_potential(q0, dq0_drho, dq0_dgradrho, rho, thetas, potential, hpot, & 472 dispersion_env, drho, grid%dvol, virial) 473 ELSE 474 CALL get_potential(q0, dq0_drho, dq0_dgradrho, rho, thetas, potential, hpot, & 475 dispersion_env) 476 END IF 477 SELECT CASE (nl_type) 478 CASE (vdw_nl_RVV10) 479 potential(:) = 0.5_dp*potential(:) + beta 480 hpot(:) = 0.5_dp*hpot(:) 481 END SELECT 482 CALL pw_pool_create_pw(pw_pool, vxc_r, use_data=REALDATA3D, in_space=REALSPACE) 483 DO i = 1, 3 484 lo(i) = LBOUND(vxc_r%cr3d, i) 485 hi(i) = UBOUND(vxc_r%cr3d, i) 486 n(i) = hi(i) - lo(i) + 1 487 END DO 488!$OMP PARALLEL DO DEFAULT(NONE) & 489!$OMP SHARED(i,n,lo,potential,vxc_r) & 490!$OMP COLLAPSE(3) 491 DO r = 0, n(3) - 1 492 DO q = 0, n(2) - 1 493 DO p = 0, n(1) - 1 494 vxc_r%cr3d(p + lo(1), q + lo(2), r + lo(3)) = potential(r*n(2)*n(1) + q*n(1) + p + 1) 495 END DO 496 END DO 497 END DO 498!$OMP END PARALLEL DO 499 DO i = 1, 3 500 lo(i) = LBOUND(tmp_r%cr3d, i) 501 hi(i) = UBOUND(tmp_r%cr3d, i) 502 n(i) = hi(i) - lo(i) + 1 503 END DO 504 DO idir = 1, 3 505!$OMP PARALLEL DO DEFAULT(NONE) & 506!$OMP SHARED(n,lo,tmp_r) & 507!$OMP COLLAPSE(3) 508 DO r = 0, n(3) - 1 509 DO q = 0, n(2) - 1 510 DO p = 0, n(1) - 1 511 tmp_r%cr3d(p + lo(1), q + lo(2), r + lo(3)) = 0.0_dp 512 END DO 513 END DO 514 END DO 515!$OMP END PARALLEL DO 516 DO ispin = 1, nspin 517!$OMP PARALLEL DO DEFAULT(NONE) & 518!$OMP SHARED(n,lo,tmp_r,hpot,drho_r,idir,ispin) & 519!$OMP COLLAPSE(3) 520 DO r = 0, n(3) - 1 521 DO q = 0, n(2) - 1 522 DO p = 0, n(1) - 1 523 tmp_r%cr3d(p + lo(1), q + lo(2), r + lo(3)) = tmp_r%cr3d(p + lo(1), q + lo(2), r + lo(3)) & 524 + hpot(r*n(2)*n(1) + q*n(1) + p + 1) & 525 *drho_r(idir, ispin)%pw%cr3d(p + lo(1), q + lo(2), r + lo(3)) 526 END DO 527 END DO 528 END DO 529!$OMP END PARALLEL DO 530 END DO 531 CALL pw_transfer(tmp_r, tmp_g) 532 CALL pw_derive(tmp_g, nd(:, idir)) 533 CALL pw_transfer(tmp_g, tmp_r) 534 CALL pw_axpy(tmp_r, vxc_r, -1._dp) 535 END DO 536 CALL pw_transfer(vxc_r, tmp_g) 537 CALL pw_pool_give_back_pw(pw_pool, vxc_r) 538 CALL pw_pool_create_pw(xc_pw_pool, vxc_r, use_data=REALDATA3D, in_space=REALSPACE) 539 CALL pw_pool_create_pw(xc_pw_pool, vxc_g, use_data=COMPLEXDATA1D, in_space=RECIPROCALSPACE) 540 CALL pw_transfer(tmp_g, vxc_g) 541 CALL pw_transfer(vxc_g, vxc_r) 542 DO ispin = 1, nspin 543 CALL pw_axpy(vxc_r, vxc_rho(ispin)%pw, 1._dp) 544 END DO 545 CALL pw_pool_give_back_pw(xc_pw_pool, vxc_r) 546 CALL pw_pool_give_back_pw(xc_pw_pool, vxc_g) 547 DEALLOCATE (q0, dq0_drho, dq0_dgradrho) 548 END IF 549 550 DEALLOCATE (thetas) 551 552 DO i = 1, dispersion_env%nqs 553 CALL pw_pool_give_back_pw(pw_pool, thetas_g(i)%pw) 554 END DO 555 DO ispin = 1, nspin 556 DO idir = 1, 3 557 CALL pw_pool_give_back_pw(pw_pool, drho_r(idir, ispin)%pw) 558 END DO 559 END DO 560 CALL pw_pool_give_back_pw(pw_pool, tmp_r) 561 CALL pw_pool_give_back_pw(pw_pool, tmp_g) 562 563 DEALLOCATE (rho, drho, drho_r, thetas_g) 564 565 CALL timestop(handle) 566 567 END SUBROUTINE calculate_dispersion_nonloc 568 569! ************************************************************************************************** 570!> \brief This routine carries out the integration of equation 8 of SOLER. It returns the non-local 571!> exchange-correlation energy and the u_alpha(k) arrays used to find the u_alpha(r) arrays via 572!> equations 11 and 12 in SOLER. 573!> energy contribution to stress is added in qs_force 574!> \param thetas_g ... 575!> \param dispersion_env ... 576!> \param vdW_xc_energy ... 577!> \param energy_only ... 578!> \param virial ... 579!> \par History 580!> OpenMP added: Aug 2016 MTucker 581! ************************************************************************************************** 582 SUBROUTINE vdW_energy(thetas_g, dispersion_env, vdW_xc_energy, energy_only, virial) 583 TYPE(pw_p_type), ALLOCATABLE, DIMENSION(:) :: thetas_g 584 TYPE(qs_dispersion_type), POINTER :: dispersion_env 585 REAL(KIND=dp), INTENT(OUT) :: vdW_xc_energy 586 LOGICAL, INTENT(IN) :: energy_only 587 TYPE(virial_type), OPTIONAL, POINTER :: virial 588 589 CHARACTER(LEN=*), PARAMETER :: routineN = 'vdW_energy', routineP = moduleN//':'//routineN 590 591 COMPLEX(KIND=dp) :: uu 592 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: theta 593 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: u_vdw(:, :) 594 INTEGER :: handle, ig, iq, l, m, nl_type, nqs, & 595 q1_i, q2_i 596 LOGICAL :: use_virial 597 REAL(KIND=dp) :: g, g2, g2_last, g_multiplier, gm 598 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dkernel_of_dk, kernel_of_k 599 REAL(KIND=dp), DIMENSION(3, 3) :: virial_thread 600 TYPE(pw_grid_type), POINTER :: grid 601 602 CALL timeset(routineN, handle) 603 nqs = dispersion_env%nqs 604 605 IF (PRESENT(virial)) THEN 606 use_virial = .TRUE. 607 virial_thread(:, :) = 0.0_dp 608 ELSE 609 use_virial = .FALSE. 610 END IF 611 612 vdW_xc_energy = 0._dp 613 grid => thetas_g(1)%pw%pw_grid 614 615 IF (grid%grid_span == HALFSPACE) THEN 616 g_multiplier = 2._dp 617 ELSE 618 g_multiplier = 1._dp 619 END IF 620 621 nl_type = dispersion_env%nl_type 622 623 IF (.NOT. energy_only) THEN 624 ALLOCATE (u_vdW(grid%ngpts_cut_local, nqs)) 625 u_vdW(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp) 626 END IF 627 628!$OMP PARALLEL DEFAULT( NONE ) & 629!$OMP SHARED( nqs, energy_only, grid, dispersion_env & 630!$OMP , use_virial, thetas_g, g_multiplier, nl_type & 631!$OMP , u_vdW & 632!$OMP ) & 633!$OMP PRIVATE( g2_last, kernel_of_k, dkernel_of_dk, theta & 634!$OMP , g2, g, iq & 635!$OMP , q2_i, uu, q1_i, gm, l, m & 636!$OMP ) & 637!$OMP REDUCTION( +: vdW_xc_energy, virial_thread & 638!$OMP ) 639 640 g2_last = HUGE(0._dp) 641 642 ALLOCATE (kernel_of_k(nqs, nqs), dkernel_of_dk(nqs, nqs)) 643 ALLOCATE (theta(nqs)) 644 645!$OMP DO 646 DO ig = 1, grid%ngpts_cut_local 647 g2 = grid%gsq(ig) 648 IF (ABS(g2 - g2_last) > 1.e-10) THEN 649 g2_last = g2 650 g = SQRT(g2) 651 CALL interpolate_kernel(g, kernel_of_k, dispersion_env) 652 IF (use_virial) CALL interpolate_dkernel_dk(g, dkernel_of_dk, dispersion_env) 653 END IF 654 DO iq = 1, nqs 655 theta(iq) = thetas_g(iq)%pw%cc(ig) 656 END DO 657 DO q2_i = 1, nqs 658 uu = CMPLX(0.0_dp, 0.0_dp, KIND=dp) 659 DO q1_i = 1, nqs 660 uu = uu + kernel_of_k(q2_i, q1_i)*theta(q1_i) 661 END DO 662 IF (ig < grid%first_gne0) THEN 663 vdW_xc_energy = vdW_xc_energy + REAL(uu*CONJG(theta(q2_i)), KIND=dp) 664 ELSE 665 vdW_xc_energy = vdW_xc_energy + g_multiplier*REAL(uu*CONJG(theta(q2_i)), KIND=dp) 666 END IF 667 IF (.NOT. energy_only) u_vdW(ig, q2_i) = uu 668 669 IF (use_virial .AND. ig >= grid%first_gne0) THEN 670 DO q1_i = 1, nqs 671 gm = 0.5_dp*g_multiplier*grid%vol*dkernel_of_dk(q1_i, q2_i)*REAL(theta(q1_i)*CONJG(theta(q2_i)), KIND=dp) 672 IF (nl_type == vdw_nl_RVV10) THEN 673 gm = 0.5_dp*gm 674 END IF 675 DO l = 1, 3 676 DO m = 1, l 677 virial_thread(l, m) = virial_thread(l, m) - gm*(grid%g(l, ig)*grid%g(m, ig))/g 678 END DO 679 END DO 680 END DO 681 END IF 682 683 END DO 684 END DO 685!$OMP END DO 686 687 DEALLOCATE (theta) 688 DEALLOCATE (kernel_of_k, dkernel_of_dk) 689 690 IF (.NOT. energy_only) THEN 691!$OMP DO 692 DO ig = 1, grid%ngpts_cut_local 693 DO iq = 1, nqs 694 thetas_g(iq)%pw%cc(ig) = u_vdW(ig, iq) 695 END DO 696 END DO 697!$OMP END DO 698 END IF 699 700!$OMP END PARALLEL 701 702 IF (.NOT. energy_only) THEN 703 DEALLOCATE (u_vdW) 704 END IF 705 706 vdW_xc_energy = vdW_xc_energy*grid%vol*0.5_dp 707 708 IF (use_virial) THEN 709 DO l = 1, 3 710 DO m = 1, (l - 1) 711 virial%pv_xc(l, m) = virial%pv_xc(l, m) + virial_thread(l, m) 712 virial%pv_xc(m, l) = virial%pv_xc(l, m) 713 END DO 714 m = l 715 virial%pv_xc(l, m) = virial%pv_xc(l, m) + virial_thread(l, m) 716 END DO 717 END IF 718 719 CALL timestop(handle) 720 721 END SUBROUTINE vdW_energy 722 723! ************************************************************************************************** 724!> \brief This routine finds the non-local correlation contribution to the potential 725!> (i.e. the derivative of the non-local piece of the energy with respect to 726!> density) given in SOLER equation 10. The u_alpha(k) functions were found 727!> while calculating the energy. They are passed in as the matrix u_vdW. 728!> Most of the required derivatives were calculated in the "get_q0_on_grid" 729!> routine, but the derivative of the interpolation polynomials, P_alpha(q), 730!> (SOLER equation 3) with respect to q is interpolated here, along with the 731!> polynomials themselves. 732!> \param q0 ... 733!> \param dq0_drho ... 734!> \param dq0_dgradrho ... 735!> \param total_rho ... 736!> \param u_vdW ... 737!> \param potential ... 738!> \param h_prefactor ... 739!> \param dispersion_env ... 740!> \param drho ... 741!> \param dvol ... 742!> \param virial ... 743!> \par History 744!> OpenMP added: Aug 2016 MTucker 745! ************************************************************************************************** 746 SUBROUTINE get_potential(q0, dq0_drho, dq0_dgradrho, total_rho, u_vdW, potential, h_prefactor, & 747 dispersion_env, drho, dvol, virial) 748 749 REAL(dp), DIMENSION(:), INTENT(in) :: q0, dq0_drho, dq0_dgradrho, total_rho 750 REAL(dp), DIMENSION(:, :), INTENT(in) :: u_vdW 751 REAL(dp), DIMENSION(:), INTENT(out) :: potential, h_prefactor 752 TYPE(qs_dispersion_type), POINTER :: dispersion_env 753 REAL(dp), DIMENSION(:, :), INTENT(in), OPTIONAL :: drho 754 REAL(dp), INTENT(IN), OPTIONAL :: dvol 755 TYPE(virial_type), OPTIONAL, POINTER :: virial 756 757 CHARACTER(len=*), PARAMETER :: routineN = 'get_potential', routineP = moduleN//':'//routineN 758 759 INTEGER :: handle, i_grid, l, m, nl_type, nqs, P_i, & 760 q, q_hi, q_low 761 LOGICAL :: use_virial 762 REAL(dp) :: a, b, b_value, c, const, d, dP_dq0, dq, & 763 dq_6, e, f, P, prefactor, tmp_1_2, & 764 tmp_1_4, tmp_3_4 765 REAL(dp), ALLOCATABLE, DIMENSION(:) :: y 766 REAL(dp), DIMENSION(3, 3) :: virial_thread 767 REAL(dp), DIMENSION(:), POINTER :: q_mesh 768 REAL(dp), DIMENSION(:, :), POINTER :: d2y_dx2 769 770 CALL timeset(routineN, handle) 771 772 IF (PRESENT(virial)) THEN 773 use_virial = .TRUE. 774 virial_thread(:, :) = 0.0_dp 775 CPASSERT(PRESENT(drho)) 776 CPASSERT(PRESENT(dvol)) 777 ELSE 778 use_virial = .FALSE. 779 END IF 780 781 b_value = dispersion_env%b_value 782 const = 1.0_dp/(3.0_dp*b_value**(3.0_dp/2.0_dp)*pi**(5.0_dp/4.0_dp)) 783 potential = 0.0_dp 784 h_prefactor = 0.0_dp 785 786 d2y_dx2 => dispersion_env%d2y_dx2 787 q_mesh => dispersion_env%q_mesh 788 nqs = dispersion_env%nqs 789 nl_type = dispersion_env%nl_type 790 791!$OMP PARALLEL DEFAULT( NONE ) & 792!$OMP SHARED( nqs, u_vdW, q_mesh, q0, d2y_dx2, nl_type & 793!$OMP , potential, h_prefactor & 794!$OMP , dq0_drho, dq0_dgradrho, total_rho, const & 795!$OMP , use_virial, drho, dvol, virial & 796!$OMP ) & 797!$OMP PRIVATE( y & 798!$OMP , q_low, q_hi, q, dq, dq_6, A, b, c, d, e, f & 799!$OMP , P_i, dP_dq0, P, prefactor, l, m & 800!$OMP , tmp_1_2, tmp_1_4, tmp_3_4 & 801!$OMP ) & 802!$OMP REDUCTION(+: virial_thread & 803!$OMP ) 804 805 ALLOCATE (y(nqs)) 806 807!$OMP DO 808 DO i_grid = 1, SIZE(u_vdW, 1) 809 q_low = 1 810 q_hi = nqs 811 ! Figure out which bin our value of q0 is in in the q_mesh 812 DO WHILE ((q_hi - q_low) > 1) 813 q = INT((q_hi + q_low)/2) 814 IF (q_mesh(q) > q0(i_grid)) THEN 815 q_hi = q 816 ELSE 817 q_low = q 818 END IF 819 END DO 820 IF (q_hi == q_low) CPABORT("get_potential: qhi == qlow") 821 822 dq = q_mesh(q_hi) - q_mesh(q_low) 823 dq_6 = dq/6.0_dp 824 825 a = (q_mesh(q_hi) - q0(i_grid))/dq 826 b = (q0(i_grid) - q_mesh(q_low))/dq 827 c = (a**3 - a)*dq*dq_6 828 d = (b**3 - b)*dq*dq_6 829 e = (3.0_dp*a**2 - 1.0_dp)*dq_6 830 f = (3.0_dp*b**2 - 1.0_dp)*dq_6 831 832 DO P_i = 1, nqs 833 y = 0.0_dp 834 y(P_i) = 1.0_dp 835 dP_dq0 = (y(q_hi) - y(q_low))/dq - e*d2y_dx2(P_i, q_low) + f*d2y_dx2(P_i, q_hi) 836 P = a*y(q_low) + b*y(q_hi) + c*d2y_dx2(P_i, q_low) + d*d2y_dx2(P_i, q_hi) 837 !! The first term in equation 13 of SOLER 838 SELECT CASE (nl_type) 839 CASE DEFAULT 840 CPABORT("Unknown vdW-DF functional") 841 CASE (vdw_nl_DRSLL, vdw_nl_LMKLL) 842 potential(i_grid) = potential(i_grid) + u_vdW(i_grid, P_i)*(P + dP_dq0*dq0_drho(i_grid)) 843 prefactor = u_vdW(i_grid, P_i)*dP_dq0*dq0_dgradrho(i_grid) 844 CASE (vdw_nl_RVV10) 845 IF (total_rho(i_grid) > epsr) THEN 846 tmp_1_2 = SQRT(total_rho(i_grid)) 847 tmp_1_4 = SQRT(tmp_1_2) ! == total_rho(i_grid)**(1.0_dp/4.0_dp) 848 tmp_3_4 = tmp_1_4*tmp_1_4*tmp_1_4 ! == total_rho(i_grid)**(3.0_dp/4.0_dp) 849 potential(i_grid) = potential(i_grid) & 850 + u_vdW(i_grid, P_i)*(const*0.75_dp/tmp_1_4*P & 851 + const*tmp_3_4*dP_dq0*dq0_drho(i_grid)) 852 prefactor = u_vdW(i_grid, P_i)*const*tmp_3_4*dP_dq0*dq0_dgradrho(i_grid) 853 ELSE 854 prefactor = 0.0_dp 855 ENDIF 856 END SELECT 857 IF (q0(i_grid) .NE. q_mesh(nqs)) THEN 858 h_prefactor(i_grid) = h_prefactor(i_grid) + prefactor 859 END IF 860 861 IF (use_virial .AND. ABS(prefactor) > 0.0_dp) THEN 862 SELECT CASE (nl_type) 863 CASE DEFAULT 864 ! do nothing 865 CASE (vdw_nl_RVV10) 866 prefactor = 0.5_dp*prefactor 867 END SELECT 868 prefactor = prefactor*dvol 869 DO l = 1, 3 870 DO m = 1, l 871 virial_thread(l, m) = virial_thread(l, m) - prefactor*drho(i_grid, l)*drho(i_grid, m) 872 ENDDO 873 ENDDO 874 END IF 875 876 END DO ! P_i = 1, nqs 877 END DO ! i_grid = 1, SIZE(u_vdW, 1) 878!$OMP END DO 879 880 DEALLOCATE (y) 881!$OMP END PARALLEL 882 883 IF (use_virial) THEN 884 DO l = 1, 3 885 DO m = 1, (l - 1) 886 virial%pv_xc(l, m) = virial%pv_xc(l, m) + virial_thread(l, m) 887 virial%pv_xc(m, l) = virial%pv_xc(l, m) 888 ENDDO 889 m = l 890 virial%pv_xc(l, m) = virial%pv_xc(l, m) + virial_thread(l, m) 891 ENDDO 892 END IF 893 894 CALL timestop(handle) 895 END SUBROUTINE get_potential 896 897! ************************************************************************************************** 898!> \brief calculates exponent = sum(from i=1 to hi, ((alpha)**i)/i) ) without <<< calling power >>> 899!> \param hi = upper index for sum 900!> \param alpha ... 901!> \param exponent = output value 902!> \par History 903!> Created: MTucker, Aug 2016 904! ************************************************************************************************** 905 SUBROUTINE calculate_exponent(hi, alpha, exponent) 906 INTEGER, INTENT(in) :: hi 907 REAL(dp), INTENT(in) :: alpha 908 REAL(dp), INTENT(out) :: exponent 909 910 INTEGER :: i 911 REAL(dp) :: multiplier 912 913 multiplier = alpha 914 exponent = alpha 915 916 DO i = 2, hi 917 multiplier = multiplier*alpha 918 exponent = exponent + (multiplier/i) 919 END DO 920 END SUBROUTINE calculate_exponent 921 922! ************************************************************************************************** 923!> \brief calculate exponent = sum(from i=1 to hi, ((alpha)**i)/i) ) without calling power 924!> also calculates derivative using similar series 925!> \param hi = upper index for sum 926!> \param alpha ... 927!> \param exponent = output value 928!> \param derivative ... 929!> \par History 930!> Created: MTucker, Aug 2016 931! ************************************************************************************************** 932 SUBROUTINE calculate_exponent_derivative(hi, alpha, exponent, derivative) 933 INTEGER, INTENT(in) :: hi 934 REAL(dp), INTENT(in) :: alpha 935 REAL(dp), INTENT(out) :: exponent, derivative 936 937 INTEGER :: i 938 REAL(dp) :: multiplier 939 940 derivative = 0.0d0 941 multiplier = 1.0d0 942 exponent = 0.0d0 943 944 DO i = 1, hi 945 derivative = derivative + multiplier 946 multiplier = multiplier*alpha 947 exponent = exponent + (multiplier/i) 948 END DO 949 END SUBROUTINE calculate_exponent_derivative 950 951 !! This routine first calculates the q value defined in (DION equations 11 and 12), then 952 !! saturates it according to (SOLER equation 5). 953! ************************************************************************************************** 954!> \brief This routine first calculates the q value defined in (DION equations 11 and 12), then 955!> saturates it according to (SOLER equation 5). 956!> \param total_rho ... 957!> \param gradient_rho ... 958!> \param q0 ... 959!> \param dq0_drho ... 960!> \param dq0_dgradrho ... 961!> \param dispersion_env ... 962! ************************************************************************************************** 963 SUBROUTINE get_q0_on_grid_vdw(total_rho, gradient_rho, q0, dq0_drho, dq0_dgradrho, dispersion_env) 964 !! 965 !! more specifically it calculates the following 966 !! 967 !! q0(ir) = q0 as defined above 968 !! dq0_drho(ir) = total_rho * d q0 /d rho 969 !! dq0_dgradrho = total_rho / |gradient_rho| * d q0 / d |gradient_rho| 970 !! 971 REAL(dp), INTENT(IN) :: total_rho(:), gradient_rho(:, :) 972 REAL(dp), INTENT(OUT) :: q0(:), dq0_drho(:), dq0_dgradrho(:) 973 TYPE(qs_dispersion_type), POINTER :: dispersion_env 974 975 CHARACTER(len=*), PARAMETER :: routineN = 'get_q0_on_grid_vdw', & 976 routineP = moduleN//':'//routineN 977 INTEGER, PARAMETER :: m_cut = 12 978 REAL(dp), PARAMETER :: LDA_A = 0.031091_dp, LDA_a1 = 0.2137_dp, LDA_b1 = 7.5957_dp, & 979 LDA_b2 = 3.5876_dp, LDA_b3 = 1.6382_dp, LDA_b4 = 0.49294_dp 980 981 INTEGER :: i_grid 982 REAL(dp) :: dq0_dq, exponent, gradient_correction, & 983 kF, LDA_1, LDA_2, q, q__q_cut, q_cut, & 984 q_min, r_s, sqrt_r_s, Z_ab 985 986 q_cut = dispersion_env%q_cut 987 q_min = dispersion_env%q_min 988 SELECT CASE (dispersion_env%nl_type) 989 CASE DEFAULT 990 CPABORT("Unknown vdW-DF functional") 991 CASE (vdw_nl_DRSLL) 992 Z_ab = -0.8491_dp 993 CASE (vdw_nl_LMKLL) 994 Z_ab = -1.887_dp 995 END SELECT 996 997 ! initialize q0-related arrays ... 998 q0(:) = q_cut 999 dq0_drho(:) = 0.0_dp 1000 dq0_dgradrho(:) = 0.0_dp 1001 1002 DO i_grid = 1, SIZE(total_rho) 1003 1004 !! This prevents numerical problems. If the charge density is negative (an 1005 !! unphysical situation), we simply treat it as very small. In that case, 1006 !! q0 will be very large and will be saturated. For a saturated q0 the derivative 1007 !! dq0_dq will be 0 so we set q0 = q_cut and dq0_drho = dq0_dgradrho = 0 and go on 1008 !! to the next point. 1009 !! ------------------------------------------------------------------------------------ 1010 IF (total_rho(i_grid) < epsr) CYCLE 1011 !! ------------------------------------------------------------------------------------ 1012 !! Calculate some intermediate values needed to find q 1013 !! ------------------------------------------------------------------------------------ 1014 kF = (3.0_dp*pi*pi*total_rho(i_grid))**(1.0_dp/3.0_dp) 1015 r_s = (3.0_dp/(4.0_dp*pi*total_rho(i_grid)))**(1.0_dp/3.0_dp) 1016 sqrt_r_s = SQRT(r_s) 1017 1018 gradient_correction = -Z_ab/(36.0_dp*kF*total_rho(i_grid)**2) & 1019 *(gradient_rho(i_grid, 1)**2 + gradient_rho(i_grid, 2)**2 + gradient_rho(i_grid, 3)**2) 1020 1021 LDA_1 = 8.0_dp*pi/3.0_dp*(LDA_A*(1.0_dp + LDA_a1*r_s)) 1022 LDA_2 = 2.0_dp*LDA_A*(LDA_b1*sqrt_r_s + LDA_b2*r_s + LDA_b3*r_s*sqrt_r_s + LDA_b4*r_s*r_s) 1023 !! --------------------------------------------------------------- 1024 !! This is the q value defined in equations 11 and 12 of DION 1025 !! --------------------------------------------------------------- 1026 q = kF + LDA_1*LOG(1.0_dp + 1.0_dp/LDA_2) + gradient_correction 1027 !! --------------------------------------------------------------- 1028 !! Here, we calculate q0 by saturating q according to equation 5 of SOLER. Also, we find 1029 !! the derivative dq0_dq needed for the derivatives dq0_drho and dq0_dgradrh0 discussed below. 1030 !! --------------------------------------------------------------------------------------- 1031 q__q_cut = q/q_cut 1032 CALL calculate_exponent_derivative(m_cut, q__q_cut, exponent, dq0_dq) 1033 q0(i_grid) = q_cut*(1.0_dp - EXP(-exponent)) 1034 dq0_dq = dq0_dq*EXP(-exponent) 1035 !! --------------------------------------------------------------------------------------- 1036 !! This is to handle a case with q0 too small. We simply set it to the smallest q value in 1037 !! out q_mesh. Hopefully this doesn't get used often (ever) 1038 !! --------------------------------------------------------------------------------------- 1039 IF (q0(i_grid) < q_min) THEN 1040 q0(i_grid) = q_min 1041 END IF 1042 !! --------------------------------------------------------------------------------------- 1043 !! Here we find derivatives. These are actually the density times the derivative of q0 with respect 1044 !! to rho and gradient_rho. The density factor comes in since we are really differentiating 1045 !! theta = (rho)*P(q0) with respect to density (or its gradient) which will be 1046 !! dtheta_drho = P(q0) + dP_dq0 * [rho * dq0_dq * dq_drho] and 1047 !! dtheta_dgradient_rho = dP_dq0 * [rho * dq0_dq * dq_dgradient_rho] 1048 !! The parts in square brackets are what is calculated here. The dP_dq0 term will be interpolated 1049 !! later. There should actually be a factor of the magnitude of the gradient in the gradient_rho derivative 1050 !! but that cancels out when we differentiate the magnitude of the gradient with respect to a particular 1051 !! component. 1052 !! ------------------------------------------------------------------------------------------------ 1053 1054 dq0_drho(i_grid) = dq0_dq*(kF/3.0_dp - 7.0_dp/3.0_dp*gradient_correction & 1055 - 8.0_dp*pi/9.0_dp*LDA_A*LDA_a1*r_s*LOG(1.0_dp + 1.0_dp/LDA_2) & 1056 + LDA_1/(LDA_2*(1.0_dp + LDA_2)) & 1057 *(2.0_dp*LDA_A*(LDA_b1/6.0_dp*sqrt_r_s + LDA_b2/3.0_dp*r_s + LDA_b3/2.0_dp*r_s*sqrt_r_s & 1058 + 2.0_dp*LDA_b4/3.0_dp*r_s**2))) 1059 1060 dq0_dgradrho(i_grid) = total_rho(i_grid)*dq0_dq*2.0_dp*(-Z_ab)/(36.0_dp*kF*total_rho(i_grid)**2) 1061 1062 END DO 1063 1064 END SUBROUTINE get_q0_on_grid_vdw 1065 1066! ************************************************************************************************** 1067!> \brief ... 1068!> \param total_rho ... 1069!> \param gradient_rho ... 1070!> \param q0 ... 1071!> \param dq0_drho ... 1072!> \param dq0_dgradrho ... 1073!> \param dispersion_env ... 1074! ************************************************************************************************** 1075 SUBROUTINE get_q0_on_grid_rvv10(total_rho, gradient_rho, q0, dq0_drho, dq0_dgradrho, dispersion_env) 1076 !! 1077 !! more specifically it calculates the following 1078 !! 1079 !! q0(ir) = q0 as defined above 1080 !! dq0_drho(ir) = total_rho * d q0 /d rho 1081 !! dq0_dgradrho = total_rho / |gradient_rho| * d q0 / d |gradient_rho| 1082 !! 1083 REAL(dp), INTENT(IN) :: total_rho(:), gradient_rho(:, :) 1084 REAL(dp), INTENT(OUT) :: q0(:), dq0_drho(:), dq0_dgradrho(:) 1085 TYPE(qs_dispersion_type), POINTER :: dispersion_env 1086 1087 INTEGER, PARAMETER :: m_cut = 12 1088 1089 INTEGER :: i_grid 1090 REAL(dp) :: b_value, C_value, dk_dn, dq0_dq, dw0_dn, & 1091 exponent, gmod2, k, mod_grad, q, & 1092 q__q_cut, q_cut, q_min, w0, wg2, wp2 1093 1094 q_cut = dispersion_env%q_cut 1095 q_min = dispersion_env%q_min 1096 b_value = dispersion_env%b_value 1097 C_value = dispersion_env%c_value 1098 1099 ! initialize q0-related arrays ... 1100 q0(:) = q_cut 1101 dq0_drho(:) = 0.0_dp 1102 dq0_dgradrho(:) = 0.0_dp 1103 1104 DO i_grid = 1, SIZE(total_rho) 1105 1106 gmod2 = gradient_rho(i_grid, 1)**2 + gradient_rho(i_grid, 2)**2 + gradient_rho(i_grid, 3)**2 1107 1108 !if (total_rho(i_grid) > epsr .and. gmod2 > epsr) cycle 1109 IF (total_rho(i_grid) > epsr) THEN 1110 1111 !! Calculate some intermediate values needed to find q 1112 !! ------------------------------------------------------------------------------------ 1113 mod_grad = SQRT(gmod2) 1114 1115 wp2 = 16.0_dp*pi*total_rho(i_grid) 1116 wg2 = 4_dp*C_value*(mod_grad/total_rho(i_grid))**4 1117 1118 k = b_value*3.0_dp*pi*((total_rho(i_grid)/(9.0_dp*pi))**(1.0_dp/6.0_dp)) 1119 w0 = SQRT(wg2 + wp2/3.0_dp) 1120 1121 q = w0/k 1122 1123 !! Here, we calculate q0 by saturating q according 1124 !! --------------------------------------------------------------------------------------- 1125 q__q_cut = q/q_cut 1126 CALL calculate_exponent_derivative(m_cut, q__q_cut, exponent, dq0_dq) 1127 q0(i_grid) = q_cut*(1.0_dp - EXP(-exponent)) 1128 dq0_dq = dq0_dq*EXP(-exponent) 1129 1130 !! --------------------------------------------------------------------------------------- 1131 IF (q0(i_grid) < q_min) THEN 1132 q0(i_grid) = q_min 1133 END IF 1134 1135 !!---------------------------------Final values--------------------------------- 1136 dw0_dn = 1.0_dp/(2.0_dp*w0)*(16.0_dp/3.0_dp*pi - 4.0_dp*wg2/total_rho(i_grid)) 1137 dk_dn = k/(6.0_dp*total_rho(i_grid)) 1138 1139 dq0_drho(i_grid) = dq0_dq*1.0_dp/(k**2.0)*(dw0_dn*k - dk_dn*w0) 1140 dq0_dgradrho(i_grid) = dq0_dq*1.0_dp/(2.0_dp*k*w0)*4.0_dp*wg2/gmod2 1141 ENDIF 1142 1143 END DO 1144 1145 END SUBROUTINE get_q0_on_grid_rvv10 1146 1147! ************************************************************************************************** 1148!> \brief ... 1149!> \param total_rho ... 1150!> \param gradient_rho ... 1151!> \param q0 ... 1152!> \param dispersion_env ... 1153! ************************************************************************************************** 1154 SUBROUTINE get_q0_on_grid_eo_vdw(total_rho, gradient_rho, q0, dispersion_env) 1155 1156 REAL(dp), INTENT(IN) :: total_rho(:), gradient_rho(:, :) 1157 REAL(dp), INTENT(OUT) :: q0(:) 1158 TYPE(qs_dispersion_type), POINTER :: dispersion_env 1159 1160 CHARACTER(len=*), PARAMETER :: routineN = 'get_q0_on_grid_eo_vdw', & 1161 routineP = moduleN//':'//routineN 1162 INTEGER, PARAMETER :: m_cut = 12 1163 REAL(dp), PARAMETER :: LDA_A = 0.031091_dp, LDA_a1 = 0.2137_dp, LDA_b1 = 7.5957_dp, & 1164 LDA_b2 = 3.5876_dp, LDA_b3 = 1.6382_dp, LDA_b4 = 0.49294_dp 1165 1166 INTEGER :: i_grid 1167 REAL(dp) :: exponent, gradient_correction, kF, & 1168 LDA_1, LDA_2, q, q__q_cut, q_cut, & 1169 q_min, r_s, sqrt_r_s, Z_ab 1170 1171 q_cut = dispersion_env%q_cut 1172 q_min = dispersion_env%q_min 1173 SELECT CASE (dispersion_env%nl_type) 1174 CASE DEFAULT 1175 CPABORT("Unknown vdW-DF functional") 1176 CASE (vdw_nl_DRSLL) 1177 Z_ab = -0.8491_dp 1178 CASE (vdw_nl_LMKLL) 1179 Z_ab = -1.887_dp 1180 END SELECT 1181 1182 ! initialize q0-related arrays ... 1183 q0(:) = q_cut 1184 DO i_grid = 1, SIZE(total_rho) 1185 !! This prevents numerical problems. If the charge density is negative (an 1186 !! unphysical situation), we simply treat it as very small. In that case, 1187 !! q0 will be very large and will be saturated. For a saturated q0 the derivative 1188 !! dq0_dq will be 0 so we set q0 = q_cut and dq0_drho = dq0_dgradrho = 0 and go on 1189 !! to the next point. 1190 !! ------------------------------------------------------------------------------------ 1191 IF (total_rho(i_grid) < epsr) CYCLE 1192 !! ------------------------------------------------------------------------------------ 1193 !! Calculate some intermediate values needed to find q 1194 !! ------------------------------------------------------------------------------------ 1195 kF = (3.0_dp*pi*pi*total_rho(i_grid))**(1.0_dp/3.0_dp) 1196 r_s = (3.0_dp/(4.0_dp*pi*total_rho(i_grid)))**(1.0_dp/3.0_dp) 1197 sqrt_r_s = SQRT(r_s) 1198 1199 gradient_correction = -Z_ab/(36.0_dp*kF*total_rho(i_grid)**2) & 1200 *(gradient_rho(i_grid, 1)**2 + gradient_rho(i_grid, 2)**2 + gradient_rho(i_grid, 3)**2) 1201 1202 LDA_1 = 8.0_dp*pi/3.0_dp*(LDA_A*(1.0_dp + LDA_a1*r_s)) 1203 LDA_2 = 2.0_dp*LDA_A*(LDA_b1*sqrt_r_s + LDA_b2*r_s + LDA_b3*r_s*sqrt_r_s + LDA_b4*r_s*r_s) 1204 !! ------------------------------------------------------------------------------------ 1205 !! This is the q value defined in equations 11 and 12 of DION 1206 !! --------------------------------------------------------------- 1207 q = kF + LDA_1*LOG(1.0_dp + 1.0_dp/LDA_2) + gradient_correction 1208 1209 !! --------------------------------------------------------------- 1210 !! Here, we calculate q0 by saturating q according to equation 5 of SOLER. Also, we find 1211 !! the derivative dq0_dq needed for the derivatives dq0_drho and dq0_dgradrh0 discussed below. 1212 !! --------------------------------------------------------------------------------------- 1213 q__q_cut = q/q_cut 1214 CALL calculate_exponent(m_cut, q__q_cut, exponent) 1215 q0(i_grid) = q_cut*(1.0_dp - EXP(-exponent)) 1216 1217 !! --------------------------------------------------------------------------------------- 1218 !! This is to handle a case with q0 too small. We simply set it to the smallest q value in 1219 !! out q_mesh. Hopefully this doesn't get used often (ever) 1220 !! --------------------------------------------------------------------------------------- 1221 IF (q0(i_grid) < q_min) THEN 1222 q0(i_grid) = q_min 1223 END IF 1224 END DO 1225 1226 END SUBROUTINE get_q0_on_grid_eo_vdw 1227 1228! ************************************************************************************************** 1229!> \brief ... 1230!> \param total_rho ... 1231!> \param gradient_rho ... 1232!> \param q0 ... 1233!> \param dispersion_env ... 1234! ************************************************************************************************** 1235 SUBROUTINE get_q0_on_grid_eo_rvv10(total_rho, gradient_rho, q0, dispersion_env) 1236 1237 REAL(dp), INTENT(IN) :: total_rho(:), gradient_rho(:, :) 1238 REAL(dp), INTENT(OUT) :: q0(:) 1239 TYPE(qs_dispersion_type), POINTER :: dispersion_env 1240 1241 INTEGER, PARAMETER :: m_cut = 12 1242 1243 INTEGER :: i_grid 1244 REAL(dp) :: b_value, C_value, exponent, gmod2, k, q, & 1245 q__q_cut, q_cut, q_min, w0, wg2, wp2 1246 1247 q_cut = dispersion_env%q_cut 1248 q_min = dispersion_env%q_min 1249 b_value = dispersion_env%b_value 1250 C_value = dispersion_env%c_value 1251 1252 ! initialize q0-related arrays ... 1253 q0(:) = q_cut 1254 1255 DO i_grid = 1, SIZE(total_rho) 1256 1257 gmod2 = gradient_rho(i_grid, 1)**2 + gradient_rho(i_grid, 2)**2 + gradient_rho(i_grid, 3)**2 1258 1259 !if (total_rho(i_grid) > epsr .and. gmod2 > epsr) cycle 1260 IF (total_rho(i_grid) > epsr) THEN 1261 1262 !! Calculate some intermediate values needed to find q 1263 !! ------------------------------------------------------------------------------------ 1264 wp2 = 16.0_dp*pi*total_rho(i_grid) 1265 wg2 = 4_dp*C_value*(gmod2*gmod2)/(total_rho(i_grid)**4) 1266 1267 k = b_value*3.0_dp*pi*((total_rho(i_grid)/(9.0_dp*pi))**(1.0_dp/6.0_dp)) 1268 w0 = SQRT(wg2 + wp2/3.0_dp) 1269 1270 q = w0/k 1271 1272 !! Here, we calculate q0 by saturating q according 1273 !! --------------------------------------------------------------------------------------- 1274 q__q_cut = q/q_cut 1275 CALL calculate_exponent(m_cut, q__q_cut, exponent) 1276 q0(i_grid) = q_cut*(1.0_dp - EXP(-exponent)) 1277 1278 IF (q0(i_grid) < q_min) THEN 1279 q0(i_grid) = q_min 1280 END IF 1281 1282 ENDIF 1283 1284 END DO 1285 1286 END SUBROUTINE get_q0_on_grid_eo_rvv10 1287 1288! ************************************************************************************************** 1289!> \brief This routine is modeled after an algorithm from "Numerical Recipes in C" by Cambridge University 1290!> press, page 97. It was adapted for Fortran, of course and for the problem at hand, in that it finds 1291!> the bin a particular x value is in and then loops over all the P_i functions so we only have to find 1292!> the bin once. 1293!> \param x ... 1294!> \param d2y_dx2 ... 1295!> \param evaluation_points ... 1296!> \param values ... 1297!> \par History 1298!> OpenMP added: Aug 2016 MTucker 1299! ************************************************************************************************** 1300 SUBROUTINE spline_interpolation(x, d2y_dx2, evaluation_points, values) 1301 REAL(dp), INTENT(in) :: x(:), d2y_dx2(:, :), evaluation_points(:) 1302 REAL(dp), INTENT(out) :: values(:, :) 1303 1304 INTEGER :: i_grid, index, lower_bound, & 1305 Ngrid_points, Nx, P_i, upper_bound 1306 REAL(dp) :: a, b, c, const, d, dx 1307 REAL(dp), ALLOCATABLE :: y(:) 1308 1309 Nx = SIZE(x) 1310 Ngrid_points = SIZE(evaluation_points) 1311 1312!$OMP PARALLEL DEFAULT( NONE ) & 1313!$OMP SHARED( x, d2y_dx2, evaluation_points, values & 1314!$OMP , Nx, Ngrid_points & 1315!$OMP ) & 1316!$OMP PRIVATE( y, lower_bound, upper_bound, index & 1317!$OMP , dx, const, A, b, c, d, P_i & 1318!$OMP ) 1319 1320 ALLOCATE (y(Nx)) 1321 1322!$OMP DO 1323 DO i_grid = 1, Ngrid_points 1324 lower_bound = 1 1325 upper_bound = Nx 1326 DO WHILE ((upper_bound - lower_bound) > 1) 1327 index = (upper_bound + lower_bound)/2 1328 IF (evaluation_points(i_grid) > x(index)) THEN 1329 lower_bound = index 1330 ELSE 1331 upper_bound = index 1332 END IF 1333 END DO 1334 1335 dx = x(upper_bound) - x(lower_bound) 1336 const = dx*dx/6.0_dp 1337 1338 a = (x(upper_bound) - evaluation_points(i_grid))/dx 1339 b = (evaluation_points(i_grid) - x(lower_bound))/dx 1340 c = (a**3 - a)*const 1341 d = (b**3 - b)*const 1342 1343 DO P_i = 1, Nx 1344 y = 0 1345 y(P_i) = 1 1346 values(i_grid, P_i) = a*y(lower_bound) + b*y(upper_bound) & 1347 + (c*d2y_dx2(P_i, lower_bound) + d*d2y_dx2(P_i, upper_bound)) 1348 END DO 1349 END DO 1350!$OMP END DO 1351 1352 DEALLOCATE (y) 1353!$OMP END PARALLEL 1354 END SUBROUTINE spline_interpolation 1355 1356! ************************************************************************************************** 1357!> \brief This routine is modeled after an algorithm from "Numerical Recipes in C" by Cambridge 1358!> University Press, pages 96-97. It was adapted for Fortran and for the problem at hand. 1359!> \param x ... 1360!> \param d2y_dx2 ... 1361!> \par History 1362!> OpenMP added: Aug 2016 MTucker 1363! ************************************************************************************************** 1364 SUBROUTINE initialize_spline_interpolation(x, d2y_dx2) 1365 1366 REAL(dp), INTENT(in) :: x(:) 1367 REAL(dp), INTENT(inout) :: d2y_dx2(:, :) 1368 1369 INTEGER :: index, Nx, P_i 1370 REAL(dp) :: temp1, temp2 1371 REAL(dp), ALLOCATABLE :: temp_array(:), y(:) 1372 1373 Nx = SIZE(x) 1374 1375!$OMP PARALLEL DEFAULT( NONE ) & 1376!$OMP SHARED( x, d2y_dx2, Nx ) & 1377!$OMP PRIVATE( temp_array, y & 1378!$OMP , index, temp1, temp2 & 1379!$OMP ) 1380 1381 ALLOCATE (temp_array(Nx), y(Nx)) 1382 1383!$OMP DO 1384 DO P_i = 1, Nx 1385 !! In the Soler method, the polynomials that are interpolated are Kronecker delta funcions 1386 !! at a particular q point. So, we set all y values to 0 except the one corresponding to 1387 !! the particular function P_i. 1388 !! ---------------------------------------------------------------------------------------- 1389 y = 0.0_dp 1390 y(P_i) = 1.0_dp 1391 !! ---------------------------------------------------------------------------------------- 1392 1393 d2y_dx2(P_i, 1) = 0.0_dp 1394 temp_array(1) = 0.0_dp 1395 DO index = 2, Nx - 1 1396 temp1 = (x(index) - x(index - 1))/(x(index + 1) - x(index - 1)) 1397 temp2 = temp1*d2y_dx2(P_i, index - 1) + 2.0_dp 1398 d2y_dx2(P_i, index) = (temp1 - 1.0_dp)/temp2 1399 temp_array(index) = (y(index + 1) - y(index))/(x(index + 1) - x(index)) & 1400 - (y(index) - y(index - 1))/(x(index) - x(index - 1)) 1401 temp_array(index) = (6.0_dp*temp_array(index)/(x(index + 1) - x(index - 1)) & 1402 - temp1*temp_array(index - 1))/temp2 1403 END DO 1404 d2y_dx2(P_i, Nx) = 0.0_dp 1405 DO index = Nx - 1, 1, -1 1406 d2y_dx2(P_i, index) = d2y_dx2(P_i, index)*d2y_dx2(P_i, index + 1) + temp_array(index) 1407 END DO 1408 END DO 1409!$OMP END DO 1410 1411 DEALLOCATE (temp_array, y) 1412!$OMP END PARALLEL 1413 1414 END SUBROUTINE initialize_spline_interpolation 1415 1416 ! ************************************************************************************************** 1417 !! This routine is modeled after an algorithm from "Numerical Recipes in C" by Cambridge 1418 !! University Press, page 97. Adapted for Fortran and the problem at hand. This function is used to 1419 !! find the Phi_alpha_beta needed for equations 8 and 11 of SOLER. 1420! ************************************************************************************************** 1421!> \brief ... 1422!> \param k ... 1423!> \param kernel_of_k ... 1424!> \param dispersion_env ... 1425!> \par History 1426!> Optimised when adding OpenMP to vdw_energy (which calls this routine): Aug 2016 MTucker 1427! ************************************************************************************************** 1428 SUBROUTINE interpolate_kernel(k, kernel_of_k, dispersion_env) 1429 REAL(dp), INTENT(in) :: k 1430 REAL(dp), INTENT(out) :: kernel_of_k(:, :) 1431 TYPE(qs_dispersion_type), POINTER :: dispersion_env 1432 1433 CHARACTER(len=*), PARAMETER :: routineN = 'interpolate_kernel', & 1434 routineP = moduleN//':'//routineN 1435 1436 INTEGER :: k_i, Nr_points, q1_i, q2_i 1437 REAL(dp) :: A, B, C, const, D, dk, r_max 1438 REAL(dp), DIMENSION(:, :, :), POINTER :: d2phi_dk2, kernel 1439 1440 Nr_points = dispersion_env%nr_points 1441 r_max = dispersion_env%r_max 1442 dk = dispersion_env%dk 1443 kernel => dispersion_env%kernel 1444 d2phi_dk2 => dispersion_env%d2phi_dk2 1445 1446 !! Check to make sure that the kernel table we have is capable of dealing with this 1447 !! value of k. If k is larger than Nr_points*2*pi/r_max then we can't perform the 1448 !! interpolation. In that case, a kernel file should be generated with a larger number 1449 !! of radial points. 1450 !! ------------------------------------------------------------------------------------- 1451 CPASSERT(k < Nr_points*dk) 1452 !! ------------------------------------------------------------------------------------- 1453 kernel_of_k = 0.0_dp 1454 !! This integer division figures out which bin k is in since the kernel 1455 !! is set on a uniform grid. 1456 k_i = INT(k/dk) 1457 1458 !! Test to see if we are trying to interpolate a k that is one of the actual 1459 !! function points we have. The value is just the value of the function in that 1460 !! case. 1461 !! ---------------------------------------------------------------------------------------- 1462 IF (MOD(k, dk) == 0) THEN 1463 DO q1_i = 1, dispersion_env%Nqs 1464 DO q2_i = 1, q1_i 1465 kernel_of_k(q1_i, q2_i) = kernel(k_i, q1_i, q2_i) 1466 kernel_of_k(q2_i, q1_i) = kernel(k_i, q2_i, q1_i) 1467 END DO 1468 END DO 1469 RETURN 1470 END IF 1471 !! ---------------------------------------------------------------------------------------- 1472 !! If we are not on a function point then we carry out the interpolation 1473 !! ---------------------------------------------------------------------------------------- 1474 const = dk*dk/6.0_dp 1475 A = (dk*(k_i + 1.0_dp) - k)/dk 1476 B = (k - dk*k_i)/dk 1477 C = (A**3 - A)*const 1478 D = (B**3 - B)*const 1479 DO q1_i = 1, dispersion_env%Nqs 1480 DO q2_i = 1, q1_i 1481 kernel_of_k(q1_i, q2_i) = A*kernel(k_i, q1_i, q2_i) + B*kernel(k_i + 1, q1_i, q2_i) & 1482 + (C*d2phi_dk2(k_i, q1_i, q2_i) + D*d2phi_dk2(k_i + 1, q1_i, q2_i)) 1483 kernel_of_k(q2_i, q1_i) = kernel_of_k(q1_i, q2_i) 1484 END DO 1485 END DO 1486 1487 END SUBROUTINE interpolate_kernel 1488 1489! ************************************************************************************************** 1490!> \brief ... 1491!> \param k ... 1492!> \param dkernel_of_dk ... 1493!> \param dispersion_env ... 1494!> \par History 1495!> Optimised when adding OpenMP to vdw_energy (which calls this routine): Aug 2016 MTucker 1496! ************************************************************************************************** 1497 SUBROUTINE interpolate_dkernel_dk(k, dkernel_of_dk, dispersion_env) 1498 REAL(dp), INTENT(in) :: k 1499 REAL(dp), INTENT(out) :: dkernel_of_dk(:, :) 1500 TYPE(qs_dispersion_type), POINTER :: dispersion_env 1501 1502 CHARACTER(len=*), PARAMETER :: routineN = 'interpolate_dkernel_dk', & 1503 routineP = moduleN//':'//routineN 1504 1505 INTEGER :: k_i, Nr_points, q1_i, q2_i 1506 REAL(dp) :: A, B, dAdk, dBdk, dCdk, dDdk, dk, dk_6 1507 REAL(dp), DIMENSION(:, :, :), POINTER :: d2phi_dk2, kernel 1508 1509 Nr_points = dispersion_env%nr_points 1510 dk = dispersion_env%dk 1511 kernel => dispersion_env%kernel 1512 d2phi_dk2 => dispersion_env%d2phi_dk2 1513 dk_6 = dk/6.0_dp 1514 1515 CPASSERT(k < Nr_points*dk) 1516 1517 dkernel_of_dk = 0.0_dp 1518 k_i = INT(k/dk) 1519 1520 A = (dk*(k_i + 1.0_dp) - k)/dk 1521 B = (k - dk*k_i)/dk 1522 dAdk = -1.0_dp/dk 1523 dBdk = 1.0_dp/dk 1524 dCdk = -(3*A**2 - 1.0_dp)*dk_6 1525 dDdk = (3*B**2 - 1.0_dp)*dk_6 1526 DO q1_i = 1, dispersion_env%Nqs 1527 DO q2_i = 1, q1_i 1528 dkernel_of_dk(q1_i, q2_i) = dAdk*kernel(k_i, q1_i, q2_i) + dBdk*kernel(k_i + 1, q1_i, q2_i) & 1529 + dCdk*d2phi_dk2(k_i, q1_i, q2_i) + dDdk*d2phi_dk2(k_i + 1, q1_i, q2_i) 1530 dkernel_of_dk(q2_i, q1_i) = dkernel_of_dk(q1_i, q2_i) 1531 END DO 1532 END DO 1533 1534 END SUBROUTINE interpolate_dkernel_dk 1535 1536! ************************************************************************************************** 1537 1538END MODULE qs_dispersion_nonloc 1539