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