Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Calculation of non local dispersion functionals
10 : !> Some routines adapted from:
11 : !> Copyright (C) 2001-2009 Quantum ESPRESSO group
12 : !> Copyright (C) 2009 Brian Kolb, Timo Thonhauser - Wake Forest University
13 : !> This file is distributed under the terms of the
14 : !> GNU General Public License. See the file `License'
15 : !> in the root directory of the present distribution,
16 : !> or http://www.gnu.org/copyleft/gpl.txt .
17 : !> \author JGH
18 : ! **************************************************************************************************
19 : MODULE qs_dispersion_nonloc
20 : USE bibliography, ONLY: Dion2004,&
21 : Romanperez2009,&
22 : Sabatini2013,&
23 : cite_reference
24 : USE cp_files, ONLY: close_file,&
25 : open_file
26 : USE input_constants, ONLY: vdw_nl_DRSLL,&
27 : vdw_nl_LMKLL,&
28 : vdw_nl_RVV10,&
29 : xc_vdw_fun_nonloc
30 : USE kinds, ONLY: default_path_length,&
31 : dp
32 : USE mathconstants, ONLY: pi,&
33 : rootpi
34 : USE message_passing, ONLY: mp_para_env_type
35 : USE pw_grid_types, ONLY: HALFSPACE,&
36 : pw_grid_type
37 : USE pw_methods, ONLY: pw_axpy,&
38 : pw_derive,&
39 : pw_transfer
40 : USE pw_pool_types, ONLY: pw_pool_type
41 : USE pw_types, ONLY: pw_c1d_gs_type,&
42 : pw_r3d_rs_type
43 : USE qs_dispersion_types, ONLY: qs_dispersion_type
44 : USE virial_types, ONLY: virial_type
45 : #include "./base/base_uses.f90"
46 :
47 : IMPLICIT NONE
48 :
49 : PRIVATE
50 :
51 : REAL(KIND=dp), PARAMETER :: epsr = 1.e-12_dp
52 :
53 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_dispersion_nonloc'
54 :
55 : PUBLIC :: qs_dispersion_nonloc_init, calculate_dispersion_nonloc
56 :
57 : ! **************************************************************************************************
58 :
59 : CONTAINS
60 :
61 : ! **************************************************************************************************
62 : !> \brief ...
63 : !> \param dispersion_env ...
64 : !> \param para_env ...
65 : ! **************************************************************************************************
66 50 : SUBROUTINE qs_dispersion_nonloc_init(dispersion_env, para_env)
67 : TYPE(qs_dispersion_type), POINTER :: dispersion_env
68 : TYPE(mp_para_env_type), POINTER :: para_env
69 :
70 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_dispersion_nonloc_init'
71 :
72 : CHARACTER(LEN=default_path_length) :: filename
73 : INTEGER :: funit, handle, ipair, itable, nqs, &
74 : nr_points, vdw_type
75 :
76 50 : CALL timeset(routineN, handle)
77 :
78 50 : SELECT CASE (dispersion_env%nl_type)
79 : CASE DEFAULT
80 0 : CPABORT("Unknown vdW-DF functional")
81 : CASE (vdw_nl_DRSLL, vdw_nl_LMKLL)
82 34 : CALL cite_reference(Dion2004)
83 : CASE (vdw_nl_RVV10)
84 50 : CALL cite_reference(Sabatini2013)
85 : END SELECT
86 50 : CALL cite_reference(RomanPerez2009)
87 :
88 50 : vdw_type = dispersion_env%type
89 50 : SELECT CASE (vdw_type)
90 : CASE DEFAULT
91 : ! do nothing
92 : CASE (xc_vdw_fun_nonloc)
93 : ! setup information on non local functionals
94 50 : filename = dispersion_env%kernel_file_name
95 50 : IF (para_env%is_source()) THEN
96 : ! Read the kernel information from file "filename"
97 25 : CALL open_file(file_name=filename, unit_number=funit, file_form="FORMATTED")
98 25 : READ (funit, *) nqs, nr_points
99 25 : READ (funit, *) dispersion_env%r_max
100 : END IF
101 50 : CALL para_env%bcast(nqs)
102 50 : CALL para_env%bcast(nr_points)
103 50 : CALL para_env%bcast(dispersion_env%r_max)
104 350 : ALLOCATE (dispersion_env%q_mesh(nqs), dispersion_env%kernel_table(nqs*(nqs + 1)/2, 0:nr_points, 2))
105 50 : dispersion_env%nqs = nqs
106 50 : dispersion_env%nr_points = nr_points
107 50 : IF (para_env%is_source()) THEN
108 : !! Read in the values of the q points used to generate this kernel
109 525 : READ (funit, "(1p, 4e23.14)") dispersion_env%q_mesh
110 : ! The file stores kernel values followed by second derivatives, both in
111 : ! lower-triangular pair order. Keep this immutable table in packed form.
112 75 : DO itable = 1, 2
113 10575 : DO ipair = 1, nqs*(nqs + 1)/2
114 10773050 : READ (funit, "(1p, 4e23.14)") dispersion_env%kernel_table(ipair, 0:nr_points, itable)
115 : END DO
116 : END DO
117 25 : CALL close_file(unit_number=funit)
118 : END IF
119 2050 : CALL para_env%bcast(dispersion_env%q_mesh)
120 43255250 : CALL para_env%bcast(dispersion_env%kernel_table)
121 : ! 2nd derivates for interpolation
122 200 : ALLOCATE (dispersion_env%d2y_dx2(nqs, nqs))
123 50 : CALL initialize_spline_interpolation(dispersion_env%q_mesh, dispersion_env%d2y_dx2)
124 : !
125 50 : dispersion_env%q_cut = dispersion_env%q_mesh(nqs)
126 50 : dispersion_env%q_min = dispersion_env%q_mesh(1)
127 100 : dispersion_env%dk = 2.0_dp*pi/dispersion_env%r_max
128 :
129 : END SELECT
130 :
131 50 : CALL timestop(handle)
132 :
133 50 : END SUBROUTINE qs_dispersion_nonloc_init
134 :
135 : ! **************************************************************************************************
136 : !> \brief Calculates the non-local vdW functional using the method of Soler
137 : !> For spin polarized cases we use E(a,b) = E(a+b), i.e. total density
138 : !> \param vxc_rho ...
139 : !> \param rho_r ...
140 : !> \param rho_g ...
141 : !> \param edispersion ...
142 : !> \param dispersion_env ...
143 : !> \param energy_only ...
144 : !> \param pw_pool ...
145 : !> \param xc_pw_pool ...
146 : !> \param para_env ...
147 : !> \param virial ...
148 : ! **************************************************************************************************
149 442 : SUBROUTINE calculate_dispersion_nonloc(vxc_rho, rho_r, rho_g, edispersion, &
150 : dispersion_env, energy_only, pw_pool, xc_pw_pool, para_env, virial)
151 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: vxc_rho, rho_r
152 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
153 : REAL(KIND=dp), INTENT(OUT) :: edispersion
154 : TYPE(qs_dispersion_type), POINTER :: dispersion_env
155 : LOGICAL, INTENT(IN) :: energy_only
156 : TYPE(pw_pool_type), POINTER :: pw_pool, xc_pw_pool
157 : TYPE(mp_para_env_type), POINTER :: para_env
158 : TYPE(virial_type), OPTIONAL, POINTER :: virial
159 :
160 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_dispersion_nonloc'
161 : INTEGER, DIMENSION(3, 3), PARAMETER :: nd = RESHAPE([1, 0, 0, 0, 1, 0, 0, 0, 1], [3, 3])
162 :
163 : INTEGER :: handle, handle_fft, i, i_grid, idir, &
164 : ispin, nl_type, np, nspin, p, q, r, s
165 442 : INTEGER, ALLOCATABLE, DIMENSION(:) :: q_low
166 : INTEGER, DIMENSION(1:3) :: hi, lo, n
167 : LOGICAL :: use_virial
168 : REAL(KIND=dp) :: b_value, beta, Ec_nl, sumnp
169 442 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: dq0_dgradrho, dq0_drho, hpot, q0, rho, &
170 442 : theta_scale
171 442 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: drho, spline_coeff, u_contract
172 442 : REAL(KIND=dp), CONTIGUOUS, POINTER :: tmp_1d(:), vxc_1d(:)
173 : TYPE(pw_c1d_gs_type) :: div_g, rho_tot_g, tmp_g, vxc_g
174 442 : TYPE(pw_c1d_gs_type), ALLOCATABLE, DIMENSION(:) :: thetas_g
175 : TYPE(pw_grid_type), POINTER :: grid
176 : TYPE(pw_r3d_rs_type) :: tmp_r, vxc_r
177 :
178 442 : CALL timeset(routineN, handle)
179 :
180 442 : CPASSERT(ASSOCIATED(rho_r))
181 442 : CPASSERT(ASSOCIATED(rho_g))
182 442 : CPASSERT(ASSOCIATED(pw_pool))
183 :
184 442 : IF (PRESENT(virial)) THEN
185 436 : use_virial = virial%pv_calculate .AND. (.NOT. virial%pv_numer)
186 : ELSE
187 : use_virial = .FALSE.
188 : END IF
189 : IF (use_virial) THEN
190 136 : CPASSERT(.NOT. energy_only)
191 : END IF
192 442 : IF (.NOT. energy_only) THEN
193 436 : CPASSERT(ASSOCIATED(vxc_rho))
194 : END IF
195 :
196 442 : nl_type = dispersion_env%nl_type
197 :
198 442 : b_value = dispersion_env%b_value
199 442 : beta = 0.03125_dp*(3.0_dp/(b_value**2.0_dp))**0.75_dp
200 442 : nspin = SIZE(rho_r)
201 :
202 : ! temporary arrays for FFT
203 442 : CALL pw_pool%create_pw(tmp_g)
204 442 : CALL pw_pool%create_pw(tmp_r)
205 :
206 : ! Sum the spin densities on the vdW grid before transforming or differentiating.
207 442 : CALL pw_pool%create_pw(rho_tot_g)
208 442 : CALL pw_transfer(rho_g(1), rho_tot_g)
209 462 : DO ispin = 2, nspin
210 20 : CALL pw_transfer(rho_g(ispin), tmp_g)
211 462 : CALL pw_axpy(tmp_g, rho_tot_g, 1._dp)
212 : END DO
213 442 : CALL pw_transfer(rho_tot_g, tmp_r)
214 :
215 1768 : np = SIZE(tmp_r%array)
216 442 : tmp_1d(1:np) => tmp_r%array
217 2210 : ALLOCATE (rho(np), drho(np, 3))
218 1768 : DO i = 1, 3
219 1326 : lo(i) = LBOUND(tmp_r%array, i)
220 1326 : hi(i) = UBOUND(tmp_r%array, i)
221 1768 : n(i) = hi(i) - lo(i) + 1
222 : END DO
223 : !$OMP PARALLEL DO DEFAULT(NONE) &
224 442 : !$OMP SHARED(n, lo, rho, tmp_r) PRIVATE(s) COLLAPSE(3)
225 : DO r = 0, n(3) - 1
226 : DO q = 0, n(2) - 1
227 : DO p = 0, n(1) - 1
228 : s = r*n(2)*n(1) + q*n(1) + p + 1
229 : rho(s) = tmp_r%array(p + lo(1), q + lo(2), r + lo(3))
230 : END DO
231 : END DO
232 : END DO
233 : !$OMP END PARALLEL DO
234 1768 : DO idir = 1, 3
235 1326 : CALL pw_transfer(rho_tot_g, tmp_g)
236 1326 : CALL pw_derive(tmp_g, nd(:, idir))
237 1326 : CALL pw_transfer(tmp_g, tmp_r)
238 : !$OMP PARALLEL DO DEFAULT(NONE) &
239 1768 : !$OMP SHARED(idir, n, lo, drho, tmp_r) PRIVATE(s) COLLAPSE(3)
240 : DO r = 0, n(3) - 1
241 : DO q = 0, n(2) - 1
242 : DO p = 0, n(1) - 1
243 : s = r*n(2)*n(1) + q*n(1) + p + 1
244 : drho(s, idir) = tmp_r%array(p + lo(1), q + lo(2), r + lo(3))
245 : END DO
246 : END DO
247 : END DO
248 : !$OMP END PARALLEL DO
249 : END DO
250 442 : CALL pw_pool%give_back_pw(rho_tot_g)
251 :
252 : !! ---------------------------------------------------------------------------------
253 : !! Find the value of q0 for all assigned grid points. q is defined in equations
254 : !! 11 and 12 of DION and q0 is the saturated version of q defined in equation
255 : !! 5 of SOLER. This routine also returns the derivatives of the q0s with respect
256 : !! to the charge-density and the gradient of the charge-density. These are needed
257 : !! for the potential calculated below.
258 : !! ---------------------------------------------------------------------------------
259 :
260 442 : IF (energy_only) THEN
261 12 : ALLOCATE (q0(np))
262 0 : SELECT CASE (nl_type)
263 : CASE DEFAULT
264 0 : CPABORT("Unknown vdW-DF functional")
265 : CASE (vdw_nl_DRSLL, vdw_nl_LMKLL)
266 0 : CALL get_q0_on_grid_eo_vdw(rho, drho, q0, dispersion_env)
267 : CASE (vdw_nl_RVV10)
268 6 : CALL get_q0_on_grid_eo_rvv10(rho, drho, q0, dispersion_env)
269 : END SELECT
270 : ELSE
271 1744 : ALLOCATE (q0(np), dq0_drho(np), dq0_dgradrho(np))
272 0 : SELECT CASE (nl_type)
273 : CASE DEFAULT
274 0 : CPABORT("Unknown vdW-DF functional")
275 : CASE (vdw_nl_DRSLL, vdw_nl_LMKLL)
276 254 : CALL get_q0_on_grid_vdw(rho, drho, q0, dq0_drho, dq0_dgradrho, dispersion_env)
277 : CASE (vdw_nl_RVV10)
278 436 : CALL get_q0_on_grid_rvv10(rho, drho, q0, dq0_drho, dq0_dgradrho, dispersion_env)
279 : END SELECT
280 : END IF
281 :
282 : ! Generate one theta channel directly in the FFT workspace. Only reciprocal-space
283 : ! channels must coexist for the convolution; no real-space np-by-nqs array is needed.
284 2652 : ALLOCATE (q_low(np), spline_coeff(np, 4), theta_scale(np))
285 442 : CALL prepare_splines(q0, rho, dispersion_env, q_low, spline_coeff, theta_scale)
286 10166 : ALLOCATE (thetas_g(dispersion_env%nqs))
287 442 : CALL timeset("vdW_theta_forward", handle_fft)
288 9282 : DO i = 1, dispersion_env%nqs
289 8840 : CALL build_theta(i, q_low, spline_coeff, theta_scale, dispersion_env, tmp_1d)
290 8840 : CALL pw_pool%create_pw(thetas_g(i))
291 9282 : CALL pw_transfer(tmp_r, thetas_g(i))
292 : END DO
293 442 : CALL timestop(handle_fft)
294 442 : DEALLOCATE (spline_coeff, theta_scale)
295 442 : grid => thetas_g(1)%pw_grid
296 : !! ---------------------------------------------------------------------------------------------
297 : !! Carry out the integration in equation 8 of SOLER. This also turns the thetas array into the
298 : !! precursor to the u_i(k) array which is inverse fourier transformed to get the u_i(r) functions
299 : !! of SOLER equation 11. Add the energy we find to the output variable etxc.
300 : !! --------------------------------------------------------------------------------------------------
301 442 : sumnp = np
302 442 : CALL para_env%sum(sumnp)
303 442 : IF (use_virial) THEN
304 : ! calculates kernel contribution to stress
305 136 : CALL vdW_energy(thetas_g, dispersion_env, Ec_nl, energy_only, virial)
306 54 : SELECT CASE (nl_type)
307 : CASE (vdw_nl_RVV10)
308 373384 : Ec_nl = 0.5_dp*Ec_nl + beta*SUM(rho(:))*grid%vol/sumnp
309 : END SELECT
310 : ! calculates energy contribution to stress
311 : ! potential contribution to stress is calculated together with other potentials (Hxc)
312 544 : DO idir = 1, 3
313 544 : virial%pv_xc(idir, idir) = virial%pv_xc(idir, idir) + Ec_nl
314 : END DO
315 : ELSE
316 306 : CALL vdW_energy(thetas_g, dispersion_env, Ec_nl, energy_only)
317 134 : SELECT CASE (nl_type)
318 : CASE (vdw_nl_RVV10)
319 662162 : Ec_nl = 0.5_dp*Ec_nl + beta*SUM(rho(:))*grid%vol/sumnp
320 : END SELECT
321 : END IF
322 442 : CALL para_env%sum(Ec_nl)
323 442 : IF (nl_type == vdw_nl_RVV10) Ec_nl = Ec_nl*dispersion_env%scale_rvv10
324 442 : edispersion = Ec_nl
325 :
326 442 : IF (energy_only) THEN
327 6 : DEALLOCATE (q0, q_low)
328 : ELSE
329 : ! Accumulate the two spline contractions and their endpoint values as each
330 : ! inverse FFT finishes. Keep the original potential-side knot convention.
331 1308 : ALLOCATE (u_contract(np, 4), hpot(np))
332 436 : !$OMP PARALLEL DO DEFAULT(NONE) SHARED(np, q_low, q0, dispersion_env, u_contract) PRIVATE(s)
333 : DO i_grid = 1, np
334 : s = q_low(i_grid)
335 : IF (s > 0 .AND. s < dispersion_env%nqs - 1) THEN
336 : IF (q0(i_grid) == dispersion_env%q_mesh(s + 1)) q_low(i_grid) = s + 1
337 : END IF
338 : u_contract(i_grid, :) = 0.0_dp
339 : END DO
340 : !$OMP END PARALLEL DO
341 436 : CALL timeset("vdW_theta_inverse", handle_fft)
342 9156 : DO i = 1, dispersion_env%nqs
343 8720 : CALL pw_transfer(thetas_g(i), tmp_r)
344 9156 : CALL accumulate_potential(i, q_low, dispersion_env, tmp_1d, u_contract)
345 : END DO
346 436 : CALL timestop(handle_fft)
347 :
348 : ! Write the local potential directly into its PW object.
349 436 : CALL pw_pool%create_pw(vxc_r)
350 436 : vxc_1d(1:np) => vxc_r%array
351 436 : IF (use_virial) THEN
352 136 : grid => tmp_g%pw_grid
353 : CALL get_potential(q0, dq0_drho, dq0_dgradrho, rho, q_low, u_contract, vxc_1d, hpot, &
354 136 : dispersion_env, drho, grid%dvol, virial)
355 : ELSE
356 : CALL get_potential(q0, dq0_drho, dq0_dgradrho, rho, q_low, u_contract, vxc_1d, hpot, &
357 300 : dispersion_env)
358 : END IF
359 436 : DEALLOCATE (u_contract, q_low, q0, dq0_drho, dq0_dgradrho)
360 182 : SELECT CASE (nl_type)
361 : CASE (vdw_nl_RVV10)
362 436 : !$OMP PARALLEL DO DEFAULT(NONE) SHARED(np, vxc_1d, hpot, beta, dispersion_env)
363 : DO i_grid = 1, np
364 : vxc_1d(i_grid) = (0.5_dp*vxc_1d(i_grid) + beta)*dispersion_env%scale_rvv10
365 : hpot(i_grid) = 0.5_dp*dispersion_env%scale_rvv10*hpot(i_grid)
366 : END DO
367 : !$OMP END PARALLEL DO
368 : END SELECT
369 436 : NULLIFY (vxc_1d)
370 : ! Sum the derivatives before the inverse FFT. Keep the real-space projection
371 : ! before transferring the potential to a possibly different XC grid.
372 436 : CALL pw_pool%create_pw(div_g)
373 1744 : DO idir = 1, 3
374 : !$OMP PARALLEL DO DEFAULT(NONE) &
375 1308 : !$OMP SHARED(n, lo, tmp_r, hpot, drho, idir) PRIVATE(s) COLLAPSE(3)
376 : DO r = 0, n(3) - 1
377 : DO q = 0, n(2) - 1
378 : DO p = 0, n(1) - 1
379 : s = r*n(2)*n(1) + q*n(1) + p + 1
380 : tmp_r%array(p + lo(1), q + lo(2), r + lo(3)) = hpot(s)*drho(s, idir)
381 : END DO
382 : END DO
383 : END DO
384 : !$OMP END PARALLEL DO
385 1308 : CALL pw_transfer(tmp_r, tmp_g)
386 1308 : CALL pw_derive(tmp_g, nd(:, idir))
387 1744 : IF (idir == 1) THEN
388 436 : CALL pw_transfer(tmp_g, div_g)
389 : ELSE
390 872 : CALL pw_axpy(tmp_g, div_g, 1._dp)
391 : END IF
392 : END DO
393 436 : CALL pw_transfer(div_g, tmp_r)
394 436 : CALL pw_pool%give_back_pw(div_g)
395 436 : CALL pw_axpy(tmp_r, vxc_r, -1._dp)
396 436 : CALL pw_transfer(vxc_r, tmp_g)
397 436 : CALL pw_pool%give_back_pw(vxc_r)
398 436 : CALL xc_pw_pool%create_pw(vxc_r)
399 436 : CALL xc_pw_pool%create_pw(vxc_g)
400 436 : CALL pw_transfer(tmp_g, vxc_g)
401 436 : CALL pw_transfer(vxc_g, vxc_r)
402 892 : DO ispin = 1, nspin
403 892 : CALL pw_axpy(vxc_r, vxc_rho(ispin), 1._dp)
404 : END DO
405 436 : CALL xc_pw_pool%give_back_pw(vxc_r)
406 872 : CALL xc_pw_pool%give_back_pw(vxc_g)
407 : END IF
408 :
409 : NULLIFY (tmp_1d)
410 :
411 9282 : DO i = 1, dispersion_env%nqs
412 9282 : CALL pw_pool%give_back_pw(thetas_g(i))
413 : END DO
414 442 : CALL pw_pool%give_back_pw(tmp_r)
415 442 : CALL pw_pool%give_back_pw(tmp_g)
416 :
417 442 : DEALLOCATE (rho, drho, thetas_g)
418 :
419 442 : CALL timestop(handle)
420 :
421 1326 : END SUBROUTINE calculate_dispersion_nonloc
422 :
423 : ! **************************************************************************************************
424 : !> \brief This routine carries out the integration of equation 8 of SOLER. It returns the non-local
425 : !> exchange-correlation energy and the u_alpha(k) arrays used to find the u_alpha(r) arrays via
426 : !> equations 11 and 12 in SOLER.
427 : !> energy contribution to stress is added in qs_force
428 : !> \param thetas_g ...
429 : !> \param dispersion_env ...
430 : !> \param vdW_xc_energy ...
431 : !> \param energy_only ...
432 : !> \param virial ...
433 : !> \par History
434 : !> OpenMP added: Aug 2016 MTucker
435 : ! **************************************************************************************************
436 442 : SUBROUTINE vdW_energy(thetas_g, dispersion_env, vdW_xc_energy, energy_only, virial)
437 : TYPE(pw_c1d_gs_type), DIMENSION(:), INTENT(IN) :: thetas_g
438 : TYPE(qs_dispersion_type), POINTER :: dispersion_env
439 : REAL(KIND=dp), INTENT(OUT) :: vdW_xc_energy
440 : LOGICAL, INTENT(IN) :: energy_only
441 : TYPE(virial_type), OPTIONAL, POINTER :: virial
442 :
443 : CHARACTER(LEN=*), PARAMETER :: routineN = 'vdW_energy'
444 : REAL(KIND=dp), PARAMETER :: eps_g2 = 1.0E-10_dp
445 :
446 : INTEGER :: handle, ig, iq, l, m, nl_type, nqs, &
447 : q1_i, q2_i
448 : LOGICAL :: use_virial
449 : REAL(KIND=dp) :: g, g2, g2_last, g_multiplier, gm
450 442 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: theta_im, theta_re, u_im, u_re
451 442 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dkernel_of_dk, kernel_of_k
452 : REAL(KIND=dp), DIMENSION(3, 3) :: virial_thread
453 : TYPE(pw_grid_type), POINTER :: grid
454 :
455 442 : CALL timeset(routineN, handle)
456 442 : nqs = dispersion_env%nqs
457 :
458 442 : use_virial = PRESENT(virial)
459 442 : virial_thread(:, :) = 0.0_dp ! always initialize to avoid floating point exceptions in OMP REDUCTION
460 :
461 442 : vdW_xc_energy = 0._dp
462 442 : grid => thetas_g(1)%pw_grid
463 :
464 442 : IF (grid%grid_span == HALFSPACE) THEN
465 : g_multiplier = 2._dp
466 : ELSE
467 442 : g_multiplier = 1._dp
468 : END IF
469 :
470 442 : nl_type = dispersion_env%nl_type
471 :
472 : !$OMP PARALLEL DEFAULT(NONE) &
473 : !$OMP SHARED(nqs, energy_only, grid, dispersion_env, use_virial, thetas_g, &
474 : !$OMP g_multiplier, nl_type) &
475 : !$OMP PRIVATE(g2_last, kernel_of_k, dkernel_of_dk, theta_re, theta_im, &
476 : !$OMP g2, g, iq, q2_i, u_re, u_im, q1_i, gm, l, m) &
477 442 : !$OMP REDUCTION(+:vdW_xc_energy, virial_thread)
478 :
479 : g2_last = HUGE(0._dp)
480 :
481 : ALLOCATE (kernel_of_k(nqs, nqs))
482 : IF (use_virial) ALLOCATE (dkernel_of_dk(nqs, nqs))
483 : ALLOCATE (theta_re(nqs), theta_im(nqs), u_re(nqs), u_im(nqs))
484 :
485 : !$OMP DO
486 : DO ig = 1, grid%ngpts_cut_local
487 : g2 = grid%gsq(ig)
488 : IF (ABS(g2 - g2_last) > eps_g2) THEN
489 : g2_last = g2
490 : g = SQRT(g2)
491 : IF (use_virial) THEN
492 : CALL interpolate_kernel(g, kernel_of_k, dispersion_env, dispersion_env%kernel_table, dkernel_of_dk)
493 : ELSE
494 : CALL interpolate_kernel(g, kernel_of_k, dispersion_env, dispersion_env%kernel_table)
495 : END IF
496 : END IF
497 : ! Save all inputs before in-place output. Vectorize over output channels;
498 : ! the input-channel accumulation order remains unchanged for every output.
499 : DO iq = 1, nqs
500 : theta_re(iq) = REAL(thetas_g(iq)%array(ig), KIND=dp)
501 : theta_im(iq) = AIMAG(thetas_g(iq)%array(ig))
502 : END DO
503 : u_re(:) = 0.0_dp
504 : u_im(:) = 0.0_dp
505 : DO q1_i = 1, nqs
506 : !$OMP SIMD
507 : DO q2_i = 1, nqs
508 : u_re(q2_i) = u_re(q2_i) + kernel_of_k(q2_i, q1_i)*theta_re(q1_i)
509 : u_im(q2_i) = u_im(q2_i) + kernel_of_k(q2_i, q1_i)*theta_im(q1_i)
510 : END DO
511 : !$OMP END SIMD
512 : END DO
513 : DO q2_i = 1, nqs
514 : IF (ig < grid%first_gne0) THEN
515 : vdW_xc_energy = vdW_xc_energy + (u_re(q2_i)*theta_re(q2_i) + u_im(q2_i)*theta_im(q2_i))
516 : ELSE
517 : vdW_xc_energy = vdW_xc_energy &
518 : + g_multiplier*(u_re(q2_i)*theta_re(q2_i) + u_im(q2_i)*theta_im(q2_i))
519 : END IF
520 : IF (.NOT. energy_only) thetas_g(q2_i)%array(ig) = CMPLX(u_re(q2_i), u_im(q2_i), KIND=dp)
521 : END DO
522 :
523 : IF (use_virial .AND. ig >= grid%first_gne0) THEN
524 : ! Reduce over all channel pairs before assembling the stress tensor.
525 : gm = 0.0_dp
526 : DO q2_i = 1, nqs
527 : !$OMP SIMD REDUCTION(+:gm)
528 : DO q1_i = 1, nqs
529 : gm = gm + dkernel_of_dk(q1_i, q2_i) &
530 : *(theta_re(q1_i)*theta_re(q2_i) + theta_im(q1_i)*theta_im(q2_i))
531 : END DO
532 : !$OMP END SIMD
533 : END DO
534 : gm = 0.5_dp*g_multiplier*grid%vol*gm
535 : IF (nl_type == vdw_nl_RVV10) gm = 0.5_dp*gm
536 : DO l = 1, 3
537 : DO m = 1, l
538 : virial_thread(l, m) = virial_thread(l, m) - gm*(grid%g(l, ig)*grid%g(m, ig))/g
539 : END DO
540 : END DO
541 : END IF
542 : END DO
543 : !$OMP END DO
544 :
545 : DEALLOCATE (theta_re, theta_im, u_re, u_im, kernel_of_k)
546 : IF (use_virial) DEALLOCATE (dkernel_of_dk)
547 :
548 : !$OMP END PARALLEL
549 :
550 442 : vdW_xc_energy = vdW_xc_energy*grid%vol*0.5_dp
551 :
552 442 : IF (use_virial) THEN
553 544 : DO l = 1, 3
554 816 : DO m = 1, (l - 1)
555 408 : virial%pv_xc(l, m) = virial%pv_xc(l, m) + virial_thread(l, m)
556 816 : virial%pv_xc(m, l) = virial%pv_xc(l, m)
557 : END DO
558 408 : m = l
559 544 : virial%pv_xc(l, m) = virial%pv_xc(l, m) + virial_thread(l, m)
560 : END DO
561 : END IF
562 :
563 442 : CALL timestop(handle)
564 :
565 442 : END SUBROUTINE vdW_energy
566 :
567 : ! **************************************************************************************************
568 : !> \brief This routine finds the non-local correlation contribution to the potential
569 : !> (i.e. the derivative of the non-local piece of the energy with respect to
570 : !> density) given in SOLER equation 10. The u_alpha(k) functions were found
571 : !> while calculating the energy. Their spline contractions are accumulated during the inverse FFTs.
572 : !> Most of the required derivatives were calculated in the "get_q0_on_grid"
573 : !> routine, but the derivative of the interpolation polynomials, P_alpha(q),
574 : !> (SOLER equation 3) with respect to q is interpolated here, along with the
575 : !> polynomials themselves.
576 : !> \param q0 ...
577 : !> \param dq0_drho ...
578 : !> \param dq0_dgradrho ...
579 : !> \param total_rho ...
580 : !> \param q_low_grid Lower spline endpoint at each grid point
581 : !> \param u_contract Second-derivative contractions and endpoint values
582 : !> \param potential ...
583 : !> \param h_prefactor ...
584 : !> \param dispersion_env ...
585 : !> \param drho ...
586 : !> \param dvol ...
587 : !> \param virial ...
588 : !> \par History
589 : !> OpenMP added: Aug 2016 MTucker
590 : ! **************************************************************************************************
591 436 : SUBROUTINE get_potential(q0, dq0_drho, dq0_dgradrho, total_rho, q_low_grid, u_contract, potential, h_prefactor, &
592 436 : dispersion_env, drho, dvol, virial)
593 :
594 : REAL(dp), DIMENSION(:), INTENT(in) :: q0, dq0_drho, dq0_dgradrho, total_rho
595 : INTEGER, DIMENSION(:), INTENT(IN) :: q_low_grid
596 : REAL(dp), DIMENSION(:, :), INTENT(in) :: u_contract
597 : REAL(dp), DIMENSION(:), INTENT(out) :: potential, h_prefactor
598 : TYPE(qs_dispersion_type), POINTER :: dispersion_env
599 : REAL(dp), DIMENSION(:, :), INTENT(in), OPTIONAL :: drho
600 : REAL(dp), INTENT(IN), OPTIONAL :: dvol
601 : TYPE(virial_type), OPTIONAL, POINTER :: virial
602 :
603 : CHARACTER(len=*), PARAMETER :: routineN = 'get_potential'
604 :
605 : INTEGER :: handle, i_grid, l, m, nl_type, nqs, &
606 : q_hi, q_low
607 : LOGICAL :: use_virial
608 : REAL(dp) :: a, b, b_value, c, const, d, dq, dq_6, e, &
609 : f, prefactor, tmp_1_2, tmp_1_4, &
610 : tmp_3_4, u_dp_dq0, u_p
611 : REAL(dp), DIMENSION(3, 3) :: virial_thread
612 436 : REAL(dp), DIMENSION(:), POINTER :: q_mesh
613 :
614 436 : CALL timeset(routineN, handle)
615 :
616 436 : use_virial = PRESENT(virial)
617 436 : CPASSERT(.NOT. use_virial .OR. PRESENT(drho))
618 436 : CPASSERT(.NOT. use_virial .OR. PRESENT(dvol))
619 :
620 436 : virial_thread(:, :) = 0.0_dp ! always initialize to avoid floating point exceptions in OMP REDUCTION
621 436 : b_value = dispersion_env%b_value
622 436 : const = 1.0_dp/(3.0_dp*b_value**(3.0_dp/2.0_dp)*pi**(5.0_dp/4.0_dp))
623 :
624 436 : q_mesh => dispersion_env%q_mesh
625 436 : nqs = dispersion_env%nqs
626 436 : nl_type = dispersion_env%nl_type
627 :
628 : !$OMP PARALLEL DEFAULT(NONE) &
629 : !$OMP SHARED(nqs, u_contract, q_low_grid, q_mesh, q0, nl_type, potential, h_prefactor, &
630 : !$OMP dq0_drho, dq0_dgradrho, total_rho, const, use_virial, drho, dvol, virial) &
631 : !$OMP PRIVATE(q_low, q_hi, dq, dq_6, A, b, c, d, e, f, u_p, u_dp_dq0, &
632 : !$OMP prefactor, l, m, tmp_1_2, tmp_1_4, tmp_3_4) &
633 436 : !$OMP REDUCTION(+:virial_thread)
634 :
635 : !$OMP DO
636 : DO i_grid = 1, SIZE(q0)
637 : potential(i_grid) = 0.0_dp
638 : h_prefactor(i_grid) = 0.0_dp
639 : IF (nl_type == vdw_nl_RVV10 .AND. total_rho(i_grid) <= epsr) CYCLE
640 : q_low = q_low_grid(i_grid)
641 : q_hi = q_low + 1
642 :
643 : dq = q_mesh(q_hi) - q_mesh(q_low)
644 : dq_6 = dq/6.0_dp
645 :
646 : a = (q_mesh(q_hi) - q0(i_grid))/dq
647 : b = (q0(i_grid) - q_mesh(q_low))/dq
648 : c = (a**3 - a)*dq*dq_6
649 : d = (b**3 - b)*dq*dq_6
650 : e = (3.0_dp*a**2 - 1.0_dp)*dq_6
651 : f = (3.0_dp*b**2 - 1.0_dp)*dq_6
652 :
653 : u_p = a*u_contract(i_grid, 3) + b*u_contract(i_grid, 4) &
654 : + c*u_contract(i_grid, 1) + d*u_contract(i_grid, 2)
655 : u_dp_dq0 = (u_contract(i_grid, 4) - u_contract(i_grid, 3))/dq &
656 : - e*u_contract(i_grid, 1) + f*u_contract(i_grid, 2)
657 :
658 : !! The first term in equation 13 of SOLER
659 : SELECT CASE (nl_type)
660 : CASE DEFAULT
661 : CPABORT("Unknown vdW-DF functional")
662 : CASE (vdw_nl_DRSLL, vdw_nl_LMKLL)
663 : potential(i_grid) = u_p + u_dp_dq0*dq0_drho(i_grid)
664 : prefactor = u_dp_dq0*dq0_dgradrho(i_grid)
665 : CASE (vdw_nl_RVV10)
666 : tmp_1_2 = SQRT(total_rho(i_grid))
667 : tmp_1_4 = SQRT(tmp_1_2)
668 : tmp_3_4 = tmp_1_4*tmp_1_4*tmp_1_4
669 : potential(i_grid) = const*0.75_dp/tmp_1_4*u_p + const*tmp_3_4*u_dp_dq0*dq0_drho(i_grid)
670 : prefactor = const*tmp_3_4*u_dp_dq0*dq0_dgradrho(i_grid)
671 : END SELECT
672 : IF (q0(i_grid) /= q_mesh(nqs)) THEN
673 : h_prefactor(i_grid) = prefactor
674 : END IF
675 :
676 : ! The saturation guard applies only to h_prefactor, not to the virial.
677 : IF (use_virial .AND. ABS(prefactor) > 0.0_dp) THEN
678 : IF (nl_type == vdw_nl_RVV10) prefactor = 0.5_dp*prefactor
679 : prefactor = prefactor*dvol
680 : DO l = 1, 3
681 : DO m = 1, l
682 : virial_thread(l, m) = virial_thread(l, m) - prefactor*drho(i_grid, l)*drho(i_grid, m)
683 : END DO
684 : END DO
685 : END IF
686 : END DO ! i_grid = 1, SIZE(q0)
687 : !$OMP END DO
688 :
689 : !$OMP END PARALLEL
690 :
691 436 : IF (use_virial) THEN
692 544 : DO l = 1, 3
693 816 : DO m = 1, (l - 1)
694 408 : virial%pv_xc(l, m) = virial%pv_xc(l, m) + virial_thread(l, m)
695 816 : virial%pv_xc(m, l) = virial%pv_xc(l, m)
696 : END DO
697 408 : m = l
698 544 : virial%pv_xc(l, m) = virial%pv_xc(l, m) + virial_thread(l, m)
699 : END DO
700 : END IF
701 :
702 436 : CALL timestop(handle)
703 436 : END SUBROUTINE get_potential
704 :
705 : ! **************************************************************************************************
706 : !> \brief calculates exponent = sum(from i=1 to hi, ((alpha)**i)/i) ) without <<< calling power >>>
707 : !> \param hi = upper index for sum
708 : !> \param alpha ...
709 : !> \param exponent = output value
710 : !> \par History
711 : !> Created: MTucker, Aug 2016
712 : ! **************************************************************************************************
713 72326 : ELEMENTAL SUBROUTINE calculate_exponent(hi, alpha, exponent)
714 : INTEGER, INTENT(in) :: hi
715 : REAL(dp), INTENT(in) :: alpha
716 : REAL(dp), INTENT(out) :: exponent
717 :
718 : INTEGER :: i
719 : REAL(dp) :: multiplier
720 :
721 72326 : multiplier = alpha
722 72326 : exponent = alpha
723 :
724 867912 : DO i = 2, hi
725 795586 : multiplier = multiplier*alpha
726 867912 : exponent = exponent + (multiplier/i)
727 : END DO
728 72326 : END SUBROUTINE calculate_exponent
729 :
730 : ! **************************************************************************************************
731 : !> \brief calculate exponent = sum(from i=1 to hi, ((alpha)**i)/i) ) without calling power
732 : !> also calculates derivative using similar series
733 : !> \param hi = upper index for sum
734 : !> \param alpha ...
735 : !> \param exponent = output value
736 : !> \param derivative ...
737 : !> \par History
738 : !> Created: MTucker, Aug 2016
739 : ! **************************************************************************************************
740 5214022 : ELEMENTAL SUBROUTINE calculate_exponent_derivative(hi, alpha, exponent, derivative)
741 : INTEGER, INTENT(in) :: hi
742 : REAL(dp), INTENT(in) :: alpha
743 : REAL(dp), INTENT(out) :: exponent, derivative
744 :
745 : INTEGER :: i
746 : REAL(dp) :: multiplier
747 :
748 5214022 : derivative = 0.0d0
749 5214022 : multiplier = 1.0d0
750 5214022 : exponent = 0.0d0
751 :
752 67782286 : DO i = 1, hi
753 62568264 : derivative = derivative + multiplier
754 62568264 : multiplier = multiplier*alpha
755 67782286 : exponent = exponent + (multiplier/i)
756 : END DO
757 5214022 : END SUBROUTINE calculate_exponent_derivative
758 :
759 : !! This routine first calculates the q value defined in (DION equations 11 and 12), then
760 : !! saturates it according to (SOLER equation 5).
761 : ! **************************************************************************************************
762 : !> \brief This routine first calculates the q value defined in (DION equations 11 and 12), then
763 : !> saturates it according to (SOLER equation 5).
764 : !> \param total_rho ...
765 : !> \param gradient_rho ...
766 : !> \param q0 ...
767 : !> \param dq0_drho ...
768 : !> \param dq0_dgradrho ...
769 : !> \param dispersion_env ...
770 : ! **************************************************************************************************
771 254 : SUBROUTINE get_q0_on_grid_vdw(total_rho, gradient_rho, q0, dq0_drho, dq0_dgradrho, dispersion_env)
772 : !!
773 : !! more specifically it calculates the following
774 : !!
775 : !! q0(ir) = q0 as defined above
776 : !! dq0_drho(ir) = total_rho * d q0 /d rho
777 : !! dq0_dgradrho = total_rho / |gradient_rho| * d q0 / d |gradient_rho|
778 : !!
779 : REAL(dp), INTENT(IN) :: total_rho(:), gradient_rho(:, :)
780 : REAL(dp), INTENT(OUT) :: q0(:), dq0_drho(:), dq0_dgradrho(:)
781 : TYPE(qs_dispersion_type), POINTER :: dispersion_env
782 :
783 : INTEGER, PARAMETER :: m_cut = 12
784 : REAL(dp), PARAMETER :: LDA_A = 0.031091_dp, LDA_a1 = 0.2137_dp, LDA_b1 = 7.5957_dp, &
785 : LDA_b2 = 3.5876_dp, LDA_b3 = 1.6382_dp, LDA_b4 = 0.49294_dp
786 :
787 : INTEGER :: i_grid
788 : REAL(dp) :: dq0_dq, exponent, gradient_correction, &
789 : kF, LDA_1, LDA_2, q, q__q_cut, q_cut, &
790 : q_min, r_s, sqrt_r_s, Z_ab
791 :
792 254 : q_cut = dispersion_env%q_cut
793 254 : q_min = dispersion_env%q_min
794 254 : SELECT CASE (dispersion_env%nl_type)
795 : CASE DEFAULT
796 0 : CPABORT("Unknown vdW-DF functional")
797 : CASE (vdw_nl_DRSLL)
798 48 : Z_ab = -0.8491_dp
799 : CASE (vdw_nl_LMKLL)
800 254 : Z_ab = -1.887_dp
801 : END SELECT
802 :
803 : !$OMP PARALLEL DO DEFAULT(NONE) &
804 : !$OMP SHARED(total_rho, gradient_rho, q0, dq0_drho, dq0_dgradrho, q_cut, q_min, Z_ab) &
805 : !$OMP PRIVATE(dq0_dq, exponent, gradient_correction, kF, LDA_1, LDA_2, q, q__q_cut, r_s, sqrt_r_s) &
806 254 : !$OMP SCHEDULE(STATIC)
807 : DO i_grid = 1, SIZE(total_rho)
808 : q0(i_grid) = q_cut
809 : dq0_drho(i_grid) = 0.0_dp
810 : dq0_dgradrho(i_grid) = 0.0_dp
811 :
812 : !! This prevents numerical problems. If the charge density is negative (an
813 : !! unphysical situation), we simply treat it as very small. In that case,
814 : !! q0 will be very large and will be saturated. For a saturated q0 the derivative
815 : !! dq0_dq will be 0 so we set q0 = q_cut and dq0_drho = dq0_dgradrho = 0 and go on
816 : !! to the next point.
817 : !! ------------------------------------------------------------------------------------
818 : IF (total_rho(i_grid) < epsr) CYCLE
819 : !! ------------------------------------------------------------------------------------
820 : !! Calculate some intermediate values needed to find q
821 : !! ------------------------------------------------------------------------------------
822 : kF = (3.0_dp*pi*pi*total_rho(i_grid))**(1.0_dp/3.0_dp)
823 : r_s = (3.0_dp/(4.0_dp*pi*total_rho(i_grid)))**(1.0_dp/3.0_dp)
824 : sqrt_r_s = SQRT(r_s)
825 :
826 : gradient_correction = -Z_ab/(36.0_dp*kF*total_rho(i_grid)**2) &
827 : *(gradient_rho(i_grid, 1)**2 + gradient_rho(i_grid, 2)**2 + gradient_rho(i_grid, 3)**2)
828 :
829 : LDA_1 = 8.0_dp*pi/3.0_dp*(LDA_A*(1.0_dp + LDA_a1*r_s))
830 : 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)
831 : !! ---------------------------------------------------------------
832 : !! This is the q value defined in equations 11 and 12 of DION
833 : !! ---------------------------------------------------------------
834 : q = kF + LDA_1*LOG(1.0_dp + 1.0_dp/LDA_2) + gradient_correction
835 : !! ---------------------------------------------------------------
836 : !! Here, we calculate q0 by saturating q according to equation 5 of SOLER. Also, we find
837 : !! the derivative dq0_dq needed for the derivatives dq0_drho and dq0_dgradrh0 discussed below.
838 : !! ---------------------------------------------------------------------------------------
839 : q__q_cut = q/q_cut
840 : CALL calculate_exponent_derivative(m_cut, q__q_cut, exponent, dq0_dq)
841 : q0(i_grid) = q_cut*(1.0_dp - EXP(-exponent))
842 : dq0_dq = dq0_dq*EXP(-exponent)
843 : !! ---------------------------------------------------------------------------------------
844 : !! This is to handle a case with q0 too small. We simply set it to the smallest q value in
845 : !! out q_mesh. Hopefully this doesn't get used often (ever)
846 : !! ---------------------------------------------------------------------------------------
847 : IF (q0(i_grid) < q_min) THEN
848 : q0(i_grid) = q_min
849 : END IF
850 : !! ---------------------------------------------------------------------------------------
851 : !! Here we find derivatives. These are actually the density times the derivative of q0 with respect
852 : !! to rho and gradient_rho. The density factor comes in since we are really differentiating
853 : !! theta = (rho)*P(q0) with respect to density (or its gradient) which will be
854 : !! dtheta_drho = P(q0) + dP_dq0 * [rho * dq0_dq * dq_drho] and
855 : !! dtheta_dgradient_rho = dP_dq0 * [rho * dq0_dq * dq_dgradient_rho]
856 : !! The parts in square brackets are what is calculated here. The dP_dq0 term will be interpolated
857 : !! later. There should actually be a factor of the magnitude of the gradient in the gradient_rho derivative
858 : !! but that cancels out when we differentiate the magnitude of the gradient with respect to a particular
859 : !! component.
860 : !! ------------------------------------------------------------------------------------------------
861 :
862 : dq0_drho(i_grid) = dq0_dq*(kF/3.0_dp - 7.0_dp/3.0_dp*gradient_correction &
863 : - 8.0_dp*pi/9.0_dp*LDA_A*LDA_a1*r_s*LOG(1.0_dp + 1.0_dp/LDA_2) &
864 : + LDA_1/(LDA_2*(1.0_dp + LDA_2)) &
865 : *(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 &
866 : + 2.0_dp*LDA_b4/3.0_dp*r_s**2)))
867 :
868 : dq0_dgradrho(i_grid) = total_rho(i_grid)*dq0_dq*2.0_dp*(-Z_ab)/(36.0_dp*kF*total_rho(i_grid)**2)
869 :
870 : END DO
871 : !$OMP END PARALLEL DO
872 :
873 254 : END SUBROUTINE get_q0_on_grid_vdw
874 :
875 : ! **************************************************************************************************
876 : !> \brief ...
877 : !> \param total_rho ...
878 : !> \param gradient_rho ...
879 : !> \param q0 ...
880 : !> \param dq0_drho ...
881 : !> \param dq0_dgradrho ...
882 : !> \param dispersion_env ...
883 : ! **************************************************************************************************
884 182 : SUBROUTINE get_q0_on_grid_rvv10(total_rho, gradient_rho, q0, dq0_drho, dq0_dgradrho, dispersion_env)
885 : !!
886 : !! more specifically it calculates the following
887 : !!
888 : !! q0(ir) = q0 as defined above
889 : !! dq0_drho(ir) = total_rho * d q0 /d rho
890 : !! dq0_dgradrho = total_rho / |gradient_rho| * d q0 / d |gradient_rho|
891 : !!
892 : REAL(dp), INTENT(IN) :: total_rho(:), gradient_rho(:, :)
893 : REAL(dp), INTENT(OUT) :: q0(:), dq0_drho(:), dq0_dgradrho(:)
894 : TYPE(qs_dispersion_type), POINTER :: dispersion_env
895 :
896 : INTEGER, PARAMETER :: m_cut = 12
897 :
898 : INTEGER :: i_grid
899 : REAL(dp) :: b_value, C_value, dk_dn, dq0_dq, dw0_dn, &
900 : exponent, gmod2, k, mod_grad, q, &
901 : q__q_cut, q_cut, q_min, w0, wg2, wp2
902 :
903 182 : q_cut = dispersion_env%q_cut
904 182 : q_min = dispersion_env%q_min
905 182 : b_value = dispersion_env%b_value
906 182 : C_value = dispersion_env%c_value
907 :
908 : !$OMP PARALLEL DO DEFAULT(NONE) &
909 : !$OMP SHARED(total_rho, gradient_rho, q0, dq0_drho, dq0_dgradrho, q_cut, q_min, b_value, C_value) &
910 : !$OMP PRIVATE(dk_dn, dq0_dq, dw0_dn, exponent, gmod2, k, mod_grad, q, q__q_cut, w0, wg2, wp2) &
911 182 : !$OMP SCHEDULE(STATIC)
912 : DO i_grid = 1, SIZE(total_rho)
913 : q0(i_grid) = q_cut
914 : dq0_drho(i_grid) = 0.0_dp
915 : dq0_dgradrho(i_grid) = 0.0_dp
916 :
917 : gmod2 = gradient_rho(i_grid, 1)**2 + gradient_rho(i_grid, 2)**2 + gradient_rho(i_grid, 3)**2
918 :
919 : !if (total_rho(i_grid) > epsr .and. gmod2 > epsr) cycle
920 : IF (total_rho(i_grid) > epsr) THEN
921 :
922 : !! Calculate some intermediate values needed to find q
923 : !! ------------------------------------------------------------------------------------
924 : mod_grad = SQRT(gmod2)
925 :
926 : wp2 = 16.0_dp*pi*total_rho(i_grid)
927 : wg2 = 4_dp*C_value*(mod_grad/total_rho(i_grid))**4
928 :
929 : k = b_value*3.0_dp*pi*((total_rho(i_grid)/(9.0_dp*pi))**(1.0_dp/6.0_dp))
930 : w0 = SQRT(wg2 + wp2/3.0_dp)
931 :
932 : q = w0/k
933 :
934 : !! Here, we calculate q0 by saturating q according
935 : !! ---------------------------------------------------------------------------------------
936 : q__q_cut = q/q_cut
937 : CALL calculate_exponent_derivative(m_cut, q__q_cut, exponent, dq0_dq)
938 : q0(i_grid) = q_cut*(1.0_dp - EXP(-exponent))
939 : dq0_dq = dq0_dq*EXP(-exponent)
940 :
941 : !! ---------------------------------------------------------------------------------------
942 : IF (q0(i_grid) < q_min) THEN
943 : q0(i_grid) = q_min
944 : END IF
945 :
946 : !!---------------------------------Final values---------------------------------
947 : dw0_dn = 1.0_dp/(2.0_dp*w0)*(16.0_dp/3.0_dp*pi - 4.0_dp*wg2/total_rho(i_grid))
948 : dk_dn = k/(6.0_dp*total_rho(i_grid))
949 :
950 : dq0_drho(i_grid) = dq0_dq*1.0_dp/(k**2)*(dw0_dn*k - dk_dn*w0)
951 : ! wg2 is proportional to |gradient_rho|**4, so this limit is zero.
952 : IF (gmod2 > 0.0_dp) THEN
953 : dq0_dgradrho(i_grid) = dq0_dq*1.0_dp/(2.0_dp*k*w0)*4.0_dp*wg2/gmod2
954 : END IF
955 : END IF
956 :
957 : END DO
958 : !$OMP END PARALLEL DO
959 :
960 182 : END SUBROUTINE get_q0_on_grid_rvv10
961 :
962 : ! **************************************************************************************************
963 : !> \brief ...
964 : !> \param total_rho ...
965 : !> \param gradient_rho ...
966 : !> \param q0 ...
967 : !> \param dispersion_env ...
968 : ! **************************************************************************************************
969 0 : SUBROUTINE get_q0_on_grid_eo_vdw(total_rho, gradient_rho, q0, dispersion_env)
970 :
971 : REAL(dp), INTENT(IN) :: total_rho(:), gradient_rho(:, :)
972 : REAL(dp), INTENT(OUT) :: q0(:)
973 : TYPE(qs_dispersion_type), POINTER :: dispersion_env
974 :
975 : INTEGER, PARAMETER :: m_cut = 12
976 : REAL(dp), PARAMETER :: LDA_A = 0.031091_dp, LDA_a1 = 0.2137_dp, LDA_b1 = 7.5957_dp, &
977 : LDA_b2 = 3.5876_dp, LDA_b3 = 1.6382_dp, LDA_b4 = 0.49294_dp
978 :
979 : INTEGER :: i_grid
980 : REAL(dp) :: exponent, gradient_correction, kF, &
981 : LDA_1, LDA_2, q, q__q_cut, q_cut, &
982 : q_min, r_s, sqrt_r_s, Z_ab
983 :
984 0 : q_cut = dispersion_env%q_cut
985 0 : q_min = dispersion_env%q_min
986 0 : SELECT CASE (dispersion_env%nl_type)
987 : CASE DEFAULT
988 0 : CPABORT("Unknown vdW-DF functional")
989 : CASE (vdw_nl_DRSLL)
990 0 : Z_ab = -0.8491_dp
991 : CASE (vdw_nl_LMKLL)
992 0 : Z_ab = -1.887_dp
993 : END SELECT
994 :
995 : !$OMP PARALLEL DO DEFAULT(NONE) &
996 : !$OMP SHARED(total_rho, gradient_rho, q0, q_cut, q_min, Z_ab) &
997 : !$OMP PRIVATE(exponent, gradient_correction, kF, LDA_1, LDA_2, q, q__q_cut, r_s, sqrt_r_s) &
998 0 : !$OMP SCHEDULE(STATIC)
999 : DO i_grid = 1, SIZE(total_rho)
1000 : q0(i_grid) = q_cut
1001 : !! This prevents numerical problems. If the charge density is negative (an
1002 : !! unphysical situation), we simply treat it as very small. In that case,
1003 : !! q0 will be very large and will be saturated. For a saturated q0 the derivative
1004 : !! dq0_dq will be 0 so we set q0 = q_cut and dq0_drho = dq0_dgradrho = 0 and go on
1005 : !! to the next point.
1006 : !! ------------------------------------------------------------------------------------
1007 : IF (total_rho(i_grid) < epsr) CYCLE
1008 : !! ------------------------------------------------------------------------------------
1009 : !! Calculate some intermediate values needed to find q
1010 : !! ------------------------------------------------------------------------------------
1011 : kF = (3.0_dp*pi*pi*total_rho(i_grid))**(1.0_dp/3.0_dp)
1012 : r_s = (3.0_dp/(4.0_dp*pi*total_rho(i_grid)))**(1.0_dp/3.0_dp)
1013 : sqrt_r_s = SQRT(r_s)
1014 :
1015 : gradient_correction = -Z_ab/(36.0_dp*kF*total_rho(i_grid)**2) &
1016 : *(gradient_rho(i_grid, 1)**2 + gradient_rho(i_grid, 2)**2 + gradient_rho(i_grid, 3)**2)
1017 :
1018 : LDA_1 = 8.0_dp*pi/3.0_dp*(LDA_A*(1.0_dp + LDA_a1*r_s))
1019 : 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)
1020 : !! ------------------------------------------------------------------------------------
1021 : !! This is the q value defined in equations 11 and 12 of DION
1022 : !! ---------------------------------------------------------------
1023 : q = kF + LDA_1*LOG(1.0_dp + 1.0_dp/LDA_2) + gradient_correction
1024 :
1025 : !! ---------------------------------------------------------------
1026 : !! Here, we calculate q0 by saturating q according to equation 5 of SOLER. Also, we find
1027 : !! the derivative dq0_dq needed for the derivatives dq0_drho and dq0_dgradrh0 discussed below.
1028 : !! ---------------------------------------------------------------------------------------
1029 : q__q_cut = q/q_cut
1030 : CALL calculate_exponent(m_cut, q__q_cut, exponent)
1031 : q0(i_grid) = q_cut*(1.0_dp - EXP(-exponent))
1032 :
1033 : !! ---------------------------------------------------------------------------------------
1034 : !! This is to handle a case with q0 too small. We simply set it to the smallest q value in
1035 : !! out q_mesh. Hopefully this doesn't get used often (ever)
1036 : !! ---------------------------------------------------------------------------------------
1037 : IF (q0(i_grid) < q_min) THEN
1038 : q0(i_grid) = q_min
1039 : END IF
1040 : END DO
1041 : !$OMP END PARALLEL DO
1042 :
1043 0 : END SUBROUTINE get_q0_on_grid_eo_vdw
1044 :
1045 : ! **************************************************************************************************
1046 : !> \brief ...
1047 : !> \param total_rho ...
1048 : !> \param gradient_rho ...
1049 : !> \param q0 ...
1050 : !> \param dispersion_env ...
1051 : ! **************************************************************************************************
1052 6 : SUBROUTINE get_q0_on_grid_eo_rvv10(total_rho, gradient_rho, q0, dispersion_env)
1053 :
1054 : REAL(dp), INTENT(IN) :: total_rho(:), gradient_rho(:, :)
1055 : REAL(dp), INTENT(OUT) :: q0(:)
1056 : TYPE(qs_dispersion_type), POINTER :: dispersion_env
1057 :
1058 : INTEGER, PARAMETER :: m_cut = 12
1059 :
1060 : INTEGER :: i_grid
1061 : REAL(dp) :: b_value, C_value, exponent, gmod2, k, q, &
1062 : q__q_cut, q_cut, q_min, w0, wg2, wp2
1063 :
1064 6 : q_cut = dispersion_env%q_cut
1065 6 : q_min = dispersion_env%q_min
1066 6 : b_value = dispersion_env%b_value
1067 6 : C_value = dispersion_env%c_value
1068 :
1069 : !$OMP PARALLEL DO DEFAULT(NONE) &
1070 : !$OMP SHARED(total_rho, gradient_rho, q0, q_cut, q_min, b_value, C_value) &
1071 : !$OMP PRIVATE(exponent, gmod2, k, q, q__q_cut, w0, wg2, wp2) &
1072 6 : !$OMP SCHEDULE(STATIC)
1073 : DO i_grid = 1, SIZE(total_rho)
1074 : q0(i_grid) = q_cut
1075 :
1076 : gmod2 = gradient_rho(i_grid, 1)**2 + gradient_rho(i_grid, 2)**2 + gradient_rho(i_grid, 3)**2
1077 :
1078 : !if (total_rho(i_grid) > epsr .and. gmod2 > epsr) cycle
1079 : IF (total_rho(i_grid) > epsr) THEN
1080 :
1081 : !! Calculate some intermediate values needed to find q
1082 : !! ------------------------------------------------------------------------------------
1083 : wp2 = 16.0_dp*pi*total_rho(i_grid)
1084 : wg2 = 4_dp*C_value*(gmod2*gmod2)/(total_rho(i_grid)**4)
1085 :
1086 : k = b_value*3.0_dp*pi*((total_rho(i_grid)/(9.0_dp*pi))**(1.0_dp/6.0_dp))
1087 : w0 = SQRT(wg2 + wp2/3.0_dp)
1088 :
1089 : q = w0/k
1090 :
1091 : !! Here, we calculate q0 by saturating q according
1092 : !! ---------------------------------------------------------------------------------------
1093 : q__q_cut = q/q_cut
1094 : CALL calculate_exponent(m_cut, q__q_cut, exponent)
1095 : q0(i_grid) = q_cut*(1.0_dp - EXP(-exponent))
1096 :
1097 : IF (q0(i_grid) < q_min) THEN
1098 : q0(i_grid) = q_min
1099 : END IF
1100 :
1101 : END IF
1102 :
1103 : END DO
1104 : !$OMP END PARALLEL DO
1105 :
1106 6 : END SUBROUTINE get_q0_on_grid_eo_rvv10
1107 :
1108 : ! **************************************************************************************************
1109 : !> \brief Prepare the spline intervals, weights and density factors for streamed theta construction.
1110 : !> \param q0 Saturated interpolation points
1111 : !> \param rho Total density
1112 : !> \param dispersion_env Non-local functional parameters
1113 : !> \param q_low Lower spline endpoint; zero denotes an inactive rVV10 point
1114 : !> \param coeff Cubic spline coefficients a, b, c, d
1115 : !> \param theta_scale Density factor multiplying the interpolated spline
1116 : ! **************************************************************************************************
1117 442 : SUBROUTINE prepare_splines(q0, rho, dispersion_env, q_low, coeff, theta_scale)
1118 : REAL(dp), INTENT(IN) :: q0(:), rho(:)
1119 : TYPE(qs_dispersion_type), POINTER :: dispersion_env
1120 : INTEGER, INTENT(OUT) :: q_low(:)
1121 : REAL(dp), INTENT(OUT) :: coeff(:, :), theta_scale(:)
1122 :
1123 : INTEGER :: i, j, lower, nqs, upper
1124 : LOGICAL :: rvv10
1125 : REAL(dp) :: a, b, const, dx, dx2_6
1126 : REAL(dp), POINTER :: q_mesh(:)
1127 :
1128 442 : q_mesh => dispersion_env%q_mesh
1129 442 : nqs = dispersion_env%nqs
1130 442 : CPASSERT(nqs >= 2)
1131 442 : rvv10 = dispersion_env%nl_type == vdw_nl_RVV10
1132 442 : const = 1.0_dp/(3.0_dp*rootpi*dispersion_env%b_value**1.5_dp)/(pi**0.75_dp)
1133 : !$OMP PARALLEL DO DEFAULT(NONE) &
1134 : !$OMP SHARED(q0, rho, q_mesh, nqs, rvv10, const, q_low, coeff, theta_scale) &
1135 442 : !$OMP PRIVATE(j, lower, upper, A, b, dx, dx2_6) SCHEDULE(STATIC)
1136 : DO i = 1, SIZE(q0)
1137 : IF (rvv10 .AND. rho(i) <= epsr) THEN
1138 : q_low(i) = 0
1139 : coeff(i, :) = 0.0_dp
1140 : theta_scale(i) = 0.0_dp
1141 : CYCLE
1142 : END IF
1143 : lower = 1
1144 : upper = nqs
1145 : DO WHILE (upper - lower > 1)
1146 : j = (upper + lower)/2
1147 : IF (q0(i) > q_mesh(j)) THEN
1148 : lower = j
1149 : ELSE
1150 : upper = j
1151 : END IF
1152 : END DO
1153 : q_low(i) = lower
1154 : dx = q_mesh(upper) - q_mesh(lower)
1155 : dx2_6 = dx*dx/6.0_dp
1156 : a = (q_mesh(upper) - q0(i))/dx
1157 : b = (q0(i) - q_mesh(lower))/dx
1158 : coeff(i, 1) = a
1159 : coeff(i, 2) = b
1160 : coeff(i, 3) = (a**3 - a)*dx2_6
1161 : coeff(i, 4) = (b**3 - b)*dx2_6
1162 : IF (rvv10) THEN
1163 : theta_scale(i) = const*rho(i)**0.75_dp
1164 : ELSE
1165 : theta_scale(i) = rho(i)
1166 : END IF
1167 : END DO
1168 : !$OMP END PARALLEL DO
1169 442 : END SUBROUTINE prepare_splines
1170 :
1171 : ! **************************************************************************************************
1172 : !> \brief Generate one real-space theta channel directly in the existing FFT workspace.
1173 : !> \param iq Channel index
1174 : !> \param q_low Lower spline endpoint
1175 : !> \param coeff Spline coefficients
1176 : !> \param theta_scale Density factor
1177 : !> \param dispersion_env Non-local functional parameters
1178 : !> \param theta FFT input, viewed as a contiguous one-dimensional array
1179 : ! **************************************************************************************************
1180 8840 : SUBROUTINE build_theta(iq, q_low, coeff, theta_scale, dispersion_env, theta)
1181 : INTEGER, INTENT(IN) :: iq, q_low(:)
1182 : REAL(dp), INTENT(IN) :: coeff(:, :), theta_scale(:)
1183 : TYPE(qs_dispersion_type), POINTER :: dispersion_env
1184 : REAL(dp), INTENT(OUT) :: theta(:)
1185 :
1186 : INTEGER :: i, lower
1187 : REAL(dp) :: p
1188 : REAL(dp), POINTER :: d2y(:, :)
1189 :
1190 8840 : d2y => dispersion_env%d2y_dx2
1191 : !$OMP PARALLEL DO SIMD DEFAULT(NONE) &
1192 8840 : !$OMP SHARED(iq, q_low, coeff, theta_scale, d2y, theta) PRIVATE(lower, p) SCHEDULE(STATIC)
1193 : DO i = 1, SIZE(theta)
1194 : lower = q_low(i)
1195 : theta(i) = 0.0_dp
1196 : IF (lower == 0) CYCLE
1197 : p = coeff(i, 1)*MERGE(1.0_dp, 0.0_dp, iq == lower) &
1198 : + coeff(i, 2)*MERGE(1.0_dp, 0.0_dp, iq == lower + 1) &
1199 : + (coeff(i, 3)*d2y(iq, lower) + coeff(i, 4)*d2y(iq, lower + 1))
1200 : theta(i) = p*theta_scale(i)
1201 : END DO
1202 : !$OMP END PARALLEL DO SIMD
1203 8840 : END SUBROUTINE build_theta
1204 :
1205 : ! **************************************************************************************************
1206 : !> \brief Accumulate one inverse-transformed channel without storing a real-space channel matrix.
1207 : !> \param iq Channel index
1208 : !> \param q_low Lower spline endpoint, using the potential-side knot convention
1209 : !> \param dispersion_env Non-local functional parameters
1210 : !> \param u Inverse FFT output
1211 : !> \param u_contract Sums against the two second-derivative columns, then the two endpoint values
1212 : ! **************************************************************************************************
1213 8720 : SUBROUTINE accumulate_potential(iq, q_low, dispersion_env, u, u_contract)
1214 : INTEGER, INTENT(IN) :: iq, q_low(:)
1215 : TYPE(qs_dispersion_type), POINTER :: dispersion_env
1216 : REAL(dp), INTENT(IN) :: u(:)
1217 : REAL(dp), INTENT(INOUT) :: u_contract(:, :)
1218 :
1219 : INTEGER :: i, lower
1220 : REAL(dp), POINTER :: d2y(:, :)
1221 :
1222 8720 : d2y => dispersion_env%d2y_dx2
1223 : !$OMP PARALLEL DO SIMD DEFAULT(NONE) &
1224 8720 : !$OMP SHARED(iq, q_low, d2y, u, u_contract) PRIVATE(lower) SCHEDULE(STATIC)
1225 : DO i = 1, SIZE(u)
1226 : lower = q_low(i)
1227 : IF (lower == 0) CYCLE
1228 : u_contract(i, 1) = u_contract(i, 1) + u(i)*d2y(iq, lower)
1229 : u_contract(i, 2) = u_contract(i, 2) + u(i)*d2y(iq, lower + 1)
1230 : IF (iq == lower) u_contract(i, 3) = u(i)
1231 : IF (iq == lower + 1) u_contract(i, 4) = u(i)
1232 : END DO
1233 : !$OMP END PARALLEL DO SIMD
1234 8720 : END SUBROUTINE accumulate_potential
1235 :
1236 : ! **************************************************************************************************
1237 : !> \brief This routine is modeled after an algorithm from "Numerical Recipes in C" by Cambridge
1238 : !> University Press, pages 96-97. It was adapted for Fortran and for the problem at hand.
1239 : !> \param x ...
1240 : !> \param d2y_dx2 ...
1241 : !> \par History
1242 : !> OpenMP added: Aug 2016 MTucker
1243 : ! **************************************************************************************************
1244 50 : SUBROUTINE initialize_spline_interpolation(x, d2y_dx2)
1245 :
1246 : REAL(dp), INTENT(in) :: x(:)
1247 : REAL(dp), INTENT(inout) :: d2y_dx2(:, :)
1248 :
1249 : INTEGER :: index, Nx, P_i
1250 : REAL(dp) :: temp1, temp2
1251 50 : REAL(dp), ALLOCATABLE :: temp_array(:), y(:)
1252 :
1253 50 : Nx = SIZE(x)
1254 :
1255 : !$OMP PARALLEL DEFAULT( NONE ) &
1256 : !$OMP SHARED( x, d2y_dx2, Nx ) &
1257 : !$OMP PRIVATE( temp_array, y &
1258 : !$OMP , index, temp1, temp2 &
1259 50 : !$OMP )
1260 :
1261 : ALLOCATE (temp_array(Nx), y(Nx))
1262 :
1263 : !$OMP DO
1264 : DO P_i = 1, Nx
1265 : !! In the Soler method, the polynomials that are interpolated are Kronecker delta functions
1266 : !! at a particular q point. So, we set all y values to 0 except the one corresponding to
1267 : !! the particular function P_i.
1268 : !! ----------------------------------------------------------------------------------------
1269 : y = 0.0_dp
1270 : y(P_i) = 1.0_dp
1271 : !! ----------------------------------------------------------------------------------------
1272 :
1273 : d2y_dx2(P_i, 1) = 0.0_dp
1274 : temp_array(1) = 0.0_dp
1275 : DO index = 2, Nx - 1
1276 : temp1 = (x(index) - x(index - 1))/(x(index + 1) - x(index - 1))
1277 : temp2 = temp1*d2y_dx2(P_i, index - 1) + 2.0_dp
1278 : d2y_dx2(P_i, index) = (temp1 - 1.0_dp)/temp2
1279 : temp_array(index) = (y(index + 1) - y(index))/(x(index + 1) - x(index)) &
1280 : - (y(index) - y(index - 1))/(x(index) - x(index - 1))
1281 : temp_array(index) = (6.0_dp*temp_array(index)/(x(index + 1) - x(index - 1)) &
1282 : - temp1*temp_array(index - 1))/temp2
1283 : END DO
1284 : d2y_dx2(P_i, Nx) = 0.0_dp
1285 : DO index = Nx - 1, 1, -1
1286 : d2y_dx2(P_i, index) = d2y_dx2(P_i, index)*d2y_dx2(P_i, index + 1) + temp_array(index)
1287 : END DO
1288 : END DO
1289 : !$OMP END DO
1290 :
1291 : DEALLOCATE (temp_array, y)
1292 : !$OMP END PARALLEL
1293 :
1294 50 : END SUBROUTINE initialize_spline_interpolation
1295 :
1296 : ! **************************************************************************************************
1297 : !> \brief Interpolate the symmetric kernel and optionally its radial derivative from packed tables.
1298 : !> \param k Reciprocal-vector length
1299 : !> \param kernel_of_k Interpolated kernel matrix
1300 : !> \param dispersion_env Radial mesh and channel count
1301 : !> \param kernel_table Packed channel pairs, radial points, and value/second derivative
1302 : !> \param dkernel_of_dk Optional derivative matrix for the kernel virial
1303 : ! **************************************************************************************************
1304 511249 : SUBROUTINE interpolate_kernel(k, kernel_of_k, dispersion_env, kernel_table, dkernel_of_dk)
1305 : REAL(dp), INTENT(IN) :: k
1306 : REAL(dp), INTENT(OUT) :: kernel_of_k(:, :)
1307 : TYPE(qs_dispersion_type), POINTER :: dispersion_env
1308 : REAL(dp), INTENT(IN) :: kernel_table(:, 0:, :)
1309 : REAL(dp), INTENT(OUT), OPTIONAL :: dkernel_of_dk(:, :)
1310 :
1311 : INTEGER :: ipair, k_i, q1_i, q2_i
1312 : LOGICAL :: on_mesh
1313 : REAL(dp) :: a, b, c, d, da, db, dc, dd, dk, dk_6, &
1314 : value
1315 :
1316 511249 : dk = dispersion_env%dk
1317 511249 : CPASSERT(k < dispersion_env%nr_points*dk)
1318 511249 : k_i = INT(k/dk)
1319 511249 : on_mesh = MOD(k, dk) == 0.0_dp
1320 511249 : a = (dk*(k_i + 1.0_dp) - k)/dk
1321 511249 : b = (k - dk*k_i)/dk
1322 511249 : c = (a**3 - a)*(dk*dk/6.0_dp)
1323 511249 : d = (b**3 - b)*(dk*dk/6.0_dp)
1324 10736229 : DO q1_i = 1, dispersion_env%nqs
1325 10736229 : !$OMP SIMD PRIVATE(ipair, value)
1326 : DO q2_i = 1, q1_i
1327 107362290 : ipair = q1_i*(q1_i - 1)/2 + q2_i
1328 107362290 : IF (on_mesh) THEN
1329 46410 : value = kernel_table(ipair, k_i, 1)
1330 : ELSE
1331 : value = a*kernel_table(ipair, k_i, 1) + b*kernel_table(ipair, k_i + 1, 1) &
1332 107315880 : + (c*kernel_table(ipair, k_i, 2) + d*kernel_table(ipair, k_i + 1, 2))
1333 : END IF
1334 107362290 : kernel_of_k(q2_i, q1_i) = value
1335 107362290 : kernel_of_k(q1_i, q2_i) = value
1336 : END DO
1337 : !$OMP END SIMD
1338 : END DO
1339 511249 : IF (PRESENT(dkernel_of_dk)) THEN
1340 405804 : dk_6 = dk/6.0_dp
1341 405804 : da = -1.0_dp/dk
1342 405804 : db = 1.0_dp/dk
1343 405804 : dc = -(3*a**2 - 1.0_dp)*dk_6
1344 405804 : dd = (3*b**2 - 1.0_dp)*dk_6
1345 8521884 : DO q1_i = 1, dispersion_env%nqs
1346 8521884 : !$OMP SIMD PRIVATE(ipair, value)
1347 : DO q2_i = 1, q1_i
1348 85218840 : ipair = q1_i*(q1_i - 1)/2 + q2_i
1349 : value = da*kernel_table(ipair, k_i, 1) + db*kernel_table(ipair, k_i + 1, 1) &
1350 85218840 : + dc*kernel_table(ipair, k_i, 2) + dd*kernel_table(ipair, k_i + 1, 2)
1351 85218840 : dkernel_of_dk(q2_i, q1_i) = value
1352 85218840 : dkernel_of_dk(q1_i, q2_i) = value
1353 : END DO
1354 : !$OMP END SIMD
1355 : END DO
1356 : END IF
1357 511249 : END SUBROUTINE interpolate_kernel
1358 :
1359 : ! **************************************************************************************************
1360 :
1361 : END MODULE qs_dispersion_nonloc
|
