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 422 : 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 422 : 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 422 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: dq0_dgradrho, dq0_drho, hpot, q0, rho, &
170 422 : theta_scale
171 422 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: drho, spline_coeff, u_contract
172 422 : 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 422 : 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 422 : CALL timeset(routineN, handle)
179 :
180 422 : CPASSERT(ASSOCIATED(rho_r))
181 422 : CPASSERT(ASSOCIATED(rho_g))
182 422 : CPASSERT(ASSOCIATED(pw_pool))
183 :
184 422 : IF (PRESENT(virial)) THEN
185 416 : 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 112 : CPASSERT(.NOT. energy_only)
191 : END IF
192 422 : IF (.NOT. energy_only) THEN
193 416 : CPASSERT(ASSOCIATED(vxc_rho))
194 : END IF
195 :
196 422 : nl_type = dispersion_env%nl_type
197 :
198 422 : b_value = dispersion_env%b_value
199 422 : beta = 0.03125_dp*(3.0_dp/(b_value**2.0_dp))**0.75_dp
200 422 : nspin = SIZE(rho_r)
201 :
202 : ! temporary arrays for FFT
203 422 : CALL pw_pool%create_pw(tmp_g)
204 422 : CALL pw_pool%create_pw(tmp_r)
205 :
206 : ! Sum the spin densities on the vdW grid before transforming or differentiating.
207 422 : CALL pw_pool%create_pw(rho_tot_g)
208 422 : CALL pw_transfer(rho_g(1), rho_tot_g)
209 442 : DO ispin = 2, nspin
210 20 : CALL pw_transfer(rho_g(ispin), tmp_g)
211 442 : CALL pw_axpy(tmp_g, rho_tot_g, 1._dp)
212 : END DO
213 422 : CALL pw_transfer(rho_tot_g, tmp_r)
214 :
215 1688 : np = SIZE(tmp_r%array)
216 422 : tmp_1d(1:np) => tmp_r%array
217 2110 : ALLOCATE (rho(np), drho(np, 3))
218 1688 : DO i = 1, 3
219 1266 : lo(i) = LBOUND(tmp_r%array, i)
220 1266 : hi(i) = UBOUND(tmp_r%array, i)
221 1688 : n(i) = hi(i) - lo(i) + 1
222 : END DO
223 : !$OMP PARALLEL DO DEFAULT(NONE) &
224 422 : !$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 1688 : DO idir = 1, 3
235 1266 : CALL pw_transfer(rho_tot_g, tmp_g)
236 1266 : CALL pw_derive(tmp_g, nd(:, idir))
237 1266 : CALL pw_transfer(tmp_g, tmp_r)
238 : !$OMP PARALLEL DO DEFAULT(NONE) &
239 1688 : !$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 422 : 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 422 : 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 1664 : 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 246 : CALL get_q0_on_grid_vdw(rho, drho, q0, dq0_drho, dq0_dgradrho, dispersion_env)
277 : CASE (vdw_nl_RVV10)
278 416 : 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 2532 : ALLOCATE (q_low(np), spline_coeff(np, 4), theta_scale(np))
285 422 : CALL prepare_splines(q0, rho, dispersion_env, q_low, spline_coeff, theta_scale)
286 9706 : ALLOCATE (thetas_g(dispersion_env%nqs))
287 422 : CALL timeset("vdW_theta_forward", handle_fft)
288 8862 : DO i = 1, dispersion_env%nqs
289 8440 : CALL build_theta(i, q_low, spline_coeff, theta_scale, dispersion_env, tmp_1d)
290 8440 : CALL pw_pool%create_pw(thetas_g(i))
291 8862 : CALL pw_transfer(tmp_r, thetas_g(i))
292 : END DO
293 422 : CALL timestop(handle_fft)
294 422 : DEALLOCATE (spline_coeff, theta_scale)
295 422 : 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 422 : sumnp = np
302 422 : CALL para_env%sum(sumnp)
303 422 : IF (use_virial) THEN
304 : ! calculates kernel contribution to stress
305 112 : CALL vdW_energy(thetas_g, dispersion_env, Ec_nl, energy_only, virial)
306 42 : SELECT CASE (nl_type)
307 : CASE (vdw_nl_RVV10)
308 290416 : 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 448 : DO idir = 1, 3
313 448 : virial%pv_xc(idir, idir) = virial%pv_xc(idir, idir) + Ec_nl
314 : END DO
315 : ELSE
316 310 : CALL vdW_energy(thetas_g, dispersion_env, Ec_nl, energy_only)
317 134 : SELECT CASE (nl_type)
318 : CASE (vdw_nl_RVV10)
319 662166 : Ec_nl = 0.5_dp*Ec_nl + beta*SUM(rho(:))*grid%vol/sumnp
320 : END SELECT
321 : END IF
322 422 : CALL para_env%sum(Ec_nl)
323 422 : IF (nl_type == vdw_nl_RVV10) Ec_nl = Ec_nl*dispersion_env%scale_rvv10
324 422 : edispersion = Ec_nl
325 :
326 422 : 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 1248 : ALLOCATE (u_contract(np, 4), hpot(np))
332 416 : !$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 416 : CALL timeset("vdW_theta_inverse", handle_fft)
342 8736 : DO i = 1, dispersion_env%nqs
343 8320 : CALL pw_transfer(thetas_g(i), tmp_r)
344 8736 : CALL accumulate_potential(i, q_low, dispersion_env, tmp_1d, u_contract)
345 : END DO
346 416 : CALL timestop(handle_fft)
347 :
348 : ! Write the local potential directly into its PW object.
349 416 : CALL pw_pool%create_pw(vxc_r)
350 416 : vxc_1d(1:np) => vxc_r%array
351 416 : IF (use_virial) THEN
352 112 : grid => tmp_g%pw_grid
353 : CALL get_potential(q0, dq0_drho, dq0_dgradrho, rho, q_low, u_contract, vxc_1d, hpot, &
354 112 : 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 304 : dispersion_env)
358 : END IF
359 416 : DEALLOCATE (u_contract, q_low, q0, dq0_drho, dq0_dgradrho)
360 170 : SELECT CASE (nl_type)
361 : CASE (vdw_nl_RVV10)
362 416 : !$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 416 : 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 416 : CALL pw_pool%create_pw(div_g)
373 1664 : DO idir = 1, 3
374 : !$OMP PARALLEL DO DEFAULT(NONE) &
375 1248 : !$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 1248 : CALL pw_transfer(tmp_r, tmp_g)
386 1248 : CALL pw_derive(tmp_g, nd(:, idir))
387 1664 : IF (idir == 1) THEN
388 416 : CALL pw_transfer(tmp_g, div_g)
389 : ELSE
390 832 : CALL pw_axpy(tmp_g, div_g, 1._dp)
391 : END IF
392 : END DO
393 416 : CALL pw_transfer(div_g, tmp_r)
394 416 : CALL pw_pool%give_back_pw(div_g)
395 416 : CALL pw_axpy(tmp_r, vxc_r, -1._dp)
396 416 : CALL pw_transfer(vxc_r, tmp_g)
397 416 : CALL pw_pool%give_back_pw(vxc_r)
398 416 : CALL xc_pw_pool%create_pw(vxc_r)
399 416 : CALL xc_pw_pool%create_pw(vxc_g)
400 416 : CALL pw_transfer(tmp_g, vxc_g)
401 416 : CALL pw_transfer(vxc_g, vxc_r)
402 852 : DO ispin = 1, nspin
403 852 : CALL pw_axpy(vxc_r, vxc_rho(ispin), 1._dp)
404 : END DO
405 416 : CALL xc_pw_pool%give_back_pw(vxc_r)
406 832 : CALL xc_pw_pool%give_back_pw(vxc_g)
407 : END IF
408 :
409 : NULLIFY (tmp_1d)
410 :
411 8862 : DO i = 1, dispersion_env%nqs
412 8862 : CALL pw_pool%give_back_pw(thetas_g(i))
413 : END DO
414 422 : CALL pw_pool%give_back_pw(tmp_r)
415 422 : CALL pw_pool%give_back_pw(tmp_g)
416 :
417 422 : DEALLOCATE (rho, drho, thetas_g)
418 :
419 422 : CALL timestop(handle)
420 :
421 1266 : 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 422 : 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 :
445 : INTEGER :: handle, ig, iq, l, m, nl_type, nqs, &
446 : q1_i, q2_i
447 : LOGICAL :: use_virial
448 : REAL(KIND=dp) :: g, g2, g2_last, g_multiplier, gm
449 422 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: theta_im, theta_re, u_im, u_re
450 422 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dkernel_of_dk, kernel_of_k
451 : REAL(KIND=dp), DIMENSION(3, 3) :: virial_thread
452 : TYPE(pw_grid_type), POINTER :: grid
453 :
454 422 : CALL timeset(routineN, handle)
455 422 : nqs = dispersion_env%nqs
456 :
457 422 : use_virial = PRESENT(virial)
458 422 : virial_thread(:, :) = 0.0_dp ! always initialize to avoid floating point exceptions in OMP REDUCTION
459 :
460 422 : vdW_xc_energy = 0._dp
461 422 : grid => thetas_g(1)%pw_grid
462 :
463 422 : IF (grid%grid_span == HALFSPACE) THEN
464 : g_multiplier = 2._dp
465 : ELSE
466 422 : g_multiplier = 1._dp
467 : END IF
468 :
469 422 : nl_type = dispersion_env%nl_type
470 :
471 : !$OMP PARALLEL DEFAULT(NONE) &
472 : !$OMP SHARED(nqs, energy_only, grid, dispersion_env, use_virial, thetas_g, &
473 : !$OMP g_multiplier, nl_type) &
474 : !$OMP PRIVATE(g2_last, kernel_of_k, dkernel_of_dk, theta_re, theta_im, &
475 : !$OMP g2, g, iq, q2_i, u_re, u_im, q1_i, gm, l, m) &
476 422 : !$OMP REDUCTION(+:vdW_xc_energy, virial_thread)
477 :
478 : g2_last = HUGE(0._dp)
479 :
480 : ALLOCATE (kernel_of_k(nqs, nqs))
481 : IF (use_virial) ALLOCATE (dkernel_of_dk(nqs, nqs))
482 : ALLOCATE (theta_re(nqs), theta_im(nqs), u_re(nqs), u_im(nqs))
483 :
484 : !$OMP DO
485 : DO ig = 1, grid%ngpts_cut_local
486 : g2 = grid%gsq(ig)
487 : IF (ABS(g2 - g2_last) > 1.e-10) THEN
488 : g2_last = g2
489 : g = SQRT(g2)
490 : IF (use_virial) THEN
491 : CALL interpolate_kernel(g, kernel_of_k, dispersion_env, dispersion_env%kernel_table, dkernel_of_dk)
492 : ELSE
493 : CALL interpolate_kernel(g, kernel_of_k, dispersion_env, dispersion_env%kernel_table)
494 : END IF
495 : END IF
496 : ! Save all inputs before in-place output. Vectorize over output channels;
497 : ! the input-channel accumulation order remains unchanged for every output.
498 : DO iq = 1, nqs
499 : theta_re(iq) = REAL(thetas_g(iq)%array(ig), KIND=dp)
500 : theta_im(iq) = AIMAG(thetas_g(iq)%array(ig))
501 : END DO
502 : u_re(:) = 0.0_dp
503 : u_im(:) = 0.0_dp
504 : DO q1_i = 1, nqs
505 : !$OMP SIMD
506 : DO q2_i = 1, nqs
507 : u_re(q2_i) = u_re(q2_i) + kernel_of_k(q2_i, q1_i)*theta_re(q1_i)
508 : u_im(q2_i) = u_im(q2_i) + kernel_of_k(q2_i, q1_i)*theta_im(q1_i)
509 : END DO
510 : !$OMP END SIMD
511 : END DO
512 : DO q2_i = 1, nqs
513 : IF (ig < grid%first_gne0) THEN
514 : vdW_xc_energy = vdW_xc_energy + (u_re(q2_i)*theta_re(q2_i) + u_im(q2_i)*theta_im(q2_i))
515 : ELSE
516 : vdW_xc_energy = vdW_xc_energy &
517 : + g_multiplier*(u_re(q2_i)*theta_re(q2_i) + u_im(q2_i)*theta_im(q2_i))
518 : END IF
519 : IF (.NOT. energy_only) thetas_g(q2_i)%array(ig) = CMPLX(u_re(q2_i), u_im(q2_i), KIND=dp)
520 : END DO
521 :
522 : IF (use_virial .AND. ig >= grid%first_gne0) THEN
523 : ! Reduce over all channel pairs before assembling the stress tensor.
524 : gm = 0.0_dp
525 : DO q2_i = 1, nqs
526 : !$OMP SIMD REDUCTION(+:gm)
527 : DO q1_i = 1, nqs
528 : gm = gm + dkernel_of_dk(q1_i, q2_i) &
529 : *(theta_re(q1_i)*theta_re(q2_i) + theta_im(q1_i)*theta_im(q2_i))
530 : END DO
531 : !$OMP END SIMD
532 : END DO
533 : gm = 0.5_dp*g_multiplier*grid%vol*gm
534 : IF (nl_type == vdw_nl_RVV10) gm = 0.5_dp*gm
535 : DO l = 1, 3
536 : DO m = 1, l
537 : virial_thread(l, m) = virial_thread(l, m) - gm*(grid%g(l, ig)*grid%g(m, ig))/g
538 : END DO
539 : END DO
540 : END IF
541 : END DO
542 : !$OMP END DO
543 :
544 : DEALLOCATE (theta_re, theta_im, u_re, u_im, kernel_of_k)
545 : IF (use_virial) DEALLOCATE (dkernel_of_dk)
546 :
547 : !$OMP END PARALLEL
548 :
549 422 : vdW_xc_energy = vdW_xc_energy*grid%vol*0.5_dp
550 :
551 422 : IF (use_virial) THEN
552 448 : DO l = 1, 3
553 672 : DO m = 1, (l - 1)
554 336 : virial%pv_xc(l, m) = virial%pv_xc(l, m) + virial_thread(l, m)
555 672 : virial%pv_xc(m, l) = virial%pv_xc(l, m)
556 : END DO
557 336 : m = l
558 448 : virial%pv_xc(l, m) = virial%pv_xc(l, m) + virial_thread(l, m)
559 : END DO
560 : END IF
561 :
562 422 : CALL timestop(handle)
563 :
564 422 : END SUBROUTINE vdW_energy
565 :
566 : ! **************************************************************************************************
567 : !> \brief This routine finds the non-local correlation contribution to the potential
568 : !> (i.e. the derivative of the non-local piece of the energy with respect to
569 : !> density) given in SOLER equation 10. The u_alpha(k) functions were found
570 : !> while calculating the energy. Their spline contractions are accumulated during the inverse FFTs.
571 : !> Most of the required derivatives were calculated in the "get_q0_on_grid"
572 : !> routine, but the derivative of the interpolation polynomials, P_alpha(q),
573 : !> (SOLER equation 3) with respect to q is interpolated here, along with the
574 : !> polynomials themselves.
575 : !> \param q0 ...
576 : !> \param dq0_drho ...
577 : !> \param dq0_dgradrho ...
578 : !> \param total_rho ...
579 : !> \param q_low_grid Lower spline endpoint at each grid point
580 : !> \param u_contract Second-derivative contractions and endpoint values
581 : !> \param potential ...
582 : !> \param h_prefactor ...
583 : !> \param dispersion_env ...
584 : !> \param drho ...
585 : !> \param dvol ...
586 : !> \param virial ...
587 : !> \par History
588 : !> OpenMP added: Aug 2016 MTucker
589 : ! **************************************************************************************************
590 416 : SUBROUTINE get_potential(q0, dq0_drho, dq0_dgradrho, total_rho, q_low_grid, u_contract, potential, h_prefactor, &
591 416 : dispersion_env, drho, dvol, virial)
592 :
593 : REAL(dp), DIMENSION(:), INTENT(in) :: q0, dq0_drho, dq0_dgradrho, total_rho
594 : INTEGER, DIMENSION(:), INTENT(IN) :: q_low_grid
595 : REAL(dp), DIMENSION(:, :), INTENT(in) :: u_contract
596 : REAL(dp), DIMENSION(:), INTENT(out) :: potential, h_prefactor
597 : TYPE(qs_dispersion_type), POINTER :: dispersion_env
598 : REAL(dp), DIMENSION(:, :), INTENT(in), OPTIONAL :: drho
599 : REAL(dp), INTENT(IN), OPTIONAL :: dvol
600 : TYPE(virial_type), OPTIONAL, POINTER :: virial
601 :
602 : CHARACTER(len=*), PARAMETER :: routineN = 'get_potential'
603 :
604 : INTEGER :: handle, i_grid, l, m, nl_type, nqs, &
605 : q_hi, q_low
606 : LOGICAL :: use_virial
607 : REAL(dp) :: a, b, b_value, c, const, d, dq, dq_6, e, &
608 : f, prefactor, tmp_1_2, tmp_1_4, &
609 : tmp_3_4, u_dp_dq0, u_p
610 : REAL(dp), DIMENSION(3, 3) :: virial_thread
611 416 : REAL(dp), DIMENSION(:), POINTER :: q_mesh
612 :
613 416 : CALL timeset(routineN, handle)
614 :
615 416 : use_virial = PRESENT(virial)
616 416 : CPASSERT(.NOT. use_virial .OR. PRESENT(drho))
617 416 : CPASSERT(.NOT. use_virial .OR. PRESENT(dvol))
618 :
619 416 : virial_thread(:, :) = 0.0_dp ! always initialize to avoid floating point exceptions in OMP REDUCTION
620 416 : b_value = dispersion_env%b_value
621 416 : const = 1.0_dp/(3.0_dp*b_value**(3.0_dp/2.0_dp)*pi**(5.0_dp/4.0_dp))
622 :
623 416 : q_mesh => dispersion_env%q_mesh
624 416 : nqs = dispersion_env%nqs
625 416 : nl_type = dispersion_env%nl_type
626 :
627 : !$OMP PARALLEL DEFAULT(NONE) &
628 : !$OMP SHARED(nqs, u_contract, q_low_grid, q_mesh, q0, nl_type, potential, h_prefactor, &
629 : !$OMP dq0_drho, dq0_dgradrho, total_rho, const, use_virial, drho, dvol, virial) &
630 : !$OMP PRIVATE(q_low, q_hi, dq, dq_6, A, b, c, d, e, f, u_p, u_dp_dq0, &
631 : !$OMP prefactor, l, m, tmp_1_2, tmp_1_4, tmp_3_4) &
632 416 : !$OMP REDUCTION(+:virial_thread)
633 :
634 : !$OMP DO
635 : DO i_grid = 1, SIZE(q0)
636 : potential(i_grid) = 0.0_dp
637 : h_prefactor(i_grid) = 0.0_dp
638 : IF (nl_type == vdw_nl_RVV10 .AND. total_rho(i_grid) <= epsr) CYCLE
639 : q_low = q_low_grid(i_grid)
640 : q_hi = q_low + 1
641 :
642 : dq = q_mesh(q_hi) - q_mesh(q_low)
643 : dq_6 = dq/6.0_dp
644 :
645 : a = (q_mesh(q_hi) - q0(i_grid))/dq
646 : b = (q0(i_grid) - q_mesh(q_low))/dq
647 : c = (a**3 - a)*dq*dq_6
648 : d = (b**3 - b)*dq*dq_6
649 : e = (3.0_dp*a**2 - 1.0_dp)*dq_6
650 : f = (3.0_dp*b**2 - 1.0_dp)*dq_6
651 :
652 : u_p = a*u_contract(i_grid, 3) + b*u_contract(i_grid, 4) &
653 : + c*u_contract(i_grid, 1) + d*u_contract(i_grid, 2)
654 : u_dp_dq0 = (u_contract(i_grid, 4) - u_contract(i_grid, 3))/dq &
655 : - e*u_contract(i_grid, 1) + f*u_contract(i_grid, 2)
656 :
657 : !! The first term in equation 13 of SOLER
658 : SELECT CASE (nl_type)
659 : CASE DEFAULT
660 : CPABORT("Unknown vdW-DF functional")
661 : CASE (vdw_nl_DRSLL, vdw_nl_LMKLL)
662 : potential(i_grid) = u_p + u_dp_dq0*dq0_drho(i_grid)
663 : prefactor = u_dp_dq0*dq0_dgradrho(i_grid)
664 : CASE (vdw_nl_RVV10)
665 : tmp_1_2 = SQRT(total_rho(i_grid))
666 : tmp_1_4 = SQRT(tmp_1_2)
667 : tmp_3_4 = tmp_1_4*tmp_1_4*tmp_1_4
668 : potential(i_grid) = const*0.75_dp/tmp_1_4*u_p + const*tmp_3_4*u_dp_dq0*dq0_drho(i_grid)
669 : prefactor = const*tmp_3_4*u_dp_dq0*dq0_dgradrho(i_grid)
670 : END SELECT
671 : IF (q0(i_grid) /= q_mesh(nqs)) THEN
672 : h_prefactor(i_grid) = prefactor
673 : END IF
674 :
675 : ! The saturation guard applies only to h_prefactor, not to the virial.
676 : IF (use_virial .AND. ABS(prefactor) > 0.0_dp) THEN
677 : IF (nl_type == vdw_nl_RVV10) prefactor = 0.5_dp*prefactor
678 : prefactor = prefactor*dvol
679 : DO l = 1, 3
680 : DO m = 1, l
681 : virial_thread(l, m) = virial_thread(l, m) - prefactor*drho(i_grid, l)*drho(i_grid, m)
682 : END DO
683 : END DO
684 : END IF
685 : END DO ! i_grid = 1, SIZE(q0)
686 : !$OMP END DO
687 :
688 : !$OMP END PARALLEL
689 :
690 416 : IF (use_virial) THEN
691 448 : DO l = 1, 3
692 672 : DO m = 1, (l - 1)
693 336 : virial%pv_xc(l, m) = virial%pv_xc(l, m) + virial_thread(l, m)
694 672 : virial%pv_xc(m, l) = virial%pv_xc(l, m)
695 : END DO
696 336 : m = l
697 448 : virial%pv_xc(l, m) = virial%pv_xc(l, m) + virial_thread(l, m)
698 : END DO
699 : END IF
700 :
701 416 : CALL timestop(handle)
702 416 : END SUBROUTINE get_potential
703 :
704 : ! **************************************************************************************************
705 : !> \brief calculates exponent = sum(from i=1 to hi, ((alpha)**i)/i) ) without <<< calling power >>>
706 : !> \param hi = upper index for sum
707 : !> \param alpha ...
708 : !> \param exponent = output value
709 : !> \par History
710 : !> Created: MTucker, Aug 2016
711 : ! **************************************************************************************************
712 72326 : ELEMENTAL SUBROUTINE calculate_exponent(hi, alpha, exponent)
713 : INTEGER, INTENT(in) :: hi
714 : REAL(dp), INTENT(in) :: alpha
715 : REAL(dp), INTENT(out) :: exponent
716 :
717 : INTEGER :: i
718 : REAL(dp) :: multiplier
719 :
720 72326 : multiplier = alpha
721 72326 : exponent = alpha
722 :
723 867912 : DO i = 2, hi
724 795586 : multiplier = multiplier*alpha
725 867912 : exponent = exponent + (multiplier/i)
726 : END DO
727 72326 : END SUBROUTINE calculate_exponent
728 :
729 : ! **************************************************************************************************
730 : !> \brief calculate exponent = sum(from i=1 to hi, ((alpha)**i)/i) ) without calling power
731 : !> also calculates derivative using similar series
732 : !> \param hi = upper index for sum
733 : !> \param alpha ...
734 : !> \param exponent = output value
735 : !> \param derivative ...
736 : !> \par History
737 : !> Created: MTucker, Aug 2016
738 : ! **************************************************************************************************
739 5076074 : ELEMENTAL SUBROUTINE calculate_exponent_derivative(hi, alpha, exponent, derivative)
740 : INTEGER, INTENT(in) :: hi
741 : REAL(dp), INTENT(in) :: alpha
742 : REAL(dp), INTENT(out) :: exponent, derivative
743 :
744 : INTEGER :: i
745 : REAL(dp) :: multiplier
746 :
747 5076074 : derivative = 0.0d0
748 5076074 : multiplier = 1.0d0
749 5076074 : exponent = 0.0d0
750 :
751 65988962 : DO i = 1, hi
752 60912888 : derivative = derivative + multiplier
753 60912888 : multiplier = multiplier*alpha
754 65988962 : exponent = exponent + (multiplier/i)
755 : END DO
756 5076074 : END SUBROUTINE calculate_exponent_derivative
757 :
758 : !! This routine first calculates the q value defined in (DION equations 11 and 12), then
759 : !! saturates it according to (SOLER equation 5).
760 : ! **************************************************************************************************
761 : !> \brief This routine first calculates the q value defined in (DION equations 11 and 12), then
762 : !> saturates it according to (SOLER equation 5).
763 : !> \param total_rho ...
764 : !> \param gradient_rho ...
765 : !> \param q0 ...
766 : !> \param dq0_drho ...
767 : !> \param dq0_dgradrho ...
768 : !> \param dispersion_env ...
769 : ! **************************************************************************************************
770 246 : SUBROUTINE get_q0_on_grid_vdw(total_rho, gradient_rho, q0, dq0_drho, dq0_dgradrho, dispersion_env)
771 : !!
772 : !! more specifically it calculates the following
773 : !!
774 : !! q0(ir) = q0 as defined above
775 : !! dq0_drho(ir) = total_rho * d q0 /d rho
776 : !! dq0_dgradrho = total_rho / |gradient_rho| * d q0 / d |gradient_rho|
777 : !!
778 : REAL(dp), INTENT(IN) :: total_rho(:), gradient_rho(:, :)
779 : REAL(dp), INTENT(OUT) :: q0(:), dq0_drho(:), dq0_dgradrho(:)
780 : TYPE(qs_dispersion_type), POINTER :: dispersion_env
781 :
782 : INTEGER, PARAMETER :: m_cut = 12
783 : REAL(dp), PARAMETER :: LDA_A = 0.031091_dp, LDA_a1 = 0.2137_dp, LDA_b1 = 7.5957_dp, &
784 : LDA_b2 = 3.5876_dp, LDA_b3 = 1.6382_dp, LDA_b4 = 0.49294_dp
785 :
786 : INTEGER :: i_grid
787 : REAL(dp) :: dq0_dq, exponent, gradient_correction, &
788 : kF, LDA_1, LDA_2, q, q__q_cut, q_cut, &
789 : q_min, r_s, sqrt_r_s, Z_ab
790 :
791 246 : q_cut = dispersion_env%q_cut
792 246 : q_min = dispersion_env%q_min
793 246 : SELECT CASE (dispersion_env%nl_type)
794 : CASE DEFAULT
795 0 : CPABORT("Unknown vdW-DF functional")
796 : CASE (vdw_nl_DRSLL)
797 48 : Z_ab = -0.8491_dp
798 : CASE (vdw_nl_LMKLL)
799 246 : Z_ab = -1.887_dp
800 : END SELECT
801 :
802 : !$OMP PARALLEL DO DEFAULT(NONE) &
803 : !$OMP SHARED(total_rho, gradient_rho, q0, dq0_drho, dq0_dgradrho, q_cut, q_min, Z_ab) &
804 : !$OMP PRIVATE(dq0_dq, exponent, gradient_correction, kF, LDA_1, LDA_2, q, q__q_cut, r_s, sqrt_r_s) &
805 246 : !$OMP SCHEDULE(STATIC)
806 : DO i_grid = 1, SIZE(total_rho)
807 : q0(i_grid) = q_cut
808 : dq0_drho(i_grid) = 0.0_dp
809 : dq0_dgradrho(i_grid) = 0.0_dp
810 :
811 : !! This prevents numerical problems. If the charge density is negative (an
812 : !! unphysical situation), we simply treat it as very small. In that case,
813 : !! q0 will be very large and will be saturated. For a saturated q0 the derivative
814 : !! dq0_dq will be 0 so we set q0 = q_cut and dq0_drho = dq0_dgradrho = 0 and go on
815 : !! to the next point.
816 : !! ------------------------------------------------------------------------------------
817 : IF (total_rho(i_grid) < epsr) CYCLE
818 : !! ------------------------------------------------------------------------------------
819 : !! Calculate some intermediate values needed to find q
820 : !! ------------------------------------------------------------------------------------
821 : kF = (3.0_dp*pi*pi*total_rho(i_grid))**(1.0_dp/3.0_dp)
822 : r_s = (3.0_dp/(4.0_dp*pi*total_rho(i_grid)))**(1.0_dp/3.0_dp)
823 : sqrt_r_s = SQRT(r_s)
824 :
825 : gradient_correction = -Z_ab/(36.0_dp*kF*total_rho(i_grid)**2) &
826 : *(gradient_rho(i_grid, 1)**2 + gradient_rho(i_grid, 2)**2 + gradient_rho(i_grid, 3)**2)
827 :
828 : LDA_1 = 8.0_dp*pi/3.0_dp*(LDA_A*(1.0_dp + LDA_a1*r_s))
829 : 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)
830 : !! ---------------------------------------------------------------
831 : !! This is the q value defined in equations 11 and 12 of DION
832 : !! ---------------------------------------------------------------
833 : q = kF + LDA_1*LOG(1.0_dp + 1.0_dp/LDA_2) + gradient_correction
834 : !! ---------------------------------------------------------------
835 : !! Here, we calculate q0 by saturating q according to equation 5 of SOLER. Also, we find
836 : !! the derivative dq0_dq needed for the derivatives dq0_drho and dq0_dgradrh0 discussed below.
837 : !! ---------------------------------------------------------------------------------------
838 : q__q_cut = q/q_cut
839 : CALL calculate_exponent_derivative(m_cut, q__q_cut, exponent, dq0_dq)
840 : q0(i_grid) = q_cut*(1.0_dp - EXP(-exponent))
841 : dq0_dq = dq0_dq*EXP(-exponent)
842 : !! ---------------------------------------------------------------------------------------
843 : !! This is to handle a case with q0 too small. We simply set it to the smallest q value in
844 : !! out q_mesh. Hopefully this doesn't get used often (ever)
845 : !! ---------------------------------------------------------------------------------------
846 : IF (q0(i_grid) < q_min) THEN
847 : q0(i_grid) = q_min
848 : END IF
849 : !! ---------------------------------------------------------------------------------------
850 : !! Here we find derivatives. These are actually the density times the derivative of q0 with respect
851 : !! to rho and gradient_rho. The density factor comes in since we are really differentiating
852 : !! theta = (rho)*P(q0) with respect to density (or its gradient) which will be
853 : !! dtheta_drho = P(q0) + dP_dq0 * [rho * dq0_dq * dq_drho] and
854 : !! dtheta_dgradient_rho = dP_dq0 * [rho * dq0_dq * dq_dgradient_rho]
855 : !! The parts in square brackets are what is calculated here. The dP_dq0 term will be interpolated
856 : !! later. There should actually be a factor of the magnitude of the gradient in the gradient_rho derivative
857 : !! but that cancels out when we differentiate the magnitude of the gradient with respect to a particular
858 : !! component.
859 : !! ------------------------------------------------------------------------------------------------
860 :
861 : dq0_drho(i_grid) = dq0_dq*(kF/3.0_dp - 7.0_dp/3.0_dp*gradient_correction &
862 : - 8.0_dp*pi/9.0_dp*LDA_A*LDA_a1*r_s*LOG(1.0_dp + 1.0_dp/LDA_2) &
863 : + LDA_1/(LDA_2*(1.0_dp + LDA_2)) &
864 : *(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 &
865 : + 2.0_dp*LDA_b4/3.0_dp*r_s**2)))
866 :
867 : dq0_dgradrho(i_grid) = total_rho(i_grid)*dq0_dq*2.0_dp*(-Z_ab)/(36.0_dp*kF*total_rho(i_grid)**2)
868 :
869 : END DO
870 : !$OMP END PARALLEL DO
871 :
872 246 : END SUBROUTINE get_q0_on_grid_vdw
873 :
874 : ! **************************************************************************************************
875 : !> \brief ...
876 : !> \param total_rho ...
877 : !> \param gradient_rho ...
878 : !> \param q0 ...
879 : !> \param dq0_drho ...
880 : !> \param dq0_dgradrho ...
881 : !> \param dispersion_env ...
882 : ! **************************************************************************************************
883 170 : SUBROUTINE get_q0_on_grid_rvv10(total_rho, gradient_rho, q0, dq0_drho, dq0_dgradrho, dispersion_env)
884 : !!
885 : !! more specifically it calculates the following
886 : !!
887 : !! q0(ir) = q0 as defined above
888 : !! dq0_drho(ir) = total_rho * d q0 /d rho
889 : !! dq0_dgradrho = total_rho / |gradient_rho| * d q0 / d |gradient_rho|
890 : !!
891 : REAL(dp), INTENT(IN) :: total_rho(:), gradient_rho(:, :)
892 : REAL(dp), INTENT(OUT) :: q0(:), dq0_drho(:), dq0_dgradrho(:)
893 : TYPE(qs_dispersion_type), POINTER :: dispersion_env
894 :
895 : INTEGER, PARAMETER :: m_cut = 12
896 :
897 : INTEGER :: i_grid
898 : REAL(dp) :: b_value, C_value, dk_dn, dq0_dq, dw0_dn, &
899 : exponent, gmod2, k, mod_grad, q, &
900 : q__q_cut, q_cut, q_min, w0, wg2, wp2
901 :
902 170 : q_cut = dispersion_env%q_cut
903 170 : q_min = dispersion_env%q_min
904 170 : b_value = dispersion_env%b_value
905 170 : C_value = dispersion_env%c_value
906 :
907 : !$OMP PARALLEL DO DEFAULT(NONE) &
908 : !$OMP SHARED(total_rho, gradient_rho, q0, dq0_drho, dq0_dgradrho, q_cut, q_min, b_value, C_value) &
909 : !$OMP PRIVATE(dk_dn, dq0_dq, dw0_dn, exponent, gmod2, k, mod_grad, q, q__q_cut, w0, wg2, wp2) &
910 170 : !$OMP SCHEDULE(STATIC)
911 : DO i_grid = 1, SIZE(total_rho)
912 : q0(i_grid) = q_cut
913 : dq0_drho(i_grid) = 0.0_dp
914 : dq0_dgradrho(i_grid) = 0.0_dp
915 :
916 : gmod2 = gradient_rho(i_grid, 1)**2 + gradient_rho(i_grid, 2)**2 + gradient_rho(i_grid, 3)**2
917 :
918 : !if (total_rho(i_grid) > epsr .and. gmod2 > epsr) cycle
919 : IF (total_rho(i_grid) > epsr) THEN
920 :
921 : !! Calculate some intermediate values needed to find q
922 : !! ------------------------------------------------------------------------------------
923 : mod_grad = SQRT(gmod2)
924 :
925 : wp2 = 16.0_dp*pi*total_rho(i_grid)
926 : wg2 = 4_dp*C_value*(mod_grad/total_rho(i_grid))**4
927 :
928 : k = b_value*3.0_dp*pi*((total_rho(i_grid)/(9.0_dp*pi))**(1.0_dp/6.0_dp))
929 : w0 = SQRT(wg2 + wp2/3.0_dp)
930 :
931 : q = w0/k
932 :
933 : !! Here, we calculate q0 by saturating q according
934 : !! ---------------------------------------------------------------------------------------
935 : q__q_cut = q/q_cut
936 : CALL calculate_exponent_derivative(m_cut, q__q_cut, exponent, dq0_dq)
937 : q0(i_grid) = q_cut*(1.0_dp - EXP(-exponent))
938 : dq0_dq = dq0_dq*EXP(-exponent)
939 :
940 : !! ---------------------------------------------------------------------------------------
941 : IF (q0(i_grid) < q_min) THEN
942 : q0(i_grid) = q_min
943 : END IF
944 :
945 : !!---------------------------------Final values---------------------------------
946 : dw0_dn = 1.0_dp/(2.0_dp*w0)*(16.0_dp/3.0_dp*pi - 4.0_dp*wg2/total_rho(i_grid))
947 : dk_dn = k/(6.0_dp*total_rho(i_grid))
948 :
949 : dq0_drho(i_grid) = dq0_dq*1.0_dp/(k**2.0)*(dw0_dn*k - dk_dn*w0)
950 : ! wg2 is proportional to |gradient_rho|**4, so this limit is zero.
951 : IF (gmod2 > 0.0_dp) THEN
952 : dq0_dgradrho(i_grid) = dq0_dq*1.0_dp/(2.0_dp*k*w0)*4.0_dp*wg2/gmod2
953 : END IF
954 : END IF
955 :
956 : END DO
957 : !$OMP END PARALLEL DO
958 :
959 170 : END SUBROUTINE get_q0_on_grid_rvv10
960 :
961 : ! **************************************************************************************************
962 : !> \brief ...
963 : !> \param total_rho ...
964 : !> \param gradient_rho ...
965 : !> \param q0 ...
966 : !> \param dispersion_env ...
967 : ! **************************************************************************************************
968 0 : SUBROUTINE get_q0_on_grid_eo_vdw(total_rho, gradient_rho, q0, dispersion_env)
969 :
970 : REAL(dp), INTENT(IN) :: total_rho(:), gradient_rho(:, :)
971 : REAL(dp), INTENT(OUT) :: q0(:)
972 : TYPE(qs_dispersion_type), POINTER :: dispersion_env
973 :
974 : INTEGER, PARAMETER :: m_cut = 12
975 : REAL(dp), PARAMETER :: LDA_A = 0.031091_dp, LDA_a1 = 0.2137_dp, LDA_b1 = 7.5957_dp, &
976 : LDA_b2 = 3.5876_dp, LDA_b3 = 1.6382_dp, LDA_b4 = 0.49294_dp
977 :
978 : INTEGER :: i_grid
979 : REAL(dp) :: exponent, gradient_correction, kF, &
980 : LDA_1, LDA_2, q, q__q_cut, q_cut, &
981 : q_min, r_s, sqrt_r_s, Z_ab
982 :
983 0 : q_cut = dispersion_env%q_cut
984 0 : q_min = dispersion_env%q_min
985 0 : SELECT CASE (dispersion_env%nl_type)
986 : CASE DEFAULT
987 0 : CPABORT("Unknown vdW-DF functional")
988 : CASE (vdw_nl_DRSLL)
989 0 : Z_ab = -0.8491_dp
990 : CASE (vdw_nl_LMKLL)
991 0 : Z_ab = -1.887_dp
992 : END SELECT
993 :
994 : !$OMP PARALLEL DO DEFAULT(NONE) &
995 : !$OMP SHARED(total_rho, gradient_rho, q0, q_cut, q_min, Z_ab) &
996 : !$OMP PRIVATE(exponent, gradient_correction, kF, LDA_1, LDA_2, q, q__q_cut, r_s, sqrt_r_s) &
997 0 : !$OMP SCHEDULE(STATIC)
998 : DO i_grid = 1, SIZE(total_rho)
999 : q0(i_grid) = q_cut
1000 : !! This prevents numerical problems. If the charge density is negative (an
1001 : !! unphysical situation), we simply treat it as very small. In that case,
1002 : !! q0 will be very large and will be saturated. For a saturated q0 the derivative
1003 : !! dq0_dq will be 0 so we set q0 = q_cut and dq0_drho = dq0_dgradrho = 0 and go on
1004 : !! to the next point.
1005 : !! ------------------------------------------------------------------------------------
1006 : IF (total_rho(i_grid) < epsr) CYCLE
1007 : !! ------------------------------------------------------------------------------------
1008 : !! Calculate some intermediate values needed to find q
1009 : !! ------------------------------------------------------------------------------------
1010 : kF = (3.0_dp*pi*pi*total_rho(i_grid))**(1.0_dp/3.0_dp)
1011 : r_s = (3.0_dp/(4.0_dp*pi*total_rho(i_grid)))**(1.0_dp/3.0_dp)
1012 : sqrt_r_s = SQRT(r_s)
1013 :
1014 : gradient_correction = -Z_ab/(36.0_dp*kF*total_rho(i_grid)**2) &
1015 : *(gradient_rho(i_grid, 1)**2 + gradient_rho(i_grid, 2)**2 + gradient_rho(i_grid, 3)**2)
1016 :
1017 : LDA_1 = 8.0_dp*pi/3.0_dp*(LDA_A*(1.0_dp + LDA_a1*r_s))
1018 : 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)
1019 : !! ------------------------------------------------------------------------------------
1020 : !! This is the q value defined in equations 11 and 12 of DION
1021 : !! ---------------------------------------------------------------
1022 : q = kF + LDA_1*LOG(1.0_dp + 1.0_dp/LDA_2) + gradient_correction
1023 :
1024 : !! ---------------------------------------------------------------
1025 : !! Here, we calculate q0 by saturating q according to equation 5 of SOLER. Also, we find
1026 : !! the derivative dq0_dq needed for the derivatives dq0_drho and dq0_dgradrh0 discussed below.
1027 : !! ---------------------------------------------------------------------------------------
1028 : q__q_cut = q/q_cut
1029 : CALL calculate_exponent(m_cut, q__q_cut, exponent)
1030 : q0(i_grid) = q_cut*(1.0_dp - EXP(-exponent))
1031 :
1032 : !! ---------------------------------------------------------------------------------------
1033 : !! This is to handle a case with q0 too small. We simply set it to the smallest q value in
1034 : !! out q_mesh. Hopefully this doesn't get used often (ever)
1035 : !! ---------------------------------------------------------------------------------------
1036 : IF (q0(i_grid) < q_min) THEN
1037 : q0(i_grid) = q_min
1038 : END IF
1039 : END DO
1040 : !$OMP END PARALLEL DO
1041 :
1042 0 : END SUBROUTINE get_q0_on_grid_eo_vdw
1043 :
1044 : ! **************************************************************************************************
1045 : !> \brief ...
1046 : !> \param total_rho ...
1047 : !> \param gradient_rho ...
1048 : !> \param q0 ...
1049 : !> \param dispersion_env ...
1050 : ! **************************************************************************************************
1051 6 : SUBROUTINE get_q0_on_grid_eo_rvv10(total_rho, gradient_rho, q0, dispersion_env)
1052 :
1053 : REAL(dp), INTENT(IN) :: total_rho(:), gradient_rho(:, :)
1054 : REAL(dp), INTENT(OUT) :: q0(:)
1055 : TYPE(qs_dispersion_type), POINTER :: dispersion_env
1056 :
1057 : INTEGER, PARAMETER :: m_cut = 12
1058 :
1059 : INTEGER :: i_grid
1060 : REAL(dp) :: b_value, C_value, exponent, gmod2, k, q, &
1061 : q__q_cut, q_cut, q_min, w0, wg2, wp2
1062 :
1063 6 : q_cut = dispersion_env%q_cut
1064 6 : q_min = dispersion_env%q_min
1065 6 : b_value = dispersion_env%b_value
1066 6 : C_value = dispersion_env%c_value
1067 :
1068 : !$OMP PARALLEL DO DEFAULT(NONE) &
1069 : !$OMP SHARED(total_rho, gradient_rho, q0, q_cut, q_min, b_value, C_value) &
1070 : !$OMP PRIVATE(exponent, gmod2, k, q, q__q_cut, w0, wg2, wp2) &
1071 6 : !$OMP SCHEDULE(STATIC)
1072 : DO i_grid = 1, SIZE(total_rho)
1073 : q0(i_grid) = q_cut
1074 :
1075 : gmod2 = gradient_rho(i_grid, 1)**2 + gradient_rho(i_grid, 2)**2 + gradient_rho(i_grid, 3)**2
1076 :
1077 : !if (total_rho(i_grid) > epsr .and. gmod2 > epsr) cycle
1078 : IF (total_rho(i_grid) > epsr) THEN
1079 :
1080 : !! Calculate some intermediate values needed to find q
1081 : !! ------------------------------------------------------------------------------------
1082 : wp2 = 16.0_dp*pi*total_rho(i_grid)
1083 : wg2 = 4_dp*C_value*(gmod2*gmod2)/(total_rho(i_grid)**4)
1084 :
1085 : k = b_value*3.0_dp*pi*((total_rho(i_grid)/(9.0_dp*pi))**(1.0_dp/6.0_dp))
1086 : w0 = SQRT(wg2 + wp2/3.0_dp)
1087 :
1088 : q = w0/k
1089 :
1090 : !! Here, we calculate q0 by saturating q according
1091 : !! ---------------------------------------------------------------------------------------
1092 : q__q_cut = q/q_cut
1093 : CALL calculate_exponent(m_cut, q__q_cut, exponent)
1094 : q0(i_grid) = q_cut*(1.0_dp - EXP(-exponent))
1095 :
1096 : IF (q0(i_grid) < q_min) THEN
1097 : q0(i_grid) = q_min
1098 : END IF
1099 :
1100 : END IF
1101 :
1102 : END DO
1103 : !$OMP END PARALLEL DO
1104 :
1105 6 : END SUBROUTINE get_q0_on_grid_eo_rvv10
1106 :
1107 : ! **************************************************************************************************
1108 : !> \brief Prepare the spline intervals, weights and density factors for streamed theta construction.
1109 : !> \param q0 Saturated interpolation points
1110 : !> \param rho Total density
1111 : !> \param dispersion_env Non-local functional parameters
1112 : !> \param q_low Lower spline endpoint; zero denotes an inactive rVV10 point
1113 : !> \param coeff Cubic spline coefficients a, b, c, d
1114 : !> \param theta_scale Density factor multiplying the interpolated spline
1115 : ! **************************************************************************************************
1116 422 : SUBROUTINE prepare_splines(q0, rho, dispersion_env, q_low, coeff, theta_scale)
1117 : REAL(dp), INTENT(IN) :: q0(:), rho(:)
1118 : TYPE(qs_dispersion_type), POINTER :: dispersion_env
1119 : INTEGER, INTENT(OUT) :: q_low(:)
1120 : REAL(dp), INTENT(OUT) :: coeff(:, :), theta_scale(:)
1121 :
1122 : INTEGER :: i, j, lower, nqs, upper
1123 : LOGICAL :: rvv10
1124 : REAL(dp) :: a, b, const, dx, dx2_6
1125 : REAL(dp), POINTER :: q_mesh(:)
1126 :
1127 422 : q_mesh => dispersion_env%q_mesh
1128 422 : nqs = dispersion_env%nqs
1129 422 : CPASSERT(nqs >= 2)
1130 422 : rvv10 = dispersion_env%nl_type == vdw_nl_RVV10
1131 422 : const = 1.0_dp/(3.0_dp*rootpi*dispersion_env%b_value**1.5_dp)/(pi**0.75_dp)
1132 : !$OMP PARALLEL DO DEFAULT(NONE) &
1133 : !$OMP SHARED(q0, rho, q_mesh, nqs, rvv10, const, q_low, coeff, theta_scale) &
1134 422 : !$OMP PRIVATE(j, lower, upper, A, b, dx, dx2_6) SCHEDULE(STATIC)
1135 : DO i = 1, SIZE(q0)
1136 : IF (rvv10 .AND. rho(i) <= epsr) THEN
1137 : q_low(i) = 0
1138 : coeff(i, :) = 0.0_dp
1139 : theta_scale(i) = 0.0_dp
1140 : CYCLE
1141 : END IF
1142 : lower = 1
1143 : upper = nqs
1144 : DO WHILE (upper - lower > 1)
1145 : j = (upper + lower)/2
1146 : IF (q0(i) > q_mesh(j)) THEN
1147 : lower = j
1148 : ELSE
1149 : upper = j
1150 : END IF
1151 : END DO
1152 : q_low(i) = lower
1153 : dx = q_mesh(upper) - q_mesh(lower)
1154 : dx2_6 = dx*dx/6.0_dp
1155 : a = (q_mesh(upper) - q0(i))/dx
1156 : b = (q0(i) - q_mesh(lower))/dx
1157 : coeff(i, 1) = a
1158 : coeff(i, 2) = b
1159 : coeff(i, 3) = (a**3 - a)*dx2_6
1160 : coeff(i, 4) = (b**3 - b)*dx2_6
1161 : IF (rvv10) THEN
1162 : theta_scale(i) = const*rho(i)**0.75_dp
1163 : ELSE
1164 : theta_scale(i) = rho(i)
1165 : END IF
1166 : END DO
1167 : !$OMP END PARALLEL DO
1168 422 : END SUBROUTINE prepare_splines
1169 :
1170 : ! **************************************************************************************************
1171 : !> \brief Generate one real-space theta channel directly in the existing FFT workspace.
1172 : !> \param iq Channel index
1173 : !> \param q_low Lower spline endpoint
1174 : !> \param coeff Spline coefficients
1175 : !> \param theta_scale Density factor
1176 : !> \param dispersion_env Non-local functional parameters
1177 : !> \param theta FFT input, viewed as a contiguous one-dimensional array
1178 : ! **************************************************************************************************
1179 8440 : SUBROUTINE build_theta(iq, q_low, coeff, theta_scale, dispersion_env, theta)
1180 : INTEGER, INTENT(IN) :: iq, q_low(:)
1181 : REAL(dp), INTENT(IN) :: coeff(:, :), theta_scale(:)
1182 : TYPE(qs_dispersion_type), POINTER :: dispersion_env
1183 : REAL(dp), INTENT(OUT) :: theta(:)
1184 :
1185 : INTEGER :: i, lower
1186 : REAL(dp) :: p
1187 : REAL(dp), POINTER :: d2y(:, :)
1188 :
1189 8440 : d2y => dispersion_env%d2y_dx2
1190 : !$OMP PARALLEL DO SIMD DEFAULT(NONE) &
1191 8440 : !$OMP SHARED(iq, q_low, coeff, theta_scale, d2y, theta) PRIVATE(lower, p) SCHEDULE(STATIC)
1192 : DO i = 1, SIZE(theta)
1193 : lower = q_low(i)
1194 : theta(i) = 0.0_dp
1195 : IF (lower == 0) CYCLE
1196 : p = coeff(i, 1)*MERGE(1.0_dp, 0.0_dp, iq == lower) &
1197 : + coeff(i, 2)*MERGE(1.0_dp, 0.0_dp, iq == lower + 1) &
1198 : + (coeff(i, 3)*d2y(iq, lower) + coeff(i, 4)*d2y(iq, lower + 1))
1199 : theta(i) = p*theta_scale(i)
1200 : END DO
1201 : !$OMP END PARALLEL DO SIMD
1202 8440 : END SUBROUTINE build_theta
1203 :
1204 : ! **************************************************************************************************
1205 : !> \brief Accumulate one inverse-transformed channel without storing a real-space channel matrix.
1206 : !> \param iq Channel index
1207 : !> \param q_low Lower spline endpoint, using the potential-side knot convention
1208 : !> \param dispersion_env Non-local functional parameters
1209 : !> \param u Inverse FFT output
1210 : !> \param u_contract Sums against the two second-derivative columns, then the two endpoint values
1211 : ! **************************************************************************************************
1212 8320 : SUBROUTINE accumulate_potential(iq, q_low, dispersion_env, u, u_contract)
1213 : INTEGER, INTENT(IN) :: iq, q_low(:)
1214 : TYPE(qs_dispersion_type), POINTER :: dispersion_env
1215 : REAL(dp), INTENT(IN) :: u(:)
1216 : REAL(dp), INTENT(INOUT) :: u_contract(:, :)
1217 :
1218 : INTEGER :: i, lower
1219 : REAL(dp), POINTER :: d2y(:, :)
1220 :
1221 8320 : d2y => dispersion_env%d2y_dx2
1222 : !$OMP PARALLEL DO SIMD DEFAULT(NONE) &
1223 8320 : !$OMP SHARED(iq, q_low, d2y, u, u_contract) PRIVATE(lower) SCHEDULE(STATIC)
1224 : DO i = 1, SIZE(u)
1225 : lower = q_low(i)
1226 : IF (lower == 0) CYCLE
1227 : u_contract(i, 1) = u_contract(i, 1) + u(i)*d2y(iq, lower)
1228 : u_contract(i, 2) = u_contract(i, 2) + u(i)*d2y(iq, lower + 1)
1229 : IF (iq == lower) u_contract(i, 3) = u(i)
1230 : IF (iq == lower + 1) u_contract(i, 4) = u(i)
1231 : END DO
1232 : !$OMP END PARALLEL DO SIMD
1233 8320 : END SUBROUTINE accumulate_potential
1234 :
1235 : ! **************************************************************************************************
1236 : !> \brief This routine is modeled after an algorithm from "Numerical Recipes in C" by Cambridge
1237 : !> University Press, pages 96-97. It was adapted for Fortran and for the problem at hand.
1238 : !> \param x ...
1239 : !> \param d2y_dx2 ...
1240 : !> \par History
1241 : !> OpenMP added: Aug 2016 MTucker
1242 : ! **************************************************************************************************
1243 50 : SUBROUTINE initialize_spline_interpolation(x, d2y_dx2)
1244 :
1245 : REAL(dp), INTENT(in) :: x(:)
1246 : REAL(dp), INTENT(inout) :: d2y_dx2(:, :)
1247 :
1248 : INTEGER :: index, Nx, P_i
1249 : REAL(dp) :: temp1, temp2
1250 50 : REAL(dp), ALLOCATABLE :: temp_array(:), y(:)
1251 :
1252 50 : Nx = SIZE(x)
1253 :
1254 : !$OMP PARALLEL DEFAULT( NONE ) &
1255 : !$OMP SHARED( x, d2y_dx2, Nx ) &
1256 : !$OMP PRIVATE( temp_array, y &
1257 : !$OMP , index, temp1, temp2 &
1258 50 : !$OMP )
1259 :
1260 : ALLOCATE (temp_array(Nx), y(Nx))
1261 :
1262 : !$OMP DO
1263 : DO P_i = 1, Nx
1264 : !! In the Soler method, the polynomials that are interpolated are Kronecker delta functions
1265 : !! at a particular q point. So, we set all y values to 0 except the one corresponding to
1266 : !! the particular function P_i.
1267 : !! ----------------------------------------------------------------------------------------
1268 : y = 0.0_dp
1269 : y(P_i) = 1.0_dp
1270 : !! ----------------------------------------------------------------------------------------
1271 :
1272 : d2y_dx2(P_i, 1) = 0.0_dp
1273 : temp_array(1) = 0.0_dp
1274 : DO index = 2, Nx - 1
1275 : temp1 = (x(index) - x(index - 1))/(x(index + 1) - x(index - 1))
1276 : temp2 = temp1*d2y_dx2(P_i, index - 1) + 2.0_dp
1277 : d2y_dx2(P_i, index) = (temp1 - 1.0_dp)/temp2
1278 : temp_array(index) = (y(index + 1) - y(index))/(x(index + 1) - x(index)) &
1279 : - (y(index) - y(index - 1))/(x(index) - x(index - 1))
1280 : temp_array(index) = (6.0_dp*temp_array(index)/(x(index + 1) - x(index - 1)) &
1281 : - temp1*temp_array(index - 1))/temp2
1282 : END DO
1283 : d2y_dx2(P_i, Nx) = 0.0_dp
1284 : DO index = Nx - 1, 1, -1
1285 : d2y_dx2(P_i, index) = d2y_dx2(P_i, index)*d2y_dx2(P_i, index + 1) + temp_array(index)
1286 : END DO
1287 : END DO
1288 : !$OMP END DO
1289 :
1290 : DEALLOCATE (temp_array, y)
1291 : !$OMP END PARALLEL
1292 :
1293 50 : END SUBROUTINE initialize_spline_interpolation
1294 :
1295 : ! **************************************************************************************************
1296 : !> \brief Interpolate the symmetric kernel and optionally its radial derivative from packed tables.
1297 : !> \param k Reciprocal-vector length
1298 : !> \param kernel_of_k Interpolated kernel matrix
1299 : !> \param dispersion_env Radial mesh and channel count
1300 : !> \param kernel_table Packed channel pairs, radial points, and value/second derivative
1301 : !> \param dkernel_of_dk Optional derivative matrix for the kernel virial
1302 : ! **************************************************************************************************
1303 405916 : SUBROUTINE interpolate_kernel(k, kernel_of_k, dispersion_env, kernel_table, dkernel_of_dk)
1304 : REAL(dp), INTENT(IN) :: k
1305 : REAL(dp), INTENT(OUT) :: kernel_of_k(:, :)
1306 : TYPE(qs_dispersion_type), POINTER :: dispersion_env
1307 : REAL(dp), INTENT(IN) :: kernel_table(:, 0:, :)
1308 : REAL(dp), INTENT(OUT), OPTIONAL :: dkernel_of_dk(:, :)
1309 :
1310 : INTEGER :: ipair, k_i, q1_i, q2_i
1311 : LOGICAL :: on_mesh
1312 : REAL(dp) :: a, b, c, d, da, db, dc, dd, dk, dk_6, &
1313 : value
1314 :
1315 405916 : dk = dispersion_env%dk
1316 405916 : CPASSERT(k < dispersion_env%nr_points*dk)
1317 405916 : k_i = INT(k/dk)
1318 405916 : on_mesh = MOD(k, dk) == 0.0_dp
1319 405916 : a = (dk*(k_i + 1.0_dp) - k)/dk
1320 405916 : b = (k - dk*k_i)/dk
1321 405916 : c = (a**3 - a)*(dk*dk/6.0_dp)
1322 405916 : d = (b**3 - b)*(dk*dk/6.0_dp)
1323 8524236 : DO q1_i = 1, dispersion_env%nqs
1324 8524236 : !$OMP SIMD PRIVATE(ipair, value)
1325 : DO q2_i = 1, q1_i
1326 85242360 : ipair = q1_i*(q1_i - 1)/2 + q2_i
1327 85242360 : IF (on_mesh) THEN
1328 44310 : value = kernel_table(ipair, k_i, 1)
1329 : ELSE
1330 : value = a*kernel_table(ipair, k_i, 1) + b*kernel_table(ipair, k_i + 1, 1) &
1331 85198050 : + (c*kernel_table(ipair, k_i, 2) + d*kernel_table(ipair, k_i + 1, 2))
1332 : END IF
1333 85242360 : kernel_of_k(q2_i, q1_i) = value
1334 85242360 : kernel_of_k(q1_i, q2_i) = value
1335 : END DO
1336 : !$OMP END SIMD
1337 : END DO
1338 405916 : IF (PRESENT(dkernel_of_dk)) THEN
1339 299487 : dk_6 = dk/6.0_dp
1340 299487 : da = -1.0_dp/dk
1341 299487 : db = 1.0_dp/dk
1342 299487 : dc = -(3*a**2 - 1.0_dp)*dk_6
1343 299487 : dd = (3*b**2 - 1.0_dp)*dk_6
1344 6289227 : DO q1_i = 1, dispersion_env%nqs
1345 6289227 : !$OMP SIMD PRIVATE(ipair, value)
1346 : DO q2_i = 1, q1_i
1347 62892270 : ipair = q1_i*(q1_i - 1)/2 + q2_i
1348 : value = da*kernel_table(ipair, k_i, 1) + db*kernel_table(ipair, k_i + 1, 1) &
1349 62892270 : + dc*kernel_table(ipair, k_i, 2) + dd*kernel_table(ipair, k_i + 1, 2)
1350 62892270 : dkernel_of_dk(q2_i, q1_i) = value
1351 62892270 : dkernel_of_dk(q1_i, q2_i) = value
1352 : END DO
1353 : !$OMP END SIMD
1354 : END DO
1355 : END IF
1356 405916 : END SUBROUTINE interpolate_kernel
1357 :
1358 : ! **************************************************************************************************
1359 :
1360 : END MODULE qs_dispersion_nonloc
|
