(git:744416f)
Loading...
Searching...
No Matches
qs_dispersion_nonloc.F
Go to the documentation of this file.
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! **************************************************************************************************
19MODULE 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
59CONTAINS
60
61! **************************************************************************************************
62!> \brief ...
63!> \param dispersion_env ...
64!> \param para_env ...
65! **************************************************************************************************
66 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 CALL timeset(routinen, handle)
77
78 SELECT CASE (dispersion_env%nl_type)
79 CASE DEFAULT
80 cpabort("Unknown vdW-DF functional")
81 CASE (vdw_nl_drsll, vdw_nl_lmkll)
82 CALL cite_reference(dion2004)
83 CASE (vdw_nl_rvv10)
84 CALL cite_reference(sabatini2013)
85 END SELECT
86 CALL cite_reference(romanperez2009)
87
88 vdw_type = dispersion_env%type
89 SELECT CASE (vdw_type)
90 CASE DEFAULT
91 ! do nothing
92 CASE (xc_vdw_fun_nonloc)
93 ! setup information on non local functionals
94 filename = dispersion_env%kernel_file_name
95 IF (para_env%is_source()) THEN
96 ! Read the kernel information from file "filename"
97 CALL open_file(file_name=filename, unit_number=funit, file_form="FORMATTED")
98 READ (funit, *) nqs, nr_points
99 READ (funit, *) dispersion_env%r_max
100 END IF
101 CALL para_env%bcast(nqs)
102 CALL para_env%bcast(nr_points)
103 CALL para_env%bcast(dispersion_env%r_max)
104 ALLOCATE (dispersion_env%q_mesh(nqs), dispersion_env%kernel_table(nqs*(nqs + 1)/2, 0:nr_points, 2))
105 dispersion_env%nqs = nqs
106 dispersion_env%nr_points = nr_points
107 IF (para_env%is_source()) THEN
108 !! Read in the values of the q points used to generate this kernel
109 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 DO itable = 1, 2
113 DO ipair = 1, nqs*(nqs + 1)/2
114 READ (funit, "(1p, 4e23.14)") dispersion_env%kernel_table(ipair, 0:nr_points, itable)
115 END DO
116 END DO
117 CALL close_file(unit_number=funit)
118 END IF
119 CALL para_env%bcast(dispersion_env%q_mesh)
120 CALL para_env%bcast(dispersion_env%kernel_table)
121 ! 2nd derivates for interpolation
122 ALLOCATE (dispersion_env%d2y_dx2(nqs, nqs))
123 CALL initialize_spline_interpolation(dispersion_env%q_mesh, dispersion_env%d2y_dx2)
124 !
125 dispersion_env%q_cut = dispersion_env%q_mesh(nqs)
126 dispersion_env%q_min = dispersion_env%q_mesh(1)
127 dispersion_env%dk = 2.0_dp*pi/dispersion_env%r_max
128
129 END SELECT
130
131 CALL timestop(handle)
132
133 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 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: dq0_dgradrho, dq0_drho, hpot, q0, rho, &
170 theta_scale
171 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: drho, spline_coeff, u_contract
172 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 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 CALL timeset(routinen, handle)
179
180 cpassert(ASSOCIATED(rho_r))
181 cpassert(ASSOCIATED(rho_g))
182 cpassert(ASSOCIATED(pw_pool))
183
184 IF (PRESENT(virial)) THEN
185 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 cpassert(.NOT. energy_only)
191 END IF
192 IF (.NOT. energy_only) THEN
193 cpassert(ASSOCIATED(vxc_rho))
194 END IF
195
196 nl_type = dispersion_env%nl_type
197
198 b_value = dispersion_env%b_value
199 beta = 0.03125_dp*(3.0_dp/(b_value**2.0_dp))**0.75_dp
200 nspin = SIZE(rho_r)
201
202 ! temporary arrays for FFT
203 CALL pw_pool%create_pw(tmp_g)
204 CALL pw_pool%create_pw(tmp_r)
205
206 ! Sum the spin densities on the vdW grid before transforming or differentiating.
207 CALL pw_pool%create_pw(rho_tot_g)
208 CALL pw_transfer(rho_g(1), rho_tot_g)
209 DO ispin = 2, nspin
210 CALL pw_transfer(rho_g(ispin), tmp_g)
211 CALL pw_axpy(tmp_g, rho_tot_g, 1._dp)
212 END DO
213 CALL pw_transfer(rho_tot_g, tmp_r)
214
215 np = SIZE(tmp_r%array)
216 tmp_1d(1:np) => tmp_r%array
217 ALLOCATE (rho(np), drho(np, 3))
218 DO i = 1, 3
219 lo(i) = lbound(tmp_r%array, i)
220 hi(i) = ubound(tmp_r%array, i)
221 n(i) = hi(i) - lo(i) + 1
222 END DO
223!$OMP PARALLEL DO DEFAULT(NONE) &
224!$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 DO idir = 1, 3
235 CALL pw_transfer(rho_tot_g, tmp_g)
236 CALL pw_derive(tmp_g, nd(:, idir))
237 CALL pw_transfer(tmp_g, tmp_r)
238!$OMP PARALLEL DO DEFAULT(NONE) &
239!$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 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 IF (energy_only) THEN
261 ALLOCATE (q0(np))
262 SELECT CASE (nl_type)
263 CASE DEFAULT
264 cpabort("Unknown vdW-DF functional")
265 CASE (vdw_nl_drsll, vdw_nl_lmkll)
266 CALL get_q0_on_grid_eo_vdw(rho, drho, q0, dispersion_env)
267 CASE (vdw_nl_rvv10)
268 CALL get_q0_on_grid_eo_rvv10(rho, drho, q0, dispersion_env)
269 END SELECT
270 ELSE
271 ALLOCATE (q0(np), dq0_drho(np), dq0_dgradrho(np))
272 SELECT CASE (nl_type)
273 CASE DEFAULT
274 cpabort("Unknown vdW-DF functional")
275 CASE (vdw_nl_drsll, vdw_nl_lmkll)
276 CALL get_q0_on_grid_vdw(rho, drho, q0, dq0_drho, dq0_dgradrho, dispersion_env)
277 CASE (vdw_nl_rvv10)
278 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 ALLOCATE (q_low(np), spline_coeff(np, 4), theta_scale(np))
285 CALL prepare_splines(q0, rho, dispersion_env, q_low, spline_coeff, theta_scale)
286 ALLOCATE (thetas_g(dispersion_env%nqs))
287 CALL timeset("vdW_theta_forward", handle_fft)
288 DO i = 1, dispersion_env%nqs
289 CALL build_theta(i, q_low, spline_coeff, theta_scale, dispersion_env, tmp_1d)
290 CALL pw_pool%create_pw(thetas_g(i))
291 CALL pw_transfer(tmp_r, thetas_g(i))
292 END DO
293 CALL timestop(handle_fft)
294 DEALLOCATE (spline_coeff, theta_scale)
295 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 sumnp = np
302 CALL para_env%sum(sumnp)
303 IF (use_virial) THEN
304 ! calculates kernel contribution to stress
305 CALL vdw_energy(thetas_g, dispersion_env, ec_nl, energy_only, virial)
306 SELECT CASE (nl_type)
307 CASE (vdw_nl_rvv10)
308 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 DO idir = 1, 3
313 virial%pv_xc(idir, idir) = virial%pv_xc(idir, idir) + ec_nl
314 END DO
315 ELSE
316 CALL vdw_energy(thetas_g, dispersion_env, ec_nl, energy_only)
317 SELECT CASE (nl_type)
318 CASE (vdw_nl_rvv10)
319 ec_nl = 0.5_dp*ec_nl + beta*sum(rho(:))*grid%vol/sumnp
320 END SELECT
321 END IF
322 CALL para_env%sum(ec_nl)
323 IF (nl_type == vdw_nl_rvv10) ec_nl = ec_nl*dispersion_env%scale_rvv10
324 edispersion = ec_nl
325
326 IF (energy_only) THEN
327 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 ALLOCATE (u_contract(np, 4), hpot(np))
332!$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 CALL timeset("vdW_theta_inverse", handle_fft)
342 DO i = 1, dispersion_env%nqs
343 CALL pw_transfer(thetas_g(i), tmp_r)
344 CALL accumulate_potential(i, q_low, dispersion_env, tmp_1d, u_contract)
345 END DO
346 CALL timestop(handle_fft)
347
348 ! Write the local potential directly into its PW object.
349 CALL pw_pool%create_pw(vxc_r)
350 vxc_1d(1:np) => vxc_r%array
351 IF (use_virial) THEN
352 grid => tmp_g%pw_grid
353 CALL get_potential(q0, dq0_drho, dq0_dgradrho, rho, q_low, u_contract, vxc_1d, hpot, &
354 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 dispersion_env)
358 END IF
359 DEALLOCATE (u_contract, q_low, q0, dq0_drho, dq0_dgradrho)
360 SELECT CASE (nl_type)
361 CASE (vdw_nl_rvv10)
362!$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 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 CALL pw_pool%create_pw(div_g)
373 DO idir = 1, 3
374!$OMP PARALLEL DO DEFAULT(NONE) &
375!$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 CALL pw_transfer(tmp_r, tmp_g)
386 CALL pw_derive(tmp_g, nd(:, idir))
387 IF (idir == 1) THEN
388 CALL pw_transfer(tmp_g, div_g)
389 ELSE
390 CALL pw_axpy(tmp_g, div_g, 1._dp)
391 END IF
392 END DO
393 CALL pw_transfer(div_g, tmp_r)
394 CALL pw_pool%give_back_pw(div_g)
395 CALL pw_axpy(tmp_r, vxc_r, -1._dp)
396 CALL pw_transfer(vxc_r, tmp_g)
397 CALL pw_pool%give_back_pw(vxc_r)
398 CALL xc_pw_pool%create_pw(vxc_r)
399 CALL xc_pw_pool%create_pw(vxc_g)
400 CALL pw_transfer(tmp_g, vxc_g)
401 CALL pw_transfer(vxc_g, vxc_r)
402 DO ispin = 1, nspin
403 CALL pw_axpy(vxc_r, vxc_rho(ispin), 1._dp)
404 END DO
405 CALL xc_pw_pool%give_back_pw(vxc_r)
406 CALL xc_pw_pool%give_back_pw(vxc_g)
407 END IF
408
409 NULLIFY (tmp_1d)
410
411 DO i = 1, dispersion_env%nqs
412 CALL pw_pool%give_back_pw(thetas_g(i))
413 END DO
414 CALL pw_pool%give_back_pw(tmp_r)
415 CALL pw_pool%give_back_pw(tmp_g)
416
417 DEALLOCATE (rho, drho, thetas_g)
418
419 CALL timestop(handle)
420
421 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: theta_im, theta_re, u_im, u_re
450 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 CALL timeset(routinen, handle)
455 nqs = dispersion_env%nqs
456
457 use_virial = PRESENT(virial)
458 virial_thread(:, :) = 0.0_dp ! always initialize to avoid floating point exceptions in OMP REDUCTION
459
460 vdw_xc_energy = 0._dp
461 grid => thetas_g(1)%pw_grid
462
463 IF (grid%grid_span == halfspace) THEN
464 g_multiplier = 2._dp
465 ELSE
466 g_multiplier = 1._dp
467 END IF
468
469 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!$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 vdw_xc_energy = vdw_xc_energy*grid%vol*0.5_dp
550
551 IF (use_virial) THEN
552 DO l = 1, 3
553 DO m = 1, (l - 1)
554 virial%pv_xc(l, m) = virial%pv_xc(l, m) + virial_thread(l, m)
555 virial%pv_xc(m, l) = virial%pv_xc(l, m)
556 END DO
557 m = l
558 virial%pv_xc(l, m) = virial%pv_xc(l, m) + virial_thread(l, m)
559 END DO
560 END IF
561
562 CALL timestop(handle)
563
564 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 SUBROUTINE get_potential(q0, dq0_drho, dq0_dgradrho, total_rho, q_low_grid, u_contract, potential, h_prefactor, &
591 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 REAL(dp), DIMENSION(:), POINTER :: q_mesh
612
613 CALL timeset(routinen, handle)
614
615 use_virial = PRESENT(virial)
616 cpassert(.NOT. use_virial .OR. PRESENT(drho))
617 cpassert(.NOT. use_virial .OR. PRESENT(dvol))
618
619 virial_thread(:, :) = 0.0_dp ! always initialize to avoid floating point exceptions in OMP REDUCTION
620 b_value = dispersion_env%b_value
621 const = 1.0_dp/(3.0_dp*b_value**(3.0_dp/2.0_dp)*pi**(5.0_dp/4.0_dp))
622
623 q_mesh => dispersion_env%q_mesh
624 nqs = dispersion_env%nqs
625 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!$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 IF (use_virial) THEN
691 DO l = 1, 3
692 DO m = 1, (l - 1)
693 virial%pv_xc(l, m) = virial%pv_xc(l, m) + virial_thread(l, m)
694 virial%pv_xc(m, l) = virial%pv_xc(l, m)
695 END DO
696 m = l
697 virial%pv_xc(l, m) = virial%pv_xc(l, m) + virial_thread(l, m)
698 END DO
699 END IF
700
701 CALL timestop(handle)
702 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 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 multiplier = alpha
721 exponent = alpha
722
723 DO i = 2, hi
724 multiplier = multiplier*alpha
725 exponent = exponent + (multiplier/i)
726 END DO
727 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 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 derivative = 0.0d0
748 multiplier = 1.0d0
749 exponent = 0.0d0
750
751 DO i = 1, hi
752 derivative = derivative + multiplier
753 multiplier = multiplier*alpha
754 exponent = exponent + (multiplier/i)
755 END DO
756 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 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 q_cut = dispersion_env%q_cut
792 q_min = dispersion_env%q_min
793 SELECT CASE (dispersion_env%nl_type)
794 CASE DEFAULT
795 cpabort("Unknown vdW-DF functional")
796 CASE (vdw_nl_drsll)
797 z_ab = -0.8491_dp
798 CASE (vdw_nl_lmkll)
799 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!$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 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 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 q_cut = dispersion_env%q_cut
903 q_min = dispersion_env%q_min
904 b_value = dispersion_env%b_value
905 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!$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 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 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 q_cut = dispersion_env%q_cut
984 q_min = dispersion_env%q_min
985 SELECT CASE (dispersion_env%nl_type)
986 CASE DEFAULT
987 cpabort("Unknown vdW-DF functional")
988 CASE (vdw_nl_drsll)
989 z_ab = -0.8491_dp
990 CASE (vdw_nl_lmkll)
991 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!$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 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 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 q_cut = dispersion_env%q_cut
1064 q_min = dispersion_env%q_min
1065 b_value = dispersion_env%b_value
1066 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!$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 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 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 q_mesh => dispersion_env%q_mesh
1128 nqs = dispersion_env%nqs
1129 cpassert(nqs >= 2)
1130 rvv10 = dispersion_env%nl_type == vdw_nl_rvv10
1131 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!$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 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 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 d2y => dispersion_env%d2y_dx2
1190!$OMP PARALLEL DO SIMD DEFAULT(NONE) &
1191!$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 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 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 d2y => dispersion_env%d2y_dx2
1222!$OMP PARALLEL DO SIMD DEFAULT(NONE) &
1223!$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 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 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 REAL(dp), ALLOCATABLE :: temp_array(:), y(:)
1251
1252 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!$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 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 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 dk = dispersion_env%dk
1316 cpassert(k < dispersion_env%nr_points*dk)
1317 k_i = int(k/dk)
1318 on_mesh = mod(k, dk) == 0.0_dp
1319 a = (dk*(k_i + 1.0_dp) - k)/dk
1320 b = (k - dk*k_i)/dk
1321 c = (a**3 - a)*(dk*dk/6.0_dp)
1322 d = (b**3 - b)*(dk*dk/6.0_dp)
1323 DO q1_i = 1, dispersion_env%nqs
1324!$OMP SIMD PRIVATE(ipair, value)
1325 DO q2_i = 1, q1_i
1326 ipair = q1_i*(q1_i - 1)/2 + q2_i
1327 IF (on_mesh) THEN
1328 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 + (c*kernel_table(ipair, k_i, 2) + d*kernel_table(ipair, k_i + 1, 2))
1332 END IF
1333 kernel_of_k(q2_i, q1_i) = value
1334 kernel_of_k(q1_i, q2_i) = value
1335 END DO
1336!$OMP END SIMD
1337 END DO
1338 IF (PRESENT(dkernel_of_dk)) THEN
1339 dk_6 = dk/6.0_dp
1340 da = -1.0_dp/dk
1341 db = 1.0_dp/dk
1342 dc = -(3*a**2 - 1.0_dp)*dk_6
1343 dd = (3*b**2 - 1.0_dp)*dk_6
1344 DO q1_i = 1, dispersion_env%nqs
1345!$OMP SIMD PRIVATE(ipair, value)
1346 DO q2_i = 1, q1_i
1347 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 + dc*kernel_table(ipair, k_i, 2) + dd*kernel_table(ipair, k_i + 1, 2)
1350 dkernel_of_dk(q2_i, q1_i) = value
1351 dkernel_of_dk(q1_i, q2_i) = value
1352 END DO
1353!$OMP END SIMD
1354 END DO
1355 END IF
1356 END SUBROUTINE interpolate_kernel
1357
1358! **************************************************************************************************
1359
1360END MODULE qs_dispersion_nonloc
pw_methods::pw_axpy
Definition pw_methods.F:165
pw_methods::pw_transfer
Definition pw_methods.F:304
bibliography
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition bibliography.F:21
bibliography::dion2004
integer, save, public dion2004
Definition bibliography.F:36
bibliography::romanperez2009
integer, save, public romanperez2009
Definition bibliography.F:36
bibliography::sabatini2013
integer, save, public sabatini2013
Definition bibliography.F:36
cp_files
Utility routines to open and close files. Tracking of preconnections.
Definition cp_files.F:16
cp_files::open_file
subroutine, public open_file(file_name, file_status, file_form, file_action, file_position, file_pad, unit_number, debug, skip_get_unit_number, file_access)
Opens the requested file using a free unit number.
Definition cp_files.F:311
cp_files::close_file
subroutine, public close_file(unit_number, file_status, keep_preconnection)
Close an open file given by its logical unit number. Optionally, keep the file and unit preconnected.
Definition cp_files.F:122
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::vdw_nl_rvv10
integer, parameter, public vdw_nl_rvv10
Definition input_constants.F:644
input_constants::xc_vdw_fun_nonloc
integer, parameter, public xc_vdw_fun_nonloc
Definition input_constants.F:637
input_constants::vdw_nl_drsll
integer, parameter, public vdw_nl_drsll
Definition input_constants.F:644
input_constants::vdw_nl_lmkll
integer, parameter, public vdw_nl_lmkll
Definition input_constants.F:644
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_path_length
integer, parameter, public default_path_length
Definition kinds.F:58
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
mathconstants::rootpi
real(kind=dp), parameter, public rootpi
Definition mathconstants.F:114
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:27
pw_grid_types
Definition pw_grid_types.F:13
pw_grid_types::halfspace
integer, parameter, public halfspace
Definition pw_grid_types.F:26
pw_methods
Definition pw_methods.F:25
pw_methods::pw_derive
subroutine, public pw_derive(pw, n)
Calculate the derivative of a plane wave vector.
Definition pw_methods.F:10257
pw_pool_types
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:24
pw_types
Definition pw_types.F:24
qs_dispersion_nonloc
Calculation of non local dispersion functionals Some routines adapted from: Copyright (C) 2001-2009 Q...
Definition qs_dispersion_nonloc.F:19
qs_dispersion_nonloc::qs_dispersion_nonloc_init
subroutine, public qs_dispersion_nonloc_init(dispersion_env, para_env)
...
Definition qs_dispersion_nonloc.F:67
qs_dispersion_nonloc::calculate_dispersion_nonloc
subroutine, public calculate_dispersion_nonloc(vxc_rho, rho_r, rho_g, edispersion, dispersion_env, energy_only, pw_pool, xc_pw_pool, para_env, virial)
Calculates the non-local vdW functional using the method of Soler For spin polarized cases we use E(a...
Definition qs_dispersion_nonloc.F:151
qs_dispersion_types
Definition of disperson types for DFT calculations.
Definition qs_dispersion_types.F:12
virial_types
Definition virial_types.F:13
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:768
pw_grid_types::pw_grid_type
Definition pw_grid_types.F:53
pw_pool_types::pw_pool_type
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:83
pw_types::pw_c1d_gs_type
Definition pw_types.F:127
pw_types::pw_r3d_rs_type
Definition pw_types.F:62
qs_dispersion_types::qs_dispersion_type
Definition qs_dispersion_types.F:34
virial_types::virial_type
Definition virial_types.F:26