(git:374b731)
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-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! **************************************************************************************************
19MODULE qs_dispersion_nonloc
20 USE bibliography, ONLY: dion2004,&
21 romanperez2009,&
22 sabatini2013,&
23 cite_reference
24 USE cp_files, ONLY: close_file,&
25 open_file
26 USE input_constants, ONLY: vdw_nl_drsll,&
27 vdw_nl_lmkll,&
28 vdw_nl_rvv10,&
29 xc_vdw_fun_nonloc
30 USE kinds, ONLY: default_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
59CONTAINS
60
61! **************************************************************************************************
62!> \brief ...
63!> \param dispersion_env ...
64!> \param para_env ...
65! **************************************************************************************************
66 SUBROUTINE qs_dispersion_nonloc_init(dispersion_env, para_env)
67 TYPE(qs_dispersion_type), POINTER :: dispersion_env
68 TYPE(mp_para_env_type), POINTER :: para_env
69
70 CHARACTER(len=*), PARAMETER :: routinen = 'qs_dispersion_nonloc_init'
71
72 CHARACTER(LEN=default_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
1512END MODULE qs_dispersion_nonloc
pw_methods::pw_axpy
Definition pw_methods.F:164
pw_methods::pw_transfer
Definition pw_methods.F:303
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
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
pw_grid_types::pw_grid_type
Definition pw_grid_types.F:62
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:31
virial_types::virial_type
Definition virial_types.F:26