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