(git:b279b6b)
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-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 
8 ! **************************************************************************************************
9 !> \brief Calculation of non local dispersion functionals
10 !> Some routines adapted from:
11 !> Copyright (C) 2001-2009 Quantum ESPRESSO group
12 !> Copyright (C) 2009 Brian Kolb, Timo Thonhauser - Wake Forest University
13 !> This file is distributed under the terms of the
14 !> GNU General Public License. See the file `License'
15 !> in the root directory of the present distribution,
16 !> or http://www.gnu.org/copyleft/gpl.txt .
17 !> \author JGH
18 ! **************************************************************************************************
19 MODULE qs_dispersion_nonloc
20  USE bibliography, ONLY: dion2004,&
21  romanperez2009,&
22  sabatini2013,&
23  cite_reference
24  USE cp_files, ONLY: close_file,&
25  open_file
26  USE input_constants, ONLY: vdw_nl_drsll,&
27  vdw_nl_lmkll,&
28  vdw_nl_rvv10,&
29  xc_vdw_fun_nonloc
30  USE kinds, ONLY: default_string_length,&
31  dp
32  USE mathconstants, ONLY: pi,&
33  rootpi
34  USE message_passing, ONLY: mp_para_env_type
35  USE pw_grid_types, ONLY: halfspace,&
36  pw_grid_type
37  USE pw_methods, ONLY: pw_axpy,&
38  pw_derive,&
39  pw_transfer
40  USE pw_pool_types, ONLY: pw_pool_type
41  USE pw_types, ONLY: pw_c1d_gs_type,&
42  pw_r3d_rs_type
43  USE qs_dispersion_types, ONLY: qs_dispersion_type
44  USE virial_types, ONLY: virial_type
45 #include "./base/base_uses.f90"
46 
47  IMPLICIT NONE
48 
49  PRIVATE
50 
51  REAL(KIND=dp), PARAMETER :: epsr = 1.e-12_dp
52 
53  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_dispersion_nonloc'
54 
55  PUBLIC :: qs_dispersion_nonloc_init, calculate_dispersion_nonloc
56 
57 ! **************************************************************************************************
58 
59 CONTAINS
60 
61 ! **************************************************************************************************
62 !> \brief ...
63 !> \param dispersion_env ...
64 !> \param para_env ...
65 ! **************************************************************************************************
66  SUBROUTINE qs_dispersion_nonloc_init(dispersion_env, para_env)
67  TYPE(qs_dispersion_type), POINTER :: dispersion_env
68  TYPE(mp_para_env_type), POINTER :: para_env
69 
70  CHARACTER(len=*), PARAMETER :: routinen = 'qs_dispersion_nonloc_init'
71 
72  CHARACTER(LEN=default_string_length) :: filename
73  INTEGER :: funit, handle, nqs, nr_points, q1_i, &
74  q2_i, vdw_type
75 
76  CALL timeset(routinen, handle)
77 
78  SELECT CASE (dispersion_env%nl_type)
79  CASE DEFAULT
80  cpabort("Unknown vdW-DF functional")
81  CASE (vdw_nl_drsll, vdw_nl_lmkll)
82  CALL cite_reference(dion2004)
83  CASE (vdw_nl_rvv10)
84  CALL cite_reference(sabatini2013)
85  END SELECT
86  CALL cite_reference(romanperez2009)
87 
88  vdw_type = dispersion_env%type
89  SELECT CASE (vdw_type)
90  CASE DEFAULT
91  ! do nothing
92  CASE (xc_vdw_fun_nonloc)
93  ! setup information on non local functionals
94  filename = dispersion_env%kernel_file_name
95  IF (para_env%is_source()) THEN
96  ! Read the kernel information from file "filename"
97  CALL open_file(file_name=filename, unit_number=funit, file_form="FORMATTED")
98  READ (funit, *) nqs, nr_points
99  READ (funit, *) dispersion_env%r_max
100  END IF
101  CALL para_env%bcast(nqs)
102  CALL para_env%bcast(nr_points)
103  CALL para_env%bcast(dispersion_env%r_max)
104  ALLOCATE (dispersion_env%q_mesh(nqs), dispersion_env%kernel(0:nr_points, nqs, nqs), &
105  dispersion_env%d2phi_dk2(0:nr_points, nqs, nqs))
106  dispersion_env%nqs = nqs
107  dispersion_env%nr_points = nr_points
108  IF (para_env%is_source()) THEN
109  !! Read in the values of the q points used to generate this kernel
110  READ (funit, "(1p, 4e23.14)") dispersion_env%q_mesh
111  !! For each pair of q values, read in the function phi_q1_q2(k).
112  !! That is, the fourier transformed kernel function assuming q1 and q2
113  !! for all the values of r used.
114  DO q1_i = 1, nqs
115  DO q2_i = 1, q1_i
116  READ (funit, "(1p, 4e23.14)") dispersion_env%kernel(0:nr_points, q1_i, q2_i)
117  dispersion_env%kernel(0:nr_points, q2_i, q1_i) = dispersion_env%kernel(0:nr_points, q1_i, q2_i)
118  END DO
119  END DO
120  !! Again, for each pair of q values (q1 and q2), read in the value
121  !! of the second derivative of the above mentiond Fourier transformed
122  !! kernel function phi_alpha_beta(k). These are used for spline
123  !! interpolation of the Fourier transformed kernel.
124  DO q1_i = 1, nqs
125  DO q2_i = 1, q1_i
126  READ (funit, "(1p, 4e23.14)") dispersion_env%d2phi_dk2(0:nr_points, q1_i, q2_i)
127  dispersion_env%d2phi_dk2(0:nr_points, q2_i, q1_i) = dispersion_env%d2phi_dk2(0:nr_points, q1_i, q2_i)
128  END DO
129  END DO
130  CALL close_file(unit_number=funit)
131  END IF
132  CALL para_env%bcast(dispersion_env%q_mesh)
133  CALL para_env%bcast(dispersion_env%kernel)
134  CALL para_env%bcast(dispersion_env%d2phi_dk2)
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_rho ...
153 !> \param rho_r ...
154 !> \param rho_g ...
155 !> \param edispersion ...
156 !> \param dispersion_env ...
157 !> \param energy_only ...
158 !> \param pw_pool ...
159 !> \param xc_pw_pool ...
160 !> \param para_env ...
161 !> \param virial ...
162 ! **************************************************************************************************
163  SUBROUTINE calculate_dispersion_nonloc(vxc_rho, rho_r, rho_g, edispersion, &
164  dispersion_env, energy_only, pw_pool, xc_pw_pool, para_env, virial)
165  TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: vxc_rho, rho_r
166  TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
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, xc_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_dispersion_nonloc'
175  INTEGER, DIMENSION(3, 3), PARAMETER :: nd = reshape((/1, 0, 0, 0, 1, 0, 0, 0, 1/), (/3, 3/))
176 
177  INTEGER :: handle, i, i_grid, idir, ispin, nl_type, &
178  np, nspin, p, q, r, s
179  INTEGER, DIMENSION(1:3) :: hi, lo, n
180  LOGICAL :: use_virial
181  REAL(kind=dp) :: b_value, beta, const, ec_nl, sumnp
182  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: dq0_dgradrho, dq0_drho, hpot, potential, &
183  q0, rho
184  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: drho, thetas
185  TYPE(pw_c1d_gs_type) :: tmp_g, vxc_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, vxc_r
189  TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:, :) :: drho_r
190 
191  CALL timeset(routinen, handle)
192 
193  cpassert(ASSOCIATED(rho_r))
194  cpassert(ASSOCIATED(rho_g))
195  cpassert(ASSOCIATED(pw_pool))
196 
197  IF (PRESENT(virial)) THEN
198  use_virial = virial%pv_calculate .AND. (.NOT. virial%pv_numer)
199  ELSE
200  use_virial = .false.
201  END IF
202  IF (use_virial) THEN
203  cpassert(.NOT. energy_only)
204  END IF
205  IF (.NOT. energy_only) THEN
206  cpassert(ASSOCIATED(vxc_rho))
207  END IF
208 
209  nl_type = dispersion_env%nl_type
210 
211  b_value = dispersion_env%b_value
212  beta = 0.03125_dp*(3.0_dp/(b_value**2.0_dp))**0.75_dp
213  nspin = SIZE(rho_r)
214 
215  const = 1.0_dp/(3.0_dp*rootpi*b_value**1.5_dp)/(pi**0.75_dp)
216 
217  ! temporary arrays for FFT
218  CALL pw_pool%create_pw(tmp_g)
219  CALL pw_pool%create_pw(tmp_r)
220 
221  ! get density derivatives
222  ALLOCATE (drho_r(3, nspin))
223  DO ispin = 1, nspin
224  DO idir = 1, 3
225  CALL pw_pool%create_pw(drho_r(idir, ispin))
226  CALL pw_transfer(rho_g(ispin), tmp_g)
227  CALL pw_derive(tmp_g, nd(:, idir))
228  CALL pw_transfer(tmp_g, drho_r(idir, ispin))
229  END DO
230  END DO
231 
232  np = SIZE(tmp_r%array)
233  ALLOCATE (rho(np), drho(np, 3)) !in the following loops, rho and drho _will_ have the same bounds
234  DO i = 1, 3
235  lo(i) = lbound(tmp_r%array, i)
236  hi(i) = ubound(tmp_r%array, i)
237  n(i) = hi(i) - lo(i) + 1
238  END DO
239 !$OMP PARALLEL DO DEFAULT(NONE) &
240 !$OMP SHARED(n,rho) &
241 !$OMP COLLAPSE(3)
242  DO r = 0, n(3) - 1
243  DO q = 0, n(2) - 1
244  DO p = 0, n(1) - 1
245  rho(r*n(2)*n(1) + q*n(1) + p + 1) = 0.0_dp
246  END DO
247  END DO
248  END DO
249 !$OMP END PARALLEL DO
250  DO i = 1, 3
251 !$OMP PARALLEL DO DEFAULT(NONE) &
252 !$OMP SHARED(n,i,drho) &
253 !$OMP COLLAPSE(3)
254  DO r = 0, n(3) - 1
255  DO q = 0, n(2) - 1
256  DO p = 0, n(1) - 1
257  drho(r*n(2)*n(1) + q*n(1) + p + 1, i) = 0.0_dp
258  END DO
259  END DO
260  END DO
261 !$OMP END PARALLEL DO
262  END DO
263  DO ispin = 1, nspin
264  CALL pw_transfer(rho_g(ispin), tmp_g)
265  CALL pw_transfer(tmp_g, tmp_r)
266 !$OMP PARALLEL DO DEFAULT(NONE) &
267 !$OMP SHARED(n,lo,rho,tmp_r) &
268 !$OMP PRIVATE(s) &
269 !$OMP COLLAPSE(3)
270  DO r = 0, n(3) - 1
271  DO q = 0, n(2) - 1
272  DO p = 0, n(1) - 1
273  s = r*n(2)*n(1) + q*n(1) + p + 1
274  rho(s) = rho(s) + tmp_r%array(p + lo(1), q + lo(2), r + lo(3))
275  END DO
276  END DO
277  END DO
278 !$OMP END PARALLEL DO
279  DO i = 1, 3
280 !$OMP PARALLEL DO DEFAULT(NONE) &
281 !$OMP SHARED(ispin,i,n,lo,drho,drho_r) &
282 !$OMP PRIVATE(s) &
283 !$OMP COLLAPSE(3)
284  DO r = 0, n(3) - 1
285  DO q = 0, n(2) - 1
286  DO p = 0, n(1) - 1
287  s = r*n(2)*n(1) + q*n(1) + p + 1
288  drho(s, i) = drho(s, i) + drho_r(i, ispin)%array(p + lo(1), q + lo(2), r + lo(3))
289  END DO
290  END DO
291  END DO
292 !$OMP END PARALLEL DO
293  END DO
294  END DO
295  !! ---------------------------------------------------------------------------------
296  !! Find the value of q0 for all assigned grid points. q is defined in equations
297  !! 11 and 12 of DION and q0 is the saturated version of q defined in equation
298  !! 5 of SOLER. This routine also returns the derivatives of the q0s with respect
299  !! to the charge-density and the gradient of the charge-density. These are needed
300  !! for the potential calculated below.
301  !! ---------------------------------------------------------------------------------
302 
303  IF (energy_only) THEN
304  ALLOCATE (q0(np))
305  SELECT CASE (nl_type)
306  CASE DEFAULT
307  cpabort("Unknown vdW-DF functional")
308  CASE (vdw_nl_drsll, vdw_nl_lmkll)
309  CALL get_q0_on_grid_eo_vdw(rho, drho, q0, dispersion_env)
310  CASE (vdw_nl_rvv10)
311  CALL get_q0_on_grid_eo_rvv10(rho, drho, q0, dispersion_env)
312  END SELECT
313  ELSE
314  ALLOCATE (q0(np), dq0_drho(np), dq0_dgradrho(np))
315  SELECT CASE (nl_type)
316  CASE DEFAULT
317  cpabort("Unknown vdW-DF functional")
318  CASE (vdw_nl_drsll, vdw_nl_lmkll)
319  CALL get_q0_on_grid_vdw(rho, drho, q0, dq0_drho, dq0_dgradrho, dispersion_env)
320  CASE (vdw_nl_rvv10)
321  CALL get_q0_on_grid_rvv10(rho, drho, q0, dq0_drho, dq0_dgradrho, dispersion_env)
322  END SELECT
323  END IF
324 
325  !! ---------------------------------------------------------------------------------------------
326  !! Here we allocate and calculate the theta functions appearing in equations 8-12 of SOLER.
327  !! They are defined as rho*P_i(q0(rho, gradient_rho)) for vdW and vdW2 or
328  !! constant*rho^(3/4)*P_i(q0(rho, gradient_rho)) for rVV10 where P_i is a polynomial that
329  !! interpolates a Kronecker delta function at the point q_i (taken from the q_mesh) and q0 is
330  !! the saturated version of q.
331  !! q is defined in equations 11 and 12 of DION and the saturation procedure is defined in equation 5
332  !! of soler. This is the biggest memory consumer in the method since the thetas array is
333  !! (total # of FFT points)*Nqs complex numbers. In a parallel run, each processor will hold the
334  !! values of all the theta functions on just the points assigned to it.
335  !! --------------------------------------------------------------------------------------------------
336  !! thetas are stored in reciprocal space as theta_i(k) because this is the way they are used later
337  !! for the convolution (equation 11 of SOLER).
338  !! --------------------------------------------------------------------------------------------------
339  ALLOCATE (thetas(np, dispersion_env%nqs))
340  !! Interpolate the P_i polynomials defined in equation 3 in SOLER for the particular
341  !! q0 values we have.
342  CALL spline_interpolation(dispersion_env%q_mesh, dispersion_env%d2y_dx2, q0, thetas)
343 
344  !! Form the thetas where theta is defined as rho*p_i(q0) for vdW and vdW2 or
345  !! constant*rho^(3/4)*p_i(q0) for rVV10
346  !! ------------------------------------------------------------------------------------
347  SELECT CASE (nl_type)
348  CASE DEFAULT
349  cpabort("Unknown vdW-DF functional")
350  CASE (vdw_nl_drsll, vdw_nl_lmkll)
351 !$OMP PARALLEL DO DEFAULT( NONE ) &
352 !$OMP SHARED( dispersion_env, thetas, rho)
353  DO i = 1, dispersion_env%nqs
354  thetas(:, i) = thetas(:, i)*rho(:)
355  END DO
356 !$OMP END PARALLEL DO
357  CASE (vdw_nl_rvv10)
358 !$OMP PARALLEL DO DEFAULT( NONE ) &
359 !$OMP SHARED( np, rho, dispersion_env, thetas, const ) &
360 !$OMP PRIVATE( i ) &
361 !$OMP SCHEDULE(DYNAMIC) ! use dynamic to allow for possibility of cases having (rho(i_grid) .LE. epsr)
362  DO i_grid = 1, np
363  IF (rho(i_grid) > epsr) THEN
364  DO i = 1, dispersion_env%nqs
365  thetas(i_grid, i) = thetas(i_grid, i)*const*rho(i_grid)**(0.75_dp)
366  END DO
367  ELSE
368  thetas(i_grid, :) = 0.0_dp
369  END IF
370  END DO
371 !$OMP END PARALLEL DO
372  END SELECT
373  !! ------------------------------------------------------------------------------------
374  !! Get thetas in reciprocal space.
375  DO i = 1, 3
376  lo(i) = lbound(tmp_r%array, i)
377  hi(i) = ubound(tmp_r%array, i)
378  n(i) = hi(i) - lo(i) + 1
379  END DO
380  ALLOCATE (thetas_g(dispersion_env%nqs))
381  DO i = 1, dispersion_env%nqs
382 !$OMP PARALLEL DO DEFAULT(NONE) &
383 !$OMP SHARED(i,n,lo,thetas,tmp_r) &
384 !$OMP COLLAPSE(3)
385  DO r = 0, n(3) - 1
386  DO q = 0, n(2) - 1
387  DO p = 0, n(1) - 1
388  tmp_r%array(p + lo(1), q + lo(2), r + lo(3)) = thetas(r*n(2)*n(1) + q*n(1) + p + 1, i)
389  END DO
390  END DO
391  END DO
392 !$OMP END PARALLEL DO
393  CALL pw_pool%create_pw(thetas_g(i))
394  CALL pw_transfer(tmp_r, thetas_g(i))
395  END DO
396  grid => thetas_g(1)%pw_grid
397  !! ---------------------------------------------------------------------------------------------
398  !! Carry out the integration in equation 8 of SOLER. This also turns the thetas array into the
399  !! precursor to the u_i(k) array which is inverse fourier transformed to get the u_i(r) functions
400  !! of SOLER equation 11. Add the energy we find to the output variable etxc.
401  !! --------------------------------------------------------------------------------------------------
402  sumnp = np
403  CALL para_env%sum(sumnp)
404  IF (use_virial) THEN
405  ! calculates kernel contribution to stress
406  CALL vdw_energy(thetas_g, dispersion_env, ec_nl, energy_only, virial)
407  SELECT CASE (nl_type)
408  CASE (vdw_nl_rvv10)
409  ec_nl = 0.5_dp*ec_nl + beta*sum(rho(:))*grid%vol/sumnp
410  END SELECT
411  ! calculates energy contribution to stress
412  ! potential contribution to stress is calculated together with other potentials (Hxc)
413  DO idir = 1, 3
414  virial%pv_xc(idir, idir) = virial%pv_xc(idir, idir) + ec_nl
415  END DO
416  ELSE
417  CALL vdw_energy(thetas_g, dispersion_env, ec_nl, energy_only)
418  SELECT CASE (nl_type)
419  CASE (vdw_nl_rvv10)
420  ec_nl = 0.5_dp*ec_nl + beta*sum(rho(:))*grid%vol/sumnp
421  END SELECT
422  END IF
423  CALL para_env%sum(ec_nl)
424  IF (nl_type == vdw_nl_rvv10) ec_nl = ec_nl*dispersion_env%scale_rvv10
425  edispersion = ec_nl
426 
427  IF (energy_only) THEN
428  DEALLOCATE (q0)
429  ELSE
430  !! ----------------------------------------------------------------------------
431  !! Inverse Fourier transform the u_i(k) to get the u_i(r) of SOLER equation 11.
432  !!-----------------------------------------------------------------------------
433  DO i = 1, 3
434  lo(i) = lbound(tmp_r%array, i)
435  hi(i) = ubound(tmp_r%array, i)
436  n(i) = hi(i) - lo(i) + 1
437  END DO
438  DO i = 1, dispersion_env%nqs
439  CALL pw_transfer(thetas_g(i), tmp_r)
440 !$OMP PARALLEL DO DEFAULT(NONE) &
441 !$OMP SHARED(i,n,lo,thetas,tmp_r) &
442 !$OMP COLLAPSE(3)
443  DO r = 0, n(3) - 1
444  DO q = 0, n(2) - 1
445  DO p = 0, n(1) - 1
446  thetas(r*n(2)*n(1) + q*n(1) + p + 1, i) = tmp_r%array(p + lo(1), q + lo(2), r + lo(3))
447  END DO
448  END DO
449  END DO
450 !$OMP END PARALLEL DO
451  END DO
452  !! -------------------------------------------------------------------------
453  !! Here we allocate the array to hold the potential. This is calculated via
454  !! equation 10 of SOLER, using the u_i(r) calculated from equations 11 and
455  !! 12 of SOLER. Each processor allocates the array to be the size of the
456  !! full grid because, as can be seen in SOLER equation 9, processors need
457  !! to access grid points outside their allocated regions.
458  !! -------------------------------------------------------------------------
459  ALLOCATE (potential(np), hpot(np))
460  IF (use_virial) THEN
461  ! calculates gradient contribution to stress
462  grid => tmp_g%pw_grid
463  CALL get_potential(q0, dq0_drho, dq0_dgradrho, rho, thetas, potential, hpot, &
464  dispersion_env, drho, grid%dvol, virial)
465  ELSE
466  CALL get_potential(q0, dq0_drho, dq0_dgradrho, rho, thetas, potential, hpot, &
467  dispersion_env)
468  END IF
469  SELECT CASE (nl_type)
470  CASE (vdw_nl_rvv10)
471  potential(:) = (0.5_dp*potential(:) + beta)*dispersion_env%scale_rvv10
472  hpot(:) = 0.5_dp*dispersion_env%scale_rvv10*hpot(:)
473  END SELECT
474  CALL pw_pool%create_pw(vxc_r)
475  DO i = 1, 3
476  lo(i) = lbound(vxc_r%array, i)
477  hi(i) = ubound(vxc_r%array, i)
478  n(i) = hi(i) - lo(i) + 1
479  END DO
480 !$OMP PARALLEL DO DEFAULT(NONE) &
481 !$OMP SHARED(i,n,lo,potential,vxc_r) &
482 !$OMP COLLAPSE(3)
483  DO r = 0, n(3) - 1
484  DO q = 0, n(2) - 1
485  DO p = 0, n(1) - 1
486  vxc_r%array(p + lo(1), q + lo(2), r + lo(3)) = potential(r*n(2)*n(1) + q*n(1) + p + 1)
487  END DO
488  END DO
489  END DO
490 !$OMP END PARALLEL DO
491  DO i = 1, 3
492  lo(i) = lbound(tmp_r%array, i)
493  hi(i) = ubound(tmp_r%array, i)
494  n(i) = hi(i) - lo(i) + 1
495  END DO
496  DO idir = 1, 3
497 !$OMP PARALLEL DO DEFAULT(NONE) &
498 !$OMP SHARED(n,lo,tmp_r) &
499 !$OMP COLLAPSE(3)
500  DO r = 0, n(3) - 1
501  DO q = 0, n(2) - 1
502  DO p = 0, n(1) - 1
503  tmp_r%array(p + lo(1), q + lo(2), r + lo(3)) = 0.0_dp
504  END DO
505  END DO
506  END DO
507 !$OMP END PARALLEL DO
508  DO ispin = 1, nspin
509 !$OMP PARALLEL DO DEFAULT(NONE) &
510 !$OMP SHARED(n,lo,tmp_r,hpot,drho_r,idir,ispin) &
511 !$OMP COLLAPSE(3)
512  DO r = 0, n(3) - 1
513  DO q = 0, n(2) - 1
514  DO p = 0, n(1) - 1
515  tmp_r%array(p + lo(1), q + lo(2), r + lo(3)) = tmp_r%array(p + lo(1), q + lo(2), r + lo(3)) &
516  + hpot(r*n(2)*n(1) + q*n(1) + p + 1) &
517  *drho_r(idir, ispin)%array(p + lo(1), q + lo(2), r + lo(3))
518  END DO
519  END DO
520  END DO
521 !$OMP END PARALLEL DO
522  END DO
523  CALL pw_transfer(tmp_r, tmp_g)
524  CALL pw_derive(tmp_g, nd(:, idir))
525  CALL pw_transfer(tmp_g, tmp_r)
526  CALL pw_axpy(tmp_r, vxc_r, -1._dp)
527  END DO
528  CALL pw_transfer(vxc_r, tmp_g)
529  CALL pw_pool%give_back_pw(vxc_r)
530  CALL xc_pw_pool%create_pw(vxc_r)
531  CALL xc_pw_pool%create_pw(vxc_g)
532  CALL pw_transfer(tmp_g, vxc_g)
533  CALL pw_transfer(vxc_g, vxc_r)
534  DO ispin = 1, nspin
535  CALL pw_axpy(vxc_r, vxc_rho(ispin), 1._dp)
536  END DO
537  CALL xc_pw_pool%give_back_pw(vxc_r)
538  CALL xc_pw_pool%give_back_pw(vxc_g)
539  DEALLOCATE (q0, dq0_drho, dq0_dgradrho)
540  END IF
541 
542  DEALLOCATE (thetas)
543 
544  DO i = 1, dispersion_env%nqs
545  CALL pw_pool%give_back_pw(thetas_g(i))
546  END DO
547  DO ispin = 1, nspin
548  DO idir = 1, 3
549  CALL pw_pool%give_back_pw(drho_r(idir, ispin))
550  END DO
551  END DO
552  CALL pw_pool%give_back_pw(tmp_r)
553  CALL pw_pool%give_back_pw(tmp_g)
554 
555  DEALLOCATE (rho, drho, drho_r, thetas_g)
556 
557  CALL timestop(handle)
558 
559  END SUBROUTINE calculate_dispersion_nonloc
560 
561 ! **************************************************************************************************
562 !> \brief This routine carries out the integration of equation 8 of SOLER. It returns the non-local
563 !> exchange-correlation energy and the u_alpha(k) arrays used to find the u_alpha(r) arrays via
564 !> equations 11 and 12 in SOLER.
565 !> energy contribution to stress is added in qs_force
566 !> \param thetas_g ...
567 !> \param dispersion_env ...
568 !> \param vdW_xc_energy ...
569 !> \param energy_only ...
570 !> \param virial ...
571 !> \par History
572 !> OpenMP added: Aug 2016 MTucker
573 ! **************************************************************************************************
574  SUBROUTINE vdw_energy(thetas_g, dispersion_env, vdW_xc_energy, energy_only, virial)
575  TYPE(pw_c1d_gs_type), DIMENSION(:), INTENT(IN) :: thetas_g
576  TYPE(qs_dispersion_type), POINTER :: dispersion_env
577  REAL(kind=dp), INTENT(OUT) :: vdw_xc_energy
578  LOGICAL, INTENT(IN) :: energy_only
579  TYPE(virial_type), OPTIONAL, POINTER :: virial
580 
581  CHARACTER(LEN=*), PARAMETER :: routinen = 'vdW_energy'
582 
583  COMPLEX(KIND=dp) :: uu
584  COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: theta
585  COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: u_vdw(:, :)
586  INTEGER :: handle, ig, iq, l, m, nl_type, nqs, &
587  q1_i, q2_i
588  LOGICAL :: use_virial
589  REAL(kind=dp) :: g, g2, g2_last, g_multiplier, gm
590  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: dkernel_of_dk, kernel_of_k
591  REAL(kind=dp), DIMENSION(3, 3) :: virial_thread
592  TYPE(pw_grid_type), POINTER :: grid
593 
594  CALL timeset(routinen, handle)
595  nqs = dispersion_env%nqs
596 
597  use_virial = PRESENT(virial)
598  virial_thread(:, :) = 0.0_dp ! always initialize to avoid floating point exceptions in OMP REDUCTION
599 
600  vdw_xc_energy = 0._dp
601  grid => thetas_g(1)%pw_grid
602 
603  IF (grid%grid_span == halfspace) THEN
604  g_multiplier = 2._dp
605  ELSE
606  g_multiplier = 1._dp
607  END IF
608 
609  nl_type = dispersion_env%nl_type
610 
611  IF (.NOT. energy_only) THEN
612  ALLOCATE (u_vdw(grid%ngpts_cut_local, nqs))
613  u_vdw(:, :) = cmplx(0.0_dp, 0.0_dp, kind=dp)
614  END IF
615 
616 !$OMP PARALLEL DEFAULT( NONE ) &
617 !$OMP SHARED( nqs, energy_only, grid, dispersion_env &
618 !$OMP , use_virial, thetas_g, g_multiplier, nl_type &
619 !$OMP , u_vdW &
620 !$OMP ) &
621 !$OMP PRIVATE( g2_last, kernel_of_k, dkernel_of_dk, theta &
622 !$OMP , g2, g, iq &
623 !$OMP , q2_i, uu, q1_i, gm, l, m &
624 !$OMP ) &
625 !$OMP REDUCTION( +: vdW_xc_energy, virial_thread &
626 !$OMP )
627 
628  g2_last = huge(0._dp)
629 
630  ALLOCATE (kernel_of_k(nqs, nqs), dkernel_of_dk(nqs, nqs))
631  ALLOCATE (theta(nqs))
632 
633 !$OMP DO
634  DO ig = 1, grid%ngpts_cut_local
635  g2 = grid%gsq(ig)
636  IF (abs(g2 - g2_last) > 1.e-10) THEN
637  g2_last = g2
638  g = sqrt(g2)
639  CALL interpolate_kernel(g, kernel_of_k, dispersion_env)
640  IF (use_virial) CALL interpolate_dkernel_dk(g, dkernel_of_dk, dispersion_env)
641  END IF
642  DO iq = 1, nqs
643  theta(iq) = thetas_g(iq)%array(ig)
644  END DO
645  DO q2_i = 1, nqs
646  uu = cmplx(0.0_dp, 0.0_dp, kind=dp)
647  DO q1_i = 1, nqs
648  uu = uu + kernel_of_k(q2_i, q1_i)*theta(q1_i)
649  END DO
650  IF (ig < grid%first_gne0) THEN
651  vdw_xc_energy = vdw_xc_energy + real(uu*conjg(theta(q2_i)), kind=dp)
652  ELSE
653  vdw_xc_energy = vdw_xc_energy + g_multiplier*real(uu*conjg(theta(q2_i)), kind=dp)
654  END IF
655  IF (.NOT. energy_only) u_vdw(ig, q2_i) = uu
656 
657  IF (use_virial .AND. ig >= grid%first_gne0) THEN
658  DO q1_i = 1, nqs
659  gm = 0.5_dp*g_multiplier*grid%vol*dkernel_of_dk(q1_i, q2_i)*real(theta(q1_i)*conjg(theta(q2_i)), kind=dp)
660  IF (nl_type == vdw_nl_rvv10) THEN
661  gm = 0.5_dp*gm
662  END IF
663  DO l = 1, 3
664  DO m = 1, l
665  virial_thread(l, m) = virial_thread(l, m) - gm*(grid%g(l, ig)*grid%g(m, ig))/g
666  END DO
667  END DO
668  END DO
669  END IF
670 
671  END DO
672  END DO
673 !$OMP END DO
674 
675  DEALLOCATE (theta)
676  DEALLOCATE (kernel_of_k, dkernel_of_dk)
677 
678  IF (.NOT. energy_only) THEN
679 !$OMP DO
680  DO ig = 1, grid%ngpts_cut_local
681  DO iq = 1, nqs
682  thetas_g(iq)%array(ig) = u_vdw(ig, iq)
683  END DO
684  END DO
685 !$OMP END DO
686  END IF
687 
688 !$OMP END PARALLEL
689 
690  IF (.NOT. energy_only) THEN
691  DEALLOCATE (u_vdw)
692  END IF
693 
694  vdw_xc_energy = vdw_xc_energy*grid%vol*0.5_dp
695 
696  IF (use_virial) THEN
697  DO l = 1, 3
698  DO m = 1, (l - 1)
699  virial%pv_xc(l, m) = virial%pv_xc(l, m) + virial_thread(l, m)
700  virial%pv_xc(m, l) = virial%pv_xc(l, m)
701  END DO
702  m = l
703  virial%pv_xc(l, m) = virial%pv_xc(l, m) + virial_thread(l, m)
704  END DO
705  END IF
706 
707  CALL timestop(handle)
708 
709  END SUBROUTINE vdw_energy
710 
711 ! **************************************************************************************************
712 !> \brief This routine finds the non-local correlation contribution to the potential
713 !> (i.e. the derivative of the non-local piece of the energy with respect to
714 !> density) given in SOLER equation 10. The u_alpha(k) functions were found
715 !> while calculating the energy. They are passed in as the matrix u_vdW.
716 !> Most of the required derivatives were calculated in the "get_q0_on_grid"
717 !> routine, but the derivative of the interpolation polynomials, P_alpha(q),
718 !> (SOLER equation 3) with respect to q is interpolated here, along with the
719 !> polynomials themselves.
720 !> \param q0 ...
721 !> \param dq0_drho ...
722 !> \param dq0_dgradrho ...
723 !> \param total_rho ...
724 !> \param u_vdW ...
725 !> \param potential ...
726 !> \param h_prefactor ...
727 !> \param dispersion_env ...
728 !> \param drho ...
729 !> \param dvol ...
730 !> \param virial ...
731 !> \par History
732 !> OpenMP added: Aug 2016 MTucker
733 ! **************************************************************************************************
734  SUBROUTINE get_potential(q0, dq0_drho, dq0_dgradrho, total_rho, u_vdW, potential, h_prefactor, &
735  dispersion_env, drho, dvol, virial)
736 
737  REAL(dp), DIMENSION(:), INTENT(in) :: q0, dq0_drho, dq0_dgradrho, total_rho
738  REAL(dp), DIMENSION(:, :), INTENT(in) :: u_vdw
739  REAL(dp), DIMENSION(:), INTENT(out) :: potential, h_prefactor
740  TYPE(qs_dispersion_type), POINTER :: dispersion_env
741  REAL(dp), DIMENSION(:, :), INTENT(in), OPTIONAL :: drho
742  REAL(dp), INTENT(IN), OPTIONAL :: dvol
743  TYPE(virial_type), OPTIONAL, POINTER :: virial
744 
745  CHARACTER(len=*), PARAMETER :: routinen = 'get_potential'
746 
747  INTEGER :: handle, i_grid, l, m, nl_type, nqs, p_i, &
748  q, q_hi, q_low
749  LOGICAL :: use_virial
750  REAL(dp) :: a, b, b_value, c, const, d, dp_dq0, dq, &
751  dq_6, e, f, p, prefactor, tmp_1_2, &
752  tmp_1_4, tmp_3_4
753  REAL(dp), ALLOCATABLE, DIMENSION(:) :: y
754  REAL(dp), DIMENSION(3, 3) :: virial_thread
755  REAL(dp), DIMENSION(:), POINTER :: q_mesh
756  REAL(dp), DIMENSION(:, :), POINTER :: d2y_dx2
757 
758  CALL timeset(routinen, handle)
759 
760  use_virial = PRESENT(virial)
761  cpassert(.NOT. use_virial .OR. PRESENT(drho))
762  cpassert(.NOT. use_virial .OR. PRESENT(dvol))
763 
764  virial_thread(:, :) = 0.0_dp ! always initialize to avoid floating point exceptions in OMP REDUCTION
765  b_value = dispersion_env%b_value
766  const = 1.0_dp/(3.0_dp*b_value**(3.0_dp/2.0_dp)*pi**(5.0_dp/4.0_dp))
767  potential = 0.0_dp
768  h_prefactor = 0.0_dp
769 
770  d2y_dx2 => dispersion_env%d2y_dx2
771  q_mesh => dispersion_env%q_mesh
772  nqs = dispersion_env%nqs
773  nl_type = dispersion_env%nl_type
774 
775 !$OMP PARALLEL DEFAULT( NONE ) &
776 !$OMP SHARED( nqs, u_vdW, q_mesh, q0, d2y_dx2, nl_type &
777 !$OMP , potential, h_prefactor &
778 !$OMP , dq0_drho, dq0_dgradrho, total_rho, const &
779 !$OMP , use_virial, drho, dvol, virial &
780 !$OMP ) &
781 !$OMP PRIVATE( y &
782 !$OMP , q_low, q_hi, q, dq, dq_6, A, b, c, d, e, f &
783 !$OMP , P_i, dP_dq0, P, prefactor, l, m &
784 !$OMP , tmp_1_2, tmp_1_4, tmp_3_4 &
785 !$OMP ) &
786 !$OMP REDUCTION(+: virial_thread &
787 !$OMP )
788 
789  ALLOCATE (y(nqs))
790 
791 !$OMP DO
792  DO i_grid = 1, SIZE(u_vdw, 1)
793  q_low = 1
794  q_hi = nqs
795  ! Figure out which bin our value of q0 is in in the q_mesh
796  DO WHILE ((q_hi - q_low) > 1)
797  q = int((q_hi + q_low)/2)
798  IF (q_mesh(q) > q0(i_grid)) THEN
799  q_hi = q
800  ELSE
801  q_low = q
802  END IF
803  END DO
804  IF (q_hi == q_low) cpabort("get_potential: qhi == qlow")
805 
806  dq = q_mesh(q_hi) - q_mesh(q_low)
807  dq_6 = dq/6.0_dp
808 
809  a = (q_mesh(q_hi) - q0(i_grid))/dq
810  b = (q0(i_grid) - q_mesh(q_low))/dq
811  c = (a**3 - a)*dq*dq_6
812  d = (b**3 - b)*dq*dq_6
813  e = (3.0_dp*a**2 - 1.0_dp)*dq_6
814  f = (3.0_dp*b**2 - 1.0_dp)*dq_6
815 
816  DO p_i = 1, nqs
817  y = 0.0_dp
818  y(p_i) = 1.0_dp
819  dp_dq0 = (y(q_hi) - y(q_low))/dq - e*d2y_dx2(p_i, q_low) + f*d2y_dx2(p_i, q_hi)
820  p = a*y(q_low) + b*y(q_hi) + c*d2y_dx2(p_i, q_low) + d*d2y_dx2(p_i, q_hi)
821  !! The first term in equation 13 of SOLER
822  SELECT CASE (nl_type)
823  CASE DEFAULT
824  cpabort("Unknown vdW-DF functional")
825  CASE (vdw_nl_drsll, vdw_nl_lmkll)
826  potential(i_grid) = potential(i_grid) + u_vdw(i_grid, p_i)*(p + dp_dq0*dq0_drho(i_grid))
827  prefactor = u_vdw(i_grid, p_i)*dp_dq0*dq0_dgradrho(i_grid)
828  CASE (vdw_nl_rvv10)
829  IF (total_rho(i_grid) > epsr) THEN
830  tmp_1_2 = sqrt(total_rho(i_grid))
831  tmp_1_4 = sqrt(tmp_1_2) ! == total_rho(i_grid)**(1.0_dp/4.0_dp)
832  tmp_3_4 = tmp_1_4*tmp_1_4*tmp_1_4 ! == total_rho(i_grid)**(3.0_dp/4.0_dp)
833  potential(i_grid) = potential(i_grid) &
834  + u_vdw(i_grid, p_i)*(const*0.75_dp/tmp_1_4*p &
835  + const*tmp_3_4*dp_dq0*dq0_drho(i_grid))
836  prefactor = u_vdw(i_grid, p_i)*const*tmp_3_4*dp_dq0*dq0_dgradrho(i_grid)
837  ELSE
838  prefactor = 0.0_dp
839  END IF
840  END SELECT
841  IF (q0(i_grid) .NE. q_mesh(nqs)) THEN
842  h_prefactor(i_grid) = h_prefactor(i_grid) + prefactor
843  END IF
844 
845  IF (use_virial .AND. abs(prefactor) > 0.0_dp) THEN
846  SELECT CASE (nl_type)
847  CASE DEFAULT
848  ! do nothing
849  CASE (vdw_nl_rvv10)
850  prefactor = 0.5_dp*prefactor
851  END SELECT
852  prefactor = prefactor*dvol
853  DO l = 1, 3
854  DO m = 1, l
855  virial_thread(l, m) = virial_thread(l, m) - prefactor*drho(i_grid, l)*drho(i_grid, m)
856  END DO
857  END DO
858  END IF
859 
860  END DO ! P_i = 1, nqs
861  END DO ! i_grid = 1, SIZE(u_vdW, 1)
862 !$OMP END DO
863 
864  DEALLOCATE (y)
865 !$OMP END PARALLEL
866 
867  IF (use_virial) THEN
868  DO l = 1, 3
869  DO m = 1, (l - 1)
870  virial%pv_xc(l, m) = virial%pv_xc(l, m) + virial_thread(l, m)
871  virial%pv_xc(m, l) = virial%pv_xc(l, m)
872  END DO
873  m = l
874  virial%pv_xc(l, m) = virial%pv_xc(l, m) + virial_thread(l, m)
875  END DO
876  END IF
877 
878  CALL timestop(handle)
879  END SUBROUTINE get_potential
880 
881 ! **************************************************************************************************
882 !> \brief calculates exponent = sum(from i=1 to hi, ((alpha)**i)/i) ) without <<< calling power >>>
883 !> \param hi = upper index for sum
884 !> \param alpha ...
885 !> \param exponent = output value
886 !> \par History
887 !> Created: MTucker, Aug 2016
888 ! **************************************************************************************************
889  ELEMENTAL SUBROUTINE calculate_exponent(hi, alpha, exponent)
890  INTEGER, INTENT(in) :: hi
891  REAL(dp), INTENT(in) :: alpha
892  REAL(dp), INTENT(out) :: exponent
893 
894  INTEGER :: i
895  REAL(dp) :: multiplier
896 
897  multiplier = alpha
898  exponent = alpha
899 
900  DO i = 2, hi
901  multiplier = multiplier*alpha
902  exponent = exponent + (multiplier/i)
903  END DO
904  END SUBROUTINE calculate_exponent
905 
906 ! **************************************************************************************************
907 !> \brief calculate exponent = sum(from i=1 to hi, ((alpha)**i)/i) ) without calling power
908 !> also calculates derivative using similar series
909 !> \param hi = upper index for sum
910 !> \param alpha ...
911 !> \param exponent = output value
912 !> \param derivative ...
913 !> \par History
914 !> Created: MTucker, Aug 2016
915 ! **************************************************************************************************
916  ELEMENTAL SUBROUTINE calculate_exponent_derivative(hi, alpha, exponent, derivative)
917  INTEGER, INTENT(in) :: hi
918  REAL(dp), INTENT(in) :: alpha
919  REAL(dp), INTENT(out) :: exponent, derivative
920 
921  INTEGER :: i
922  REAL(dp) :: multiplier
923 
924  derivative = 0.0d0
925  multiplier = 1.0d0
926  exponent = 0.0d0
927 
928  DO i = 1, hi
929  derivative = derivative + multiplier
930  multiplier = multiplier*alpha
931  exponent = exponent + (multiplier/i)
932  END DO
933  END SUBROUTINE calculate_exponent_derivative
934 
935  !! This routine first calculates the q value defined in (DION equations 11 and 12), then
936  !! saturates it according to (SOLER equation 5).
937 ! **************************************************************************************************
938 !> \brief This routine first calculates the q value defined in (DION equations 11 and 12), then
939 !> saturates it according to (SOLER equation 5).
940 !> \param total_rho ...
941 !> \param gradient_rho ...
942 !> \param q0 ...
943 !> \param dq0_drho ...
944 !> \param dq0_dgradrho ...
945 !> \param dispersion_env ...
946 ! **************************************************************************************************
947  SUBROUTINE get_q0_on_grid_vdw(total_rho, gradient_rho, q0, dq0_drho, dq0_dgradrho, dispersion_env)
948  !!
949  !! more specifically it calculates the following
950  !!
951  !! q0(ir) = q0 as defined above
952  !! dq0_drho(ir) = total_rho * d q0 /d rho
953  !! dq0_dgradrho = total_rho / |gradient_rho| * d q0 / d |gradient_rho|
954  !!
955  REAL(dp), INTENT(IN) :: total_rho(:), gradient_rho(:, :)
956  REAL(dp), INTENT(OUT) :: q0(:), dq0_drho(:), dq0_dgradrho(:)
957  TYPE(qs_dispersion_type), POINTER :: dispersion_env
958 
959  INTEGER, PARAMETER :: m_cut = 12
960  REAL(dp), PARAMETER :: lda_a = 0.031091_dp, lda_a1 = 0.2137_dp, lda_b1 = 7.5957_dp, &
961  lda_b2 = 3.5876_dp, lda_b3 = 1.6382_dp, lda_b4 = 0.49294_dp
962 
963  INTEGER :: i_grid
964  REAL(dp) :: dq0_dq, exponent, gradient_correction, &
965  kf, lda_1, lda_2, q, q__q_cut, q_cut, &
966  q_min, r_s, sqrt_r_s, z_ab
967 
968  q_cut = dispersion_env%q_cut
969  q_min = dispersion_env%q_min
970  SELECT CASE (dispersion_env%nl_type)
971  CASE DEFAULT
972  cpabort("Unknown vdW-DF functional")
973  CASE (vdw_nl_drsll)
974  z_ab = -0.8491_dp
975  CASE (vdw_nl_lmkll)
976  z_ab = -1.887_dp
977  END SELECT
978 
979  ! initialize q0-related arrays ...
980  q0(:) = q_cut
981  dq0_drho(:) = 0.0_dp
982  dq0_dgradrho(:) = 0.0_dp
983 
984  DO i_grid = 1, SIZE(total_rho)
985 
986  !! This prevents numerical problems. If the charge density is negative (an
987  !! unphysical situation), we simply treat it as very small. In that case,
988  !! q0 will be very large and will be saturated. For a saturated q0 the derivative
989  !! dq0_dq will be 0 so we set q0 = q_cut and dq0_drho = dq0_dgradrho = 0 and go on
990  !! to the next point.
991  !! ------------------------------------------------------------------------------------
992  IF (total_rho(i_grid) < epsr) cycle
993  !! ------------------------------------------------------------------------------------
994  !! Calculate some intermediate values needed to find q
995  !! ------------------------------------------------------------------------------------
996  kf = (3.0_dp*pi*pi*total_rho(i_grid))**(1.0_dp/3.0_dp)
997  r_s = (3.0_dp/(4.0_dp*pi*total_rho(i_grid)))**(1.0_dp/3.0_dp)
998  sqrt_r_s = sqrt(r_s)
999 
1000  gradient_correction = -z_ab/(36.0_dp*kf*total_rho(i_grid)**2) &
1001  *(gradient_rho(i_grid, 1)**2 + gradient_rho(i_grid, 2)**2 + gradient_rho(i_grid, 3)**2)
1002 
1003  lda_1 = 8.0_dp*pi/3.0_dp*(lda_a*(1.0_dp + lda_a1*r_s))
1004  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)
1005  !! ---------------------------------------------------------------
1006  !! This is the q value defined in equations 11 and 12 of DION
1007  !! ---------------------------------------------------------------
1008  q = kf + lda_1*log(1.0_dp + 1.0_dp/lda_2) + gradient_correction
1009  !! ---------------------------------------------------------------
1010  !! Here, we calculate q0 by saturating q according to equation 5 of SOLER. Also, we find
1011  !! the derivative dq0_dq needed for the derivatives dq0_drho and dq0_dgradrh0 discussed below.
1012  !! ---------------------------------------------------------------------------------------
1013  q__q_cut = q/q_cut
1014  CALL calculate_exponent_derivative(m_cut, q__q_cut, exponent, dq0_dq)
1015  q0(i_grid) = q_cut*(1.0_dp - exp(-exponent))
1016  dq0_dq = dq0_dq*exp(-exponent)
1017  !! ---------------------------------------------------------------------------------------
1018  !! This is to handle a case with q0 too small. We simply set it to the smallest q value in
1019  !! out q_mesh. Hopefully this doesn't get used often (ever)
1020  !! ---------------------------------------------------------------------------------------
1021  IF (q0(i_grid) < q_min) THEN
1022  q0(i_grid) = q_min
1023  END IF
1024  !! ---------------------------------------------------------------------------------------
1025  !! Here we find derivatives. These are actually the density times the derivative of q0 with respect
1026  !! to rho and gradient_rho. The density factor comes in since we are really differentiating
1027  !! theta = (rho)*P(q0) with respect to density (or its gradient) which will be
1028  !! dtheta_drho = P(q0) + dP_dq0 * [rho * dq0_dq * dq_drho] and
1029  !! dtheta_dgradient_rho = dP_dq0 * [rho * dq0_dq * dq_dgradient_rho]
1030  !! The parts in square brackets are what is calculated here. The dP_dq0 term will be interpolated
1031  !! later. There should actually be a factor of the magnitude of the gradient in the gradient_rho derivative
1032  !! but that cancels out when we differentiate the magnitude of the gradient with respect to a particular
1033  !! component.
1034  !! ------------------------------------------------------------------------------------------------
1035 
1036  dq0_drho(i_grid) = dq0_dq*(kf/3.0_dp - 7.0_dp/3.0_dp*gradient_correction &
1037  - 8.0_dp*pi/9.0_dp*lda_a*lda_a1*r_s*log(1.0_dp + 1.0_dp/lda_2) &
1038  + lda_1/(lda_2*(1.0_dp + lda_2)) &
1039  *(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 &
1040  + 2.0_dp*lda_b4/3.0_dp*r_s**2)))
1041 
1042  dq0_dgradrho(i_grid) = total_rho(i_grid)*dq0_dq*2.0_dp*(-z_ab)/(36.0_dp*kf*total_rho(i_grid)**2)
1043 
1044  END DO
1045 
1046  END SUBROUTINE get_q0_on_grid_vdw
1047 
1048 ! **************************************************************************************************
1049 !> \brief ...
1050 !> \param total_rho ...
1051 !> \param gradient_rho ...
1052 !> \param q0 ...
1053 !> \param dq0_drho ...
1054 !> \param dq0_dgradrho ...
1055 !> \param dispersion_env ...
1056 ! **************************************************************************************************
1057  PURE SUBROUTINE get_q0_on_grid_rvv10(total_rho, gradient_rho, q0, dq0_drho, dq0_dgradrho, dispersion_env)
1058  !!
1059  !! more specifically it calculates the following
1060  !!
1061  !! q0(ir) = q0 as defined above
1062  !! dq0_drho(ir) = total_rho * d q0 /d rho
1063  !! dq0_dgradrho = total_rho / |gradient_rho| * d q0 / d |gradient_rho|
1064  !!
1065  REAL(dp), INTENT(IN) :: total_rho(:), gradient_rho(:, :)
1066  REAL(dp), INTENT(OUT) :: q0(:), dq0_drho(:), dq0_dgradrho(:)
1067  TYPE(qs_dispersion_type), POINTER :: dispersion_env
1068 
1069  INTEGER, PARAMETER :: m_cut = 12
1070 
1071  INTEGER :: i_grid
1072  REAL(dp) :: b_value, c_value, dk_dn, dq0_dq, dw0_dn, &
1073  exponent, gmod2, k, mod_grad, q, &
1074  q__q_cut, q_cut, q_min, w0, wg2, wp2
1075 
1076  q_cut = dispersion_env%q_cut
1077  q_min = dispersion_env%q_min
1078  b_value = dispersion_env%b_value
1079  c_value = dispersion_env%c_value
1080 
1081  ! initialize q0-related arrays ...
1082  q0(:) = q_cut
1083  dq0_drho(:) = 0.0_dp
1084  dq0_dgradrho(:) = 0.0_dp
1085 
1086  DO i_grid = 1, SIZE(total_rho)
1087 
1088  gmod2 = gradient_rho(i_grid, 1)**2 + gradient_rho(i_grid, 2)**2 + gradient_rho(i_grid, 3)**2
1089 
1090  !if (total_rho(i_grid) > epsr .and. gmod2 > epsr) cycle
1091  IF (total_rho(i_grid) > epsr) THEN
1092 
1093  !! Calculate some intermediate values needed to find q
1094  !! ------------------------------------------------------------------------------------
1095  mod_grad = sqrt(gmod2)
1096 
1097  wp2 = 16.0_dp*pi*total_rho(i_grid)
1098  wg2 = 4_dp*c_value*(mod_grad/total_rho(i_grid))**4
1099 
1100  k = b_value*3.0_dp*pi*((total_rho(i_grid)/(9.0_dp*pi))**(1.0_dp/6.0_dp))
1101  w0 = sqrt(wg2 + wp2/3.0_dp)
1102 
1103  q = w0/k
1104 
1105  !! Here, we calculate q0 by saturating q according
1106  !! ---------------------------------------------------------------------------------------
1107  q__q_cut = q/q_cut
1108  CALL calculate_exponent_derivative(m_cut, q__q_cut, exponent, dq0_dq)
1109  q0(i_grid) = q_cut*(1.0_dp - exp(-exponent))
1110  dq0_dq = dq0_dq*exp(-exponent)
1111 
1112  !! ---------------------------------------------------------------------------------------
1113  IF (q0(i_grid) < q_min) THEN
1114  q0(i_grid) = q_min
1115  END IF
1116 
1117  !!---------------------------------Final values---------------------------------
1118  dw0_dn = 1.0_dp/(2.0_dp*w0)*(16.0_dp/3.0_dp*pi - 4.0_dp*wg2/total_rho(i_grid))
1119  dk_dn = k/(6.0_dp*total_rho(i_grid))
1120 
1121  dq0_drho(i_grid) = dq0_dq*1.0_dp/(k**2.0)*(dw0_dn*k - dk_dn*w0)
1122  dq0_dgradrho(i_grid) = dq0_dq*1.0_dp/(2.0_dp*k*w0)*4.0_dp*wg2/gmod2
1123  END IF
1124 
1125  END DO
1126 
1127  END SUBROUTINE get_q0_on_grid_rvv10
1128 
1129 ! **************************************************************************************************
1130 !> \brief ...
1131 !> \param total_rho ...
1132 !> \param gradient_rho ...
1133 !> \param q0 ...
1134 !> \param dispersion_env ...
1135 ! **************************************************************************************************
1136  SUBROUTINE get_q0_on_grid_eo_vdw(total_rho, gradient_rho, q0, dispersion_env)
1137 
1138  REAL(dp), INTENT(IN) :: total_rho(:), gradient_rho(:, :)
1139  REAL(dp), INTENT(OUT) :: q0(:)
1140  TYPE(qs_dispersion_type), POINTER :: dispersion_env
1141 
1142  INTEGER, PARAMETER :: m_cut = 12
1143  REAL(dp), PARAMETER :: lda_a = 0.031091_dp, lda_a1 = 0.2137_dp, lda_b1 = 7.5957_dp, &
1144  lda_b2 = 3.5876_dp, lda_b3 = 1.6382_dp, lda_b4 = 0.49294_dp
1145 
1146  INTEGER :: i_grid
1147  REAL(dp) :: exponent, gradient_correction, kf, &
1148  lda_1, lda_2, q, q__q_cut, q_cut, &
1149  q_min, r_s, sqrt_r_s, z_ab
1150 
1151  q_cut = dispersion_env%q_cut
1152  q_min = dispersion_env%q_min
1153  SELECT CASE (dispersion_env%nl_type)
1154  CASE DEFAULT
1155  cpabort("Unknown vdW-DF functional")
1156  CASE (vdw_nl_drsll)
1157  z_ab = -0.8491_dp
1158  CASE (vdw_nl_lmkll)
1159  z_ab = -1.887_dp
1160  END SELECT
1161 
1162  ! initialize q0-related arrays ...
1163  q0(:) = q_cut
1164  DO i_grid = 1, SIZE(total_rho)
1165  !! This prevents numerical problems. If the charge density is negative (an
1166  !! unphysical situation), we simply treat it as very small. In that case,
1167  !! q0 will be very large and will be saturated. For a saturated q0 the derivative
1168  !! dq0_dq will be 0 so we set q0 = q_cut and dq0_drho = dq0_dgradrho = 0 and go on
1169  !! to the next point.
1170  !! ------------------------------------------------------------------------------------
1171  IF (total_rho(i_grid) < epsr) cycle
1172  !! ------------------------------------------------------------------------------------
1173  !! Calculate some intermediate values needed to find q
1174  !! ------------------------------------------------------------------------------------
1175  kf = (3.0_dp*pi*pi*total_rho(i_grid))**(1.0_dp/3.0_dp)
1176  r_s = (3.0_dp/(4.0_dp*pi*total_rho(i_grid)))**(1.0_dp/3.0_dp)
1177  sqrt_r_s = sqrt(r_s)
1178 
1179  gradient_correction = -z_ab/(36.0_dp*kf*total_rho(i_grid)**2) &
1180  *(gradient_rho(i_grid, 1)**2 + gradient_rho(i_grid, 2)**2 + gradient_rho(i_grid, 3)**2)
1181 
1182  lda_1 = 8.0_dp*pi/3.0_dp*(lda_a*(1.0_dp + lda_a1*r_s))
1183  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)
1184  !! ------------------------------------------------------------------------------------
1185  !! This is the q value defined in equations 11 and 12 of DION
1186  !! ---------------------------------------------------------------
1187  q = kf + lda_1*log(1.0_dp + 1.0_dp/lda_2) + gradient_correction
1188 
1189  !! ---------------------------------------------------------------
1190  !! Here, we calculate q0 by saturating q according to equation 5 of SOLER. Also, we find
1191  !! the derivative dq0_dq needed for the derivatives dq0_drho and dq0_dgradrh0 discussed below.
1192  !! ---------------------------------------------------------------------------------------
1193  q__q_cut = q/q_cut
1194  CALL calculate_exponent(m_cut, q__q_cut, exponent)
1195  q0(i_grid) = q_cut*(1.0_dp - exp(-exponent))
1196 
1197  !! ---------------------------------------------------------------------------------------
1198  !! This is to handle a case with q0 too small. We simply set it to the smallest q value in
1199  !! out q_mesh. Hopefully this doesn't get used often (ever)
1200  !! ---------------------------------------------------------------------------------------
1201  IF (q0(i_grid) < q_min) THEN
1202  q0(i_grid) = q_min
1203  END IF
1204  END DO
1205 
1206  END SUBROUTINE get_q0_on_grid_eo_vdw
1207 
1208 ! **************************************************************************************************
1209 !> \brief ...
1210 !> \param total_rho ...
1211 !> \param gradient_rho ...
1212 !> \param q0 ...
1213 !> \param dispersion_env ...
1214 ! **************************************************************************************************
1215  PURE SUBROUTINE get_q0_on_grid_eo_rvv10(total_rho, gradient_rho, q0, dispersion_env)
1216 
1217  REAL(dp), INTENT(IN) :: total_rho(:), gradient_rho(:, :)
1218  REAL(dp), INTENT(OUT) :: q0(:)
1219  TYPE(qs_dispersion_type), POINTER :: dispersion_env
1220 
1221  INTEGER, PARAMETER :: m_cut = 12
1222 
1223  INTEGER :: i_grid
1224  REAL(dp) :: b_value, c_value, exponent, gmod2, k, q, &
1225  q__q_cut, q_cut, q_min, w0, wg2, wp2
1226 
1227  q_cut = dispersion_env%q_cut
1228  q_min = dispersion_env%q_min
1229  b_value = dispersion_env%b_value
1230  c_value = dispersion_env%c_value
1231 
1232  ! initialize q0-related arrays ...
1233  q0(:) = q_cut
1234 
1235  DO i_grid = 1, SIZE(total_rho)
1236 
1237  gmod2 = gradient_rho(i_grid, 1)**2 + gradient_rho(i_grid, 2)**2 + gradient_rho(i_grid, 3)**2
1238 
1239  !if (total_rho(i_grid) > epsr .and. gmod2 > epsr) cycle
1240  IF (total_rho(i_grid) > epsr) THEN
1241 
1242  !! Calculate some intermediate values needed to find q
1243  !! ------------------------------------------------------------------------------------
1244  wp2 = 16.0_dp*pi*total_rho(i_grid)
1245  wg2 = 4_dp*c_value*(gmod2*gmod2)/(total_rho(i_grid)**4)
1246 
1247  k = b_value*3.0_dp*pi*((total_rho(i_grid)/(9.0_dp*pi))**(1.0_dp/6.0_dp))
1248  w0 = sqrt(wg2 + wp2/3.0_dp)
1249 
1250  q = w0/k
1251 
1252  !! Here, we calculate q0 by saturating q according
1253  !! ---------------------------------------------------------------------------------------
1254  q__q_cut = q/q_cut
1255  CALL calculate_exponent(m_cut, q__q_cut, exponent)
1256  q0(i_grid) = q_cut*(1.0_dp - exp(-exponent))
1257 
1258  IF (q0(i_grid) < q_min) THEN
1259  q0(i_grid) = q_min
1260  END IF
1261 
1262  END IF
1263 
1264  END DO
1265 
1266  END SUBROUTINE get_q0_on_grid_eo_rvv10
1267 
1268 ! **************************************************************************************************
1269 !> \brief This routine is modeled after an algorithm from "Numerical Recipes in C" by Cambridge University
1270 !> press, page 97. It was adapted for Fortran, of course and for the problem at hand, in that it finds
1271 !> the bin a particular x value is in and then loops over all the P_i functions so we only have to find
1272 !> the bin once.
1273 !> \param x ...
1274 !> \param d2y_dx2 ...
1275 !> \param evaluation_points ...
1276 !> \param values ...
1277 !> \par History
1278 !> OpenMP added: Aug 2016 MTucker
1279 ! **************************************************************************************************
1280  SUBROUTINE spline_interpolation(x, d2y_dx2, evaluation_points, values)
1281  REAL(dp), INTENT(in) :: x(:), d2y_dx2(:, :), evaluation_points(:)
1282  REAL(dp), INTENT(out) :: values(:, :)
1283 
1284  INTEGER :: i_grid, index, lower_bound, &
1285  ngrid_points, nx, p_i, upper_bound
1286  REAL(dp) :: a, b, c, const, d, dx
1287  REAL(dp), ALLOCATABLE :: y(:)
1288 
1289  nx = SIZE(x)
1290  ngrid_points = SIZE(evaluation_points)
1291 
1292 !$OMP PARALLEL DEFAULT( NONE ) &
1293 !$OMP SHARED( x, d2y_dx2, evaluation_points, values &
1294 !$OMP , Nx, Ngrid_points &
1295 !$OMP ) &
1296 !$OMP PRIVATE( y, lower_bound, upper_bound, index &
1297 !$OMP , dx, const, A, b, c, d, P_i &
1298 !$OMP )
1299 
1300  ALLOCATE (y(nx))
1301 
1302 !$OMP DO
1303  DO i_grid = 1, ngrid_points
1304  lower_bound = 1
1305  upper_bound = nx
1306  DO WHILE ((upper_bound - lower_bound) > 1)
1307  index = (upper_bound + lower_bound)/2
1308  IF (evaluation_points(i_grid) > x(index)) THEN
1309  lower_bound = index
1310  ELSE
1311  upper_bound = index
1312  END IF
1313  END DO
1314 
1315  dx = x(upper_bound) - x(lower_bound)
1316  const = dx*dx/6.0_dp
1317 
1318  a = (x(upper_bound) - evaluation_points(i_grid))/dx
1319  b = (evaluation_points(i_grid) - x(lower_bound))/dx
1320  c = (a**3 - a)*const
1321  d = (b**3 - b)*const
1322 
1323  DO p_i = 1, nx
1324  y = 0
1325  y(p_i) = 1
1326  values(i_grid, p_i) = a*y(lower_bound) + b*y(upper_bound) &
1327  + (c*d2y_dx2(p_i, lower_bound) + d*d2y_dx2(p_i, upper_bound))
1328  END DO
1329  END DO
1330 !$OMP END DO
1331 
1332  DEALLOCATE (y)
1333 !$OMP END PARALLEL
1334  END SUBROUTINE spline_interpolation
1335 
1336 ! **************************************************************************************************
1337 !> \brief This routine is modeled after an algorithm from "Numerical Recipes in C" by Cambridge
1338 !> University Press, pages 96-97. It was adapted for Fortran and for the problem at hand.
1339 !> \param x ...
1340 !> \param d2y_dx2 ...
1341 !> \par History
1342 !> OpenMP added: Aug 2016 MTucker
1343 ! **************************************************************************************************
1344  SUBROUTINE initialize_spline_interpolation(x, d2y_dx2)
1345 
1346  REAL(dp), INTENT(in) :: x(:)
1347  REAL(dp), INTENT(inout) :: d2y_dx2(:, :)
1348 
1349  INTEGER :: index, nx, p_i
1350  REAL(dp) :: temp1, temp2
1351  REAL(dp), ALLOCATABLE :: temp_array(:), y(:)
1352 
1353  nx = SIZE(x)
1354 
1355 !$OMP PARALLEL DEFAULT( NONE ) &
1356 !$OMP SHARED( x, d2y_dx2, Nx ) &
1357 !$OMP PRIVATE( temp_array, y &
1358 !$OMP , index, temp1, temp2 &
1359 !$OMP )
1360 
1361  ALLOCATE (temp_array(nx), y(nx))
1362 
1363 !$OMP DO
1364  DO p_i = 1, nx
1365  !! In the Soler method, the polynomials that are interpolated are Kronecker delta functions
1366  !! at a particular q point. So, we set all y values to 0 except the one corresponding to
1367  !! the particular function P_i.
1368  !! ----------------------------------------------------------------------------------------
1369  y = 0.0_dp
1370  y(p_i) = 1.0_dp
1371  !! ----------------------------------------------------------------------------------------
1372 
1373  d2y_dx2(p_i, 1) = 0.0_dp
1374  temp_array(1) = 0.0_dp
1375  DO index = 2, nx - 1
1376  temp1 = (x(index) - x(index - 1))/(x(index + 1) - x(index - 1))
1377  temp2 = temp1*d2y_dx2(p_i, index - 1) + 2.0_dp
1378  d2y_dx2(p_i, index) = (temp1 - 1.0_dp)/temp2
1379  temp_array(index) = (y(index + 1) - y(index))/(x(index + 1) - x(index)) &
1380  - (y(index) - y(index - 1))/(x(index) - x(index - 1))
1381  temp_array(index) = (6.0_dp*temp_array(index)/(x(index + 1) - x(index - 1)) &
1382  - temp1*temp_array(index - 1))/temp2
1383  END DO
1384  d2y_dx2(p_i, nx) = 0.0_dp
1385  DO index = nx - 1, 1, -1
1386  d2y_dx2(p_i, index) = d2y_dx2(p_i, index)*d2y_dx2(p_i, index + 1) + temp_array(index)
1387  END DO
1388  END DO
1389 !$OMP END DO
1390 
1391  DEALLOCATE (temp_array, y)
1392 !$OMP END PARALLEL
1393 
1394  END SUBROUTINE initialize_spline_interpolation
1395 
1396  ! **************************************************************************************************
1397  !! This routine is modeled after an algorithm from "Numerical Recipes in C" by Cambridge
1398  !! University Press, page 97. Adapted for Fortran and the problem at hand. This function is used to
1399  !! find the Phi_alpha_beta needed for equations 8 and 11 of SOLER.
1400 ! **************************************************************************************************
1401 !> \brief ...
1402 !> \param k ...
1403 !> \param kernel_of_k ...
1404 !> \param dispersion_env ...
1405 !> \par History
1406 !> Optimised when adding OpenMP to vdw_energy (which calls this routine): Aug 2016 MTucker
1407 ! **************************************************************************************************
1408  SUBROUTINE interpolate_kernel(k, kernel_of_k, dispersion_env)
1409  REAL(dp), INTENT(in) :: k
1410  REAL(dp), INTENT(out) :: kernel_of_k(:, :)
1411  TYPE(qs_dispersion_type), POINTER :: dispersion_env
1412 
1413  INTEGER :: k_i, nr_points, q1_i, q2_i
1414  REAL(dp) :: a, b, c, const, d, dk, r_max
1415  REAL(dp), DIMENSION(:, :, :), POINTER :: d2phi_dk2, kernel
1416 
1417  nr_points = dispersion_env%nr_points
1418  r_max = dispersion_env%r_max
1419  dk = dispersion_env%dk
1420  kernel => dispersion_env%kernel
1421  d2phi_dk2 => dispersion_env%d2phi_dk2
1422 
1423  !! Check to make sure that the kernel table we have is capable of dealing with this
1424  !! value of k. If k is larger than Nr_points*2*pi/r_max then we can't perform the
1425  !! interpolation. In that case, a kernel file should be generated with a larger number
1426  !! of radial points.
1427  !! -------------------------------------------------------------------------------------
1428  cpassert(k < nr_points*dk)
1429  !! -------------------------------------------------------------------------------------
1430  kernel_of_k = 0.0_dp
1431  !! This integer division figures out which bin k is in since the kernel
1432  !! is set on a uniform grid.
1433  k_i = int(k/dk)
1434 
1435  !! Test to see if we are trying to interpolate a k that is one of the actual
1436  !! function points we have. The value is just the value of the function in that
1437  !! case.
1438  !! ----------------------------------------------------------------------------------------
1439  IF (mod(k, dk) == 0) THEN
1440  DO q1_i = 1, dispersion_env%Nqs
1441  DO q2_i = 1, q1_i
1442  kernel_of_k(q1_i, q2_i) = kernel(k_i, q1_i, q2_i)
1443  kernel_of_k(q2_i, q1_i) = kernel(k_i, q2_i, q1_i)
1444  END DO
1445  END DO
1446  RETURN
1447  END IF
1448  !! ----------------------------------------------------------------------------------------
1449  !! If we are not on a function point then we carry out the interpolation
1450  !! ----------------------------------------------------------------------------------------
1451  const = dk*dk/6.0_dp
1452  a = (dk*(k_i + 1.0_dp) - k)/dk
1453  b = (k - dk*k_i)/dk
1454  c = (a**3 - a)*const
1455  d = (b**3 - b)*const
1456  DO q1_i = 1, dispersion_env%Nqs
1457  DO q2_i = 1, q1_i
1458  kernel_of_k(q1_i, q2_i) = a*kernel(k_i, q1_i, q2_i) + b*kernel(k_i + 1, q1_i, q2_i) &
1459  + (c*d2phi_dk2(k_i, q1_i, q2_i) + d*d2phi_dk2(k_i + 1, q1_i, q2_i))
1460  kernel_of_k(q2_i, q1_i) = kernel_of_k(q1_i, q2_i)
1461  END DO
1462  END DO
1463 
1464  END SUBROUTINE interpolate_kernel
1465 
1466 ! **************************************************************************************************
1467 !> \brief ...
1468 !> \param k ...
1469 !> \param dkernel_of_dk ...
1470 !> \param dispersion_env ...
1471 !> \par History
1472 !> Optimised when adding OpenMP to vdw_energy (which calls this routine): Aug 2016 MTucker
1473 ! **************************************************************************************************
1474  SUBROUTINE interpolate_dkernel_dk(k, dkernel_of_dk, dispersion_env)
1475  REAL(dp), INTENT(in) :: k
1476  REAL(dp), INTENT(out) :: dkernel_of_dk(:, :)
1477  TYPE(qs_dispersion_type), POINTER :: dispersion_env
1478 
1479  INTEGER :: k_i, nr_points, q1_i, q2_i
1480  REAL(dp) :: a, b, dadk, dbdk, dcdk, dddk, dk, dk_6
1481  REAL(dp), DIMENSION(:, :, :), POINTER :: d2phi_dk2, kernel
1482 
1483  nr_points = dispersion_env%nr_points
1484  dk = dispersion_env%dk
1485  kernel => dispersion_env%kernel
1486  d2phi_dk2 => dispersion_env%d2phi_dk2
1487  dk_6 = dk/6.0_dp
1488 
1489  cpassert(k < nr_points*dk)
1490 
1491  dkernel_of_dk = 0.0_dp
1492  k_i = int(k/dk)
1493 
1494  a = (dk*(k_i + 1.0_dp) - k)/dk
1495  b = (k - dk*k_i)/dk
1496  dadk = -1.0_dp/dk
1497  dbdk = 1.0_dp/dk
1498  dcdk = -(3*a**2 - 1.0_dp)*dk_6
1499  dddk = (3*b**2 - 1.0_dp)*dk_6
1500  DO q1_i = 1, dispersion_env%Nqs
1501  DO q2_i = 1, q1_i
1502  dkernel_of_dk(q1_i, q2_i) = dadk*kernel(k_i, q1_i, q2_i) + dbdk*kernel(k_i + 1, q1_i, q2_i) &
1503  + dcdk*d2phi_dk2(k_i, q1_i, q2_i) + dddk*d2phi_dk2(k_i + 1, q1_i, q2_i)
1504  dkernel_of_dk(q2_i, q1_i) = dkernel_of_dk(q1_i, q2_i)
1505  END DO
1506  END DO
1507 
1508  END SUBROUTINE interpolate_dkernel_dk
1509 
1510 ! **************************************************************************************************
1511 
1512 END MODULE qs_dispersion_nonloc
bibliography
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition: bibliography.F:28
bibliography::dion2004
integer, save, public dion2004
Definition: bibliography.F:43
bibliography::romanperez2009
integer, save, public romanperez2009
Definition: bibliography.F:43
bibliography::sabatini2013
integer, save, public sabatini2013
Definition: bibliography.F:43
cp_files
Utility routines to open and close files. Tracking of preconnections.
Definition: cp_files.F:16
cp_files::open_file
subroutine, public open_file(file_name, file_status, file_form, file_action, file_position, file_pad, unit_number, debug, skip_get_unit_number, file_access)
Opens the requested file using a free unit number.
Definition: cp_files.F:308
cp_files::close_file
subroutine, public close_file(unit_number, file_status, keep_preconnection)
Close an open file given by its logical unit number. Optionally, keep the file and unit preconnected.
Definition: cp_files.F:119
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition: input_constants.F:17
input_constants::vdw_nl_rvv10
integer, parameter, public vdw_nl_rvv10
Definition: input_constants.F:598
input_constants::xc_vdw_fun_nonloc
integer, parameter, public xc_vdw_fun_nonloc
Definition: input_constants.F:592
input_constants::vdw_nl_drsll
integer, parameter, public vdw_nl_drsll
Definition: input_constants.F:598
input_constants::vdw_nl_lmkll
integer, parameter, public vdw_nl_lmkll
Definition: input_constants.F:598
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition: kinds.F:57
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:23
pw_grid_types
Definition: pw_grid_types.F:13
pw_grid_types::halfspace
integer, parameter, public halfspace
Definition: pw_grid_types.F:28
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:10106
pw_pool_types
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition: pw_pool_types.F:24
pw_types
Definition: pw_types.F:24
qs_dispersion_nonloc
Calculation of non local dispersion functionals Some routines adapted from: Copyright (C) 2001-2009 Q...
Definition: qs_dispersion_nonloc.F:19
qs_dispersion_nonloc::qs_dispersion_nonloc_init
subroutine, public qs_dispersion_nonloc_init(dispersion_env, para_env)
...
Definition: qs_dispersion_nonloc.F:67
qs_dispersion_nonloc::calculate_dispersion_nonloc
subroutine, public calculate_dispersion_nonloc(vxc_rho, rho_r, rho_g, edispersion, dispersion_env, energy_only, pw_pool, xc_pw_pool, para_env, virial)
Calculates the non-local vdW functional using the method of Soler For spin polarized cases we use E(a...
Definition: qs_dispersion_nonloc.F:165
qs_dispersion_types
Definition of disperson types for DFT calculations.
Definition: qs_dispersion_types.F:12
virial_types
Definition: virial_types.F:13