(git:c641773)
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_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
73CONTAINS
74
75! **************************************************************************************************
76!> \brief ...
77!> \param dispersion_env ...
78!> \param para_env ...
79! **************************************************************************************************
80 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 CALL timeset(routinen, handle)
91
92 SELECT CASE (dispersion_env%nl_type)
93 CASE DEFAULT
94 cpabort("Unknown vdW-DF functional")
95 CASE (vdw_nl_drsll, vdw_nl_lmkll)
96 CALL cite_reference(dion2004)
97 CASE (vdw_nl_rvv10)
98 CALL cite_reference(sabatini2013)
99 END SELECT
100 CALL cite_reference(romanperez2009)
101
102 vdw_type = dispersion_env%type
103 SELECT CASE (vdw_type)
104 CASE DEFAULT
105 ! do nothing
106 CASE (xc_vdw_fun_nonloc)
107 ! setup information on non local functionals
108 filename = dispersion_env%kernel_file_name
109 IF (para_env%is_source()) THEN
110 ! Read the kernel information from file "filename"
111 CALL open_file(file_name=filename, unit_number=funit, file_form="FORMATTED")
112 READ (funit, *) nqs, nr_points
113 READ (funit, *) dispersion_env%r_max
114 END IF
115 CALL para_env%bcast(nqs)
116 CALL para_env%bcast(nr_points)
117 CALL para_env%bcast(dispersion_env%r_max)
118 ALLOCATE (dispersion_env%q_mesh(nqs), dispersion_env%kernel_table(nqs*(nqs + 1)/2, 0:nr_points, 2))
119 dispersion_env%nqs = nqs
120 dispersion_env%nr_points = nr_points
121 IF (para_env%is_source()) THEN
122 !! Read in the values of the q points used to generate this kernel
123 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 DO itable = 1, 2
127 DO ipair = 1, nqs*(nqs + 1)/2
128 READ (funit, "(1p, 4e23.14)") dispersion_env%kernel_table(ipair, 0:nr_points, itable)
129 END DO
130 END DO
131 CALL close_file(unit_number=funit)
132 END IF
133 CALL para_env%bcast(dispersion_env%q_mesh)
134 CALL para_env%bcast(dispersion_env%kernel_table)
135 ! 2nd derivates for interpolation
136 ALLOCATE (dispersion_env%d2y_dx2(nqs, nqs))
137 CALL initialize_spline_interpolation(dispersion_env%q_mesh, dispersion_env%d2y_dx2)
138 !
139 dispersion_env%q_cut = dispersion_env%q_mesh(nqs)
140 dispersion_env%q_min = dispersion_env%q_mesh(1)
141 dispersion_env%dk = 2.0_dp*pi/dispersion_env%r_max
142
143 END SELECT
144
145 CALL timestop(handle)
146
147 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 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 INTEGER, ALLOCATABLE, DIMENSION(:) :: q_low
179 LOGICAL :: use_virial
180 REAL(kind=dp) :: b_value, beta, ec_nl, sumnp
181 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: dq0_dgradrho, dq0_drho, hpot, q0, &
182 theta_scale
183 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: spline_coeff, u_contract
184 REAL(kind=dp), CONTIGUOUS, POINTER :: tmp_1d(:), vxc_1d(:)
185 TYPE(pw_c1d_gs_type) :: tmp_g
186 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 CALL timeset(routinen, handle)
191
192 IF (PRESENT(virial)) THEN
193 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 cpassert(.NOT. energy_only)
199 END IF
200
201 nl_type = dispersion_env%nl_type
202
203 b_value = dispersion_env%b_value
204 beta = 0.03125_dp*(3.0_dp/(b_value**2.0_dp))**0.75_dp
205
206 CALL pw_pool%create_pw(tmp_g)
207 CALL pw_pool%create_pw(tmp_r)
208 np = pw_pool%pw_grid%ngpts_local
209 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 IF (energy_only) THEN
220 ALLOCATE (q0(np))
221 SELECT CASE (nl_type)
222 CASE DEFAULT
223 cpabort("Unknown vdW-DF functional")
224 CASE (vdw_nl_drsll, vdw_nl_lmkll)
225 CALL get_q0_on_grid_eo_vdw(rho, drho, q0, dispersion_env)
226 CASE (vdw_nl_rvv10)
227 CALL get_q0_on_grid_eo_rvv10(rho, drho, q0, dispersion_env)
228 END SELECT
229 ELSE
230 ALLOCATE (q0(np), dq0_drho(np), dq0_dgradrho(np))
231 SELECT CASE (nl_type)
232 CASE DEFAULT
233 cpabort("Unknown vdW-DF functional")
234 CASE (vdw_nl_drsll, vdw_nl_lmkll)
235 CALL get_q0_on_grid_vdw(rho, drho, q0, dq0_drho, dq0_dgradrho, dispersion_env)
236 CASE (vdw_nl_rvv10)
237 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 ALLOCATE (q_low(np), spline_coeff(np, 4), theta_scale(np))
244 CALL prepare_splines(q0, rho, dispersion_env, q_low, spline_coeff, theta_scale)
245 ALLOCATE (thetas_g(dispersion_env%nqs))
246 CALL timeset("vdW_theta_forward", handle_fft)
247 DO i = 1, dispersion_env%nqs
248 CALL build_theta(i, q_low, spline_coeff, theta_scale, dispersion_env, tmp_1d)
249 CALL pw_pool%create_pw(thetas_g(i))
250 CALL pw_transfer(tmp_r, thetas_g(i))
251 END DO
252 CALL timestop(handle_fft)
253 DEALLOCATE (spline_coeff, theta_scale)
254 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 sumnp = np
261 CALL para_env%sum(sumnp)
262 IF (use_virial) THEN
263 ! calculates kernel contribution to stress
264 CALL vdw_energy(thetas_g, dispersion_env, ec_nl, energy_only, virial)
265 SELECT CASE (nl_type)
266 CASE (vdw_nl_rvv10)
267 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 DO idir = 1, 3
272 virial%pv_xc(idir, idir) = virial%pv_xc(idir, idir) + ec_nl
273 END DO
274 ELSE
275 CALL vdw_energy(thetas_g, dispersion_env, ec_nl, energy_only)
276 SELECT CASE (nl_type)
277 CASE (vdw_nl_rvv10)
278 ec_nl = 0.5_dp*ec_nl + beta*sum(rho(:))*grid%vol/sumnp
279 END SELECT
280 END IF
281 CALL para_env%sum(ec_nl)
282 IF (nl_type == vdw_nl_rvv10) ec_nl = ec_nl*dispersion_env%scale_rvv10
283 edispersion = ec_nl
284
285 IF (energy_only) THEN
286 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 ALLOCATE (u_contract(np, 4), hpot(np))
291!$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 CALL timeset("vdW_theta_inverse", handle_fft)
301 DO i = 1, dispersion_env%nqs
302 CALL pw_transfer(thetas_g(i), tmp_r)
303 CALL accumulate_potential(i, q_low, dispersion_env, tmp_1d, u_contract)
304 END DO
305 CALL timestop(handle_fft)
306
307 ! Write the local potential directly into its PW object.
308 vxc_1d(1:np) => vxc_r%array
309 IF (use_virial) THEN
310 grid => tmp_g%pw_grid
311 CALL get_potential(q0, dq0_drho, dq0_dgradrho, rho, q_low, u_contract, vxc_1d, hpot, &
312 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 dispersion_env)
316 END IF
317 DEALLOCATE (u_contract, q_low, q0, dq0_drho, dq0_dgradrho)
318 SELECT CASE (nl_type)
319 CASE (vdw_nl_rvv10)
320!$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 NULLIFY (vxc_1d)
328
329 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 DO i = 1, dispersion_env%nqs
336 CALL pw_pool%give_back_pw(thetas_g(i))
337 END DO
338 CALL pw_pool%give_back_pw(tmp_r)
339 CALL pw_pool%give_back_pw(tmp_g)
340
341 DEALLOCATE (thetas_g)
342
343 CALL timestop(handle)
344
345 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 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 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: fxc_vdw
371
372 CALL timeset(routinen, handle)
373
374 ! nonlocal vdW functionals
375 vdw_nl = .false.
376 CALL section_vals_val_get(xc_section, "VDW_POTENTIAL%POTENTIAL_TYPE", i_val=vdw)
377 vdw_nl = (vdw == xc_vdw_fun_nonloc)
378 IF (vdw_nl) THEN
379 nsteps = section_get_ival(xc_section, "NSTEPS")
380 order = 2*nsteps
381 eps_delta = section_get_rval(xc_section, "STEP_SIZE")
382 nspins = SIZE(rho0_r)
383 CALL get_qs_env(qs_env, pw_env=pw_env)
384 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
385 ALLOCATE (fxc_vdw(nspins))
386 DO ispin = 1, nspins
387 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 order, eps_delta, fxc_vdw)
391 DO ispin = 1, nspins
392 CALL pw_axpy(fxc_vdw(ispin), fxc_rho(ispin), 1.0_dp)
393 CALL auxbas_pw_pool%give_back_pw(fxc_vdw(ispin))
394 END DO
395 DEALLOCATE (fxc_vdw)
396
397 END IF
398
399 CALL timestop(handle)
400
401 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: rho, rho1, rhonl
443 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 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: vxc00
451
452 CALL timeset(routinen, handle)
453
454 ak = 0.0_dp
455 SELECT CASE (accuracy)
456 CASE (:4)
457 nstep = 2
458 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 nstep = 3
461 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 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 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 pw_env=pw_env, para_env=para_env)
470 CALL pw_env_get(pw_env, auxbas_pw_pool=pw_pool)
471 CALL pw_env_get(pw_env, vdw_pw_pool=vdw_pw_pool)
472
473 nspins = dft_control%nspins
474
475 ALLOCATE (vxc00(nspins))
476 DO ispin = 1, nspins
477 CALL pw_pool%create_pw(vxc00(ispin))
478 CALL pw_zero(fxc_rho(ispin))
479 END DO
480
481 np = vdw_pw_pool%pw_grid%ngpts_local
482 ALLOCATE (rho(np), drho(np, 3))
483 CALL prepare_nlvdw_density(rho0_r, rho, drho, vdw_pw_pool, pw_pool)
484 ALLOCATE (rho1(np), drho1(np, 3))
485 CALL prepare_nlvdw_density(rho1_r, rho1, drho1, vdw_pw_pool, pw_pool)
486 ALLOCATE (rhonl(np), drhonl(np, 3))
487
488 DO istep = -nstep, nstep
489
490 IF (ak(istep) /= 0.0_dp) THEN
491 beta = real(istep, kind=dp)*epsrho
492 rhonl(:) = rho(:) + beta*rho1(:)
493 drhonl(:, :) = drho(:, :) + beta*drho1(:, :)
494 CALL vdw_pw_pool%create_pw(vxcnl)
495 CALL calculate_vxc_nlvdw(vxcnl, rhonl, drhonl, edisp, dispersion_env, &
496 .false., vdw_pw_pool, para_env)
497 DO ispin = 1, nspins
498 CALL pw_zero(vxc00(ispin))
499 END DO
500 CALL return_nlvdw_potential(vxc00, vxcnl, vdw_pw_pool, pw_pool)
501 CALL vdw_pw_pool%give_back_pw(vxcnl)
502 DO ispin = 1, nspins
503 CALL pw_axpy(vxc00(ispin), fxc_rho(ispin), ak(istep))
504 END DO
505 END IF
506
507 END DO
508
509 DEALLOCATE (rho, drho)
510 DEALLOCATE (rho1, drho1)
511 DEALLOCATE (rhonl, drhonl)
512 DO ispin = 1, nspins
513 CALL pw_pool%give_back_pw(vxc00(ispin))
514 END DO
515 DEALLOCATE (vxc00)
516
517 oeps1 = 1.0_dp/epsrho
518 DO ispin = 1, nspins
519 CALL pw_scale(fxc_rho(ispin), oeps1)
520 END DO
521
522 CALL timestop(handle)
523
524 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 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 CALL timeset(routinen, handle)
549
550 nspins = SIZE(rho_r)
551
552 CALL pw_pool%create_pw(tmp_g)
553 CALL pw_pool%create_pw(tmp_r)
554
555 ! Sum the spin densities on the vdW grid
556 ! before transforming or differentiating.
557 CALL xc_pw_pool%create_pw(rho_g)
558 CALL pw_pool%create_pw(rho_tot_g)
559 CALL pw_zero(rho_tot_g)
560 DO ispin = 1, nspins
561 CALL pw_transfer(rho_r(ispin), rho_g)
562 CALL pw_transfer(rho_g, tmp_g)
563 CALL pw_axpy(tmp_g, rho_tot_g, 1._dp)
564 END DO
565 CALL pw_transfer(rho_tot_g, tmp_r)
566 CALL xc_pw_pool%give_back_pw(rho_g)
567
568 np = SIZE(tmp_r%array)
569 cpassert(np == SIZE(rho))
570 cpassert(np == SIZE(drho, 1))
571 DO i = 1, 3
572 lo(i) = lbound(tmp_r%array, i)
573 hi(i) = ubound(tmp_r%array, i)
574 n(i) = hi(i) - lo(i) + 1
575 END DO
576!$OMP PARALLEL DO DEFAULT(NONE) &
577!$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 DO idir = 1, 3
588 CALL pw_transfer(rho_tot_g, tmp_g)
589 CALL pw_derive(tmp_g, nd(:, idir))
590 CALL pw_transfer(tmp_g, tmp_r)
591!$OMP PARALLEL DO DEFAULT(NONE) &
592!$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 CALL pw_pool%give_back_pw(rho_tot_g)
604 CALL pw_pool%give_back_pw(tmp_r)
605 CALL pw_pool%give_back_pw(tmp_g)
606
607 CALL timestop(handle)
608
609 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 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 REAL(kind=dp), CONTIGUOUS, POINTER :: tmp_1d(:)
633 TYPE(pw_c1d_gs_type) :: div_g
634
635 CALL timeset(routinen, handle)
636
637 np = SIZE(tmp_r%array)
638 tmp_1d(1:np) => tmp_r%array
639 DO i = 1, 3
640 lo(i) = lbound(tmp_r%array, i)
641 hi(i) = ubound(tmp_r%array, i)
642 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 CALL pw_pool%create_pw(div_g)
648 CALL pw_zero(div_g)
649 DO idir = 1, 3
650!$OMP PARALLEL DO DEFAULT(NONE) &
651!$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 CALL pw_transfer(tmp_r, tmp_g)
662 CALL pw_derive(tmp_g, nd(:, idir))
663 CALL pw_axpy(tmp_g, div_g, 1._dp)
664 END DO
665
666 CALL pw_transfer(div_g, tmp_r)
667 CALL pw_axpy(tmp_r, vxc_r, -1._dp)
668
669 CALL pw_pool%give_back_pw(div_g)
670
671 CALL timestop(handle)
672
673 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 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 CALL timeset(routinen, handle)
694
695 nspins = SIZE(vxc_rho)
696
697 CALL pw_pool%create_pw(tmp_g)
698 CALL pw_transfer(vxc_r, tmp_g)
699 CALL xc_pw_pool%create_pw(vxc_g)
700 CALL pw_transfer(tmp_g, vxc_g)
701 CALL xc_pw_pool%create_pw(tmp_r)
702 CALL pw_transfer(vxc_g, tmp_r)
703 DO ispin = 1, nspins
704 CALL pw_axpy(tmp_r, vxc_rho(ispin), 1._dp)
705 END DO
706 CALL pw_pool%give_back_pw(tmp_g)
707 CALL xc_pw_pool%give_back_pw(vxc_g)
708 CALL xc_pw_pool%give_back_pw(tmp_r)
709
710 CALL timestop(handle)
711
712 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: theta_im, theta_re, u_im, u_re
742 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 CALL timeset(routinen, handle)
747 nqs = dispersion_env%nqs
748
749 use_virial = PRESENT(virial)
750 virial_thread(:, :) = 0.0_dp ! always initialize to avoid floating point exceptions in OMP REDUCTION
751
752 vdw_xc_energy = 0._dp
753 grid => thetas_g(1)%pw_grid
754
755 IF (grid%grid_span == halfspace) THEN
756 g_multiplier = 2._dp
757 ELSE
758 g_multiplier = 1._dp
759 END IF
760
761 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!$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 vdw_xc_energy = vdw_xc_energy*grid%vol*0.5_dp
842
843 IF (use_virial) THEN
844 DO l = 1, 3
845 DO m = 1, (l - 1)
846 virial%pv_xc(l, m) = virial%pv_xc(l, m) + virial_thread(l, m)
847 virial%pv_xc(m, l) = virial%pv_xc(l, m)
848 END DO
849 m = l
850 virial%pv_xc(l, m) = virial%pv_xc(l, m) + virial_thread(l, m)
851 END DO
852 END IF
853
854 CALL timestop(handle)
855
856 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 SUBROUTINE get_potential(q0, dq0_drho, dq0_dgradrho, total_rho, q_low_grid, u_contract, potential, h_prefactor, &
883 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 REAL(dp), DIMENSION(:), POINTER :: q_mesh
904
905 CALL timeset(routinen, handle)
906
907 use_virial = PRESENT(virial)
908 cpassert(.NOT. use_virial .OR. PRESENT(drho))
909 cpassert(.NOT. use_virial .OR. PRESENT(dvol))
910
911 virial_thread(:, :) = 0.0_dp ! always initialize to avoid floating point exceptions in OMP REDUCTION
912 b_value = dispersion_env%b_value
913 const = 1.0_dp/(3.0_dp*b_value**(3.0_dp/2.0_dp)*pi**(5.0_dp/4.0_dp))
914
915 q_mesh => dispersion_env%q_mesh
916 nqs = dispersion_env%nqs
917 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!$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 IF (use_virial) THEN
983 DO l = 1, 3
984 DO m = 1, (l - 1)
985 virial%pv_xc(l, m) = virial%pv_xc(l, m) + virial_thread(l, m)
986 virial%pv_xc(m, l) = virial%pv_xc(l, m)
987 END DO
988 m = l
989 virial%pv_xc(l, m) = virial%pv_xc(l, m) + virial_thread(l, m)
990 END DO
991 END IF
992
993 CALL timestop(handle)
994 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 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 multiplier = alpha
1013 exponent = alpha
1014
1015 DO i = 2, hi
1016 multiplier = multiplier*alpha
1017 exponent = exponent + (multiplier/i)
1018 END DO
1019 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 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 derivative = 0.0d0
1040 multiplier = 1.0d0
1041 exponent = 0.0d0
1042
1043 DO i = 1, hi
1044 derivative = derivative + multiplier
1045 multiplier = multiplier*alpha
1046 exponent = exponent + (multiplier/i)
1047 END DO
1048 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 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 q_cut = dispersion_env%q_cut
1084 q_min = dispersion_env%q_min
1085 SELECT CASE (dispersion_env%nl_type)
1086 CASE DEFAULT
1087 cpabort("Unknown vdW-DF functional")
1088 CASE (vdw_nl_drsll)
1089 z_ab = -0.8491_dp
1090 CASE (vdw_nl_lmkll)
1091 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!$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 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 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 q_cut = dispersion_env%q_cut
1195 q_min = dispersion_env%q_min
1196 b_value = dispersion_env%b_value
1197 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!$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 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 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 q_cut = dispersion_env%q_cut
1276 q_min = dispersion_env%q_min
1277 SELECT CASE (dispersion_env%nl_type)
1278 CASE DEFAULT
1279 cpabort("Unknown vdW-DF functional")
1280 CASE (vdw_nl_drsll)
1281 z_ab = -0.8491_dp
1282 CASE (vdw_nl_lmkll)
1283 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!$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 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 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 q_cut = dispersion_env%q_cut
1356 q_min = dispersion_env%q_min
1357 b_value = dispersion_env%b_value
1358 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!$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 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 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 q_mesh => dispersion_env%q_mesh
1420 nqs = dispersion_env%nqs
1421 cpassert(nqs >= 2)
1422 rvv10 = dispersion_env%nl_type == vdw_nl_rvv10
1423 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!$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 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 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 d2y => dispersion_env%d2y_dx2
1482!$OMP PARALLEL DO SIMD DEFAULT(NONE) &
1483!$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 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 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 d2y => dispersion_env%d2y_dx2
1514!$OMP PARALLEL DO SIMD DEFAULT(NONE) &
1515!$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 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 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 REAL(dp), ALLOCATABLE :: temp_array(:), y(:)
1543
1544 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!$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 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 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 dk = dispersion_env%dk
1608 cpassert(k < dispersion_env%nr_points*dk)
1609 k_i = int(k/dk)
1610 on_mesh = mod(k, dk) == 0.0_dp
1611 a = (dk*(k_i + 1.0_dp) - k)/dk
1612 b = (k - dk*k_i)/dk
1613 c = (a**3 - a)*(dk*dk/6.0_dp)
1614 d = (b**3 - b)*(dk*dk/6.0_dp)
1615 DO q1_i = 1, dispersion_env%nqs
1616!$OMP SIMD PRIVATE(ipair, value)
1617 DO q2_i = 1, q1_i
1618 ipair = q1_i*(q1_i - 1)/2 + q2_i
1619 IF (on_mesh) THEN
1620 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 + (c*kernel_table(ipair, k_i, 2) + d*kernel_table(ipair, k_i + 1, 2))
1624 END IF
1625 kernel_of_k(q2_i, q1_i) = value
1626 kernel_of_k(q1_i, q2_i) = value
1627 END DO
1628!$OMP END SIMD
1629 END DO
1630 IF (PRESENT(dkernel_of_dk)) THEN
1631 dk_6 = dk/6.0_dp
1632 da = -1.0_dp/dk
1633 db = 1.0_dp/dk
1634 dc = -(3*a**2 - 1.0_dp)*dk_6
1635 dd = (3*b**2 - 1.0_dp)*dk_6
1636 DO q1_i = 1, dispersion_env%nqs
1637!$OMP SIMD PRIVATE(ipair, value)
1638 DO q2_i = 1, q1_i
1639 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 + dc*kernel_table(ipair, k_i, 2) + dd*kernel_table(ipair, k_i + 1, 2)
1642 dkernel_of_dk(q2_i, q1_i) = value
1643 dkernel_of_dk(q1_i, q2_i) = value
1644 END DO
1645!$OMP END SIMD
1646 END DO
1647 END IF
1648 END SUBROUTINE interpolate_kernel
1649
1650! **************************************************************************************************
1651
1652END MODULE qs_dispersion_nonloc
pw_methods::pw_axpy
Definition pw_methods.F:165
pw_methods::pw_scale
Definition pw_methods.F:94
pw_methods::pw_transfer
Definition pw_methods.F:304
pw_methods::pw_zero
Definition pw_methods.F:83
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_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_files
Utility routines to open and close files.
Definition cp_files.F:22
cp_files::open_file
subroutine, public open_file(file_name, unit_number, file_status, file_form, file_action, file_position, file_pad, file_access, debug)
Opens a file and returns its logical unit number.
Definition cp_files.F:197
cp_files::close_file
subroutine, public close_file(unit_number, file_status)
Closes an open file given by its logical unit number.
Definition cp_files.F:64
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:657
input_constants::xc_vdw_fun_nonloc
integer, parameter, public xc_vdw_fun_nonloc
Definition input_constants.F:649
input_constants::vdw_nl_drsll
integer, parameter, public vdw_nl_drsll
Definition input_constants.F:657
input_constants::vdw_nl_lmkll
integer, parameter, public vdw_nl_lmkll
Definition input_constants.F:657
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_get_rval
real(kind=dp) function, public section_get_rval(section_vals, keyword_name)
...
Definition input_section_types.F:953
input_section_types::section_get_ival
integer function, public section_get_ival(section_vals, keyword_name)
...
Definition input_section_types.F:985
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1052
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_env_types
container for various plainwaves related things
Definition pw_env_types.F:14
pw_env_types::pw_env_get
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Definition pw_env_types.F:113
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:10149
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:81
qs_dispersion_nonloc::prepare_nlvdw_density
subroutine, public prepare_nlvdw_density(rho_r, rho, drho, pw_pool, xc_pw_pool)
Transform the input density to a 1d grid used by the nl-vdW routines.
Definition qs_dispersion_nonloc.F:535
qs_dispersion_nonloc::qs_fxc_nlvdw_apply
subroutine, public qs_fxc_nlvdw_apply(qs_env, xc_section, dispersion_env, rho0_r, rho1_r, fxc_rho)
Second derivative of nl-vdW potential.
Definition qs_dispersion_nonloc.F:357
qs_dispersion_nonloc::qs_fxc_nlvdw_fdiff
subroutine, public qs_fxc_nlvdw_fdiff(qs_env, dispersion_env, rho0_r, rho1_r, accuracy, epsrho, fxc_rho)
Second derivative of nl-vdW potential using finite differences.
Definition qs_dispersion_nonloc.F:430
qs_dispersion_nonloc::calculate_vxc_nlvdw
subroutine, public calculate_vxc_nlvdw(vxc_r, rho, drho, edispersion, dispersion_env, energy_only, 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:164
qs_dispersion_nonloc::return_nlvdw_potential
subroutine, public return_nlvdw_potential(vxc_rho, vxc_r, pw_pool, xc_pw_pool)
Take the nl-vdW potential and return it on the XC grid and for all spin channels.
Definition qs_dispersion_nonloc.F:683
qs_dispersion_types
Definition of disperson types for DFT calculations.
Definition qs_dispersion_types.F:12
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:544
virial_types
Definition virial_types.F:13
cp_control_types::dft_control_type
Definition cp_control_types.F:786
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:774
pw_env_types::pw_env_type
contained for different pw related things
Definition pw_env_types.F:70
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
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:223
virial_types::virial_type
Definition virial_types.F:26