Line data Source code
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 : ! **************************************************************************************************
19 : MODULE qs_dispersion_nonloc
20 : USE bibliography, ONLY: Dion2004,&
21 : Romanperez2009,&
22 : Sabatini2013,&
23 : cite_reference
24 : USE cp_files, ONLY: close_file,&
25 : open_file
26 : USE input_constants, ONLY: vdw_nl_DRSLL,&
27 : vdw_nl_LMKLL,&
28 : vdw_nl_RVV10,&
29 : xc_vdw_fun_nonloc
30 : USE kinds, ONLY: default_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 :
60 : CONTAINS
61 :
62 : ! **************************************************************************************************
63 : !> \brief ...
64 : !> \param dispersion_env ...
65 : !> \param para_env ...
66 : ! **************************************************************************************************
67 46 : 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 46 : CALL timeset(routineN, handle)
78 :
79 46 : SELECT CASE (dispersion_env%nl_type)
80 : CASE DEFAULT
81 0 : CPABORT("Unknown vdW-DF functional")
82 : CASE (vdw_nl_DRSLL, vdw_nl_LMKLL)
83 32 : CALL cite_reference(Dion2004)
84 : CASE (vdw_nl_RVV10)
85 46 : CALL cite_reference(Sabatini2013)
86 : END SELECT
87 46 : CALL cite_reference(RomanPerez2009)
88 :
89 46 : vdw_type = dispersion_env%type
90 46 : SELECT CASE (vdw_type)
91 : CASE DEFAULT
92 : ! do nothing
93 : CASE (xc_vdw_fun_nonloc)
94 : ! setup information on non local functionals
95 46 : filename = dispersion_env%kernel_file_name
96 46 : IF (para_env%is_source()) THEN
97 : ! Read the kernel information from file "filename"
98 23 : CALL open_file(file_name=filename, unit_number=funit, file_form="FORMATTED")
99 23 : READ (funit, *) nqs, nr_points
100 23 : READ (funit, *) dispersion_env%r_max
101 : END IF
102 46 : CALL para_env%bcast(nqs)
103 46 : CALL para_env%bcast(nr_points)
104 46 : CALL para_env%bcast(dispersion_env%r_max)
105 : ALLOCATE (dispersion_env%q_mesh(nqs), dispersion_env%kernel(0:nr_points, nqs, nqs), &
106 460 : dispersion_env%d2phi_dk2(0:nr_points, nqs, nqs))
107 46 : dispersion_env%nqs = nqs
108 46 : dispersion_env%nr_points = nr_points
109 46 : IF (para_env%is_source()) THEN
110 : !! Read in the values of the q points used to generate this kernel
111 483 : 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 483 : DO q1_i = 1, nqs
116 5313 : DO q2_i = 1, q1_i
117 4830 : READ (funit, "(1p, 4e23.14)") dispersion_env%kernel(0:nr_points, q1_i, q2_i)
118 4956040 : 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 483 : DO q1_i = 1, nqs
126 5313 : DO q2_i = 1, q1_i
127 4830 : READ (funit, "(1p, 4e23.14)") dispersion_env%d2phi_dk2(0:nr_points, q1_i, q2_i)
128 4956040 : 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 23 : CALL close_file(unit_number=funit)
132 : END IF
133 1886 : CALL para_env%bcast(dispersion_env%q_mesh)
134 37758686 : CALL para_env%bcast(dispersion_env%kernel)
135 37758686 : CALL para_env%bcast(dispersion_env%d2phi_dk2)
136 : ! 2nd derivates for interpolation
137 184 : ALLOCATE (dispersion_env%d2y_dx2(nqs, nqs))
138 46 : CALL initialize_spline_interpolation(dispersion_env%q_mesh, dispersion_env%d2y_dx2)
139 : !
140 46 : dispersion_env%q_cut = dispersion_env%q_mesh(nqs)
141 46 : dispersion_env%q_min = dispersion_env%q_mesh(1)
142 92 : dispersion_env%dk = 2.0_dp*pi/dispersion_env%r_max
143 :
144 : END SELECT
145 :
146 46 : CALL timestop(handle)
147 :
148 46 : 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 388 : 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 388 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: dq0_dgradrho, dq0_drho, hpot, potential, &
184 388 : q0, rho
185 388 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: drho, thetas
186 : TYPE(pw_c1d_gs_type) :: tmp_g, vxc_g
187 388 : 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 388 : TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:, :) :: drho_r
191 :
192 388 : CALL timeset(routineN, handle)
193 :
194 388 : CPASSERT(ASSOCIATED(rho_r))
195 388 : CPASSERT(ASSOCIATED(rho_g))
196 388 : CPASSERT(ASSOCIATED(pw_pool))
197 :
198 388 : IF (PRESENT(virial)) THEN
199 382 : 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 108 : CPASSERT(.NOT. energy_only)
205 : END IF
206 388 : IF (.NOT. energy_only) THEN
207 382 : CPASSERT(ASSOCIATED(vxc_rho))
208 : END IF
209 :
210 388 : nl_type = dispersion_env%nl_type
211 :
212 388 : b_value = dispersion_env%b_value
213 388 : beta = 0.03125_dp*(3.0_dp/(b_value**2.0_dp))**0.75_dp
214 388 : nspin = SIZE(rho_r)
215 :
216 388 : const = 1.0_dp/(3.0_dp*rootpi*b_value**1.5_dp)/(pi**0.75_dp)
217 :
218 : ! temporary arrays for FFT
219 388 : CALL pw_pool%create_pw(tmp_g)
220 388 : CALL pw_pool%create_pw(tmp_r)
221 :
222 : ! get density derivatives
223 2796 : ALLOCATE (drho_r(3, nspin))
224 796 : DO ispin = 1, nspin
225 2020 : DO idir = 1, 3
226 1224 : CALL pw_pool%create_pw(drho_r(idir, ispin))
227 1224 : CALL pw_transfer(rho_g(ispin), tmp_g)
228 1224 : CALL pw_derive(tmp_g, nd(:, idir))
229 1632 : CALL pw_transfer(tmp_g, drho_r(idir, ispin))
230 : END DO
231 : END DO
232 :
233 1552 : np = SIZE(tmp_r%array)
234 1940 : ALLOCATE (rho(np), drho(np, 3)) !in the following loops, rho and drho _will_ have the same bounds
235 1552 : DO i = 1, 3
236 1164 : lo(i) = LBOUND(tmp_r%array, i)
237 1164 : hi(i) = UBOUND(tmp_r%array, i)
238 1552 : n(i) = hi(i) - lo(i) + 1
239 : END DO
240 : !$OMP PARALLEL DO DEFAULT(NONE) &
241 : !$OMP SHARED(n,rho) &
242 388 : !$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 1552 : DO i = 1, 3
252 : !$OMP PARALLEL DO DEFAULT(NONE) &
253 : !$OMP SHARED(n,i,drho) &
254 1552 : !$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 796 : DO ispin = 1, nspin
265 408 : CALL pw_transfer(rho_g(ispin), tmp_g)
266 408 : 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 408 : !$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 2020 : 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 1632 : !$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 388 : IF (energy_only) THEN
305 12 : ALLOCATE (q0(np))
306 0 : SELECT CASE (nl_type)
307 : CASE DEFAULT
308 0 : CPABORT("Unknown vdW-DF functional")
309 : CASE (vdw_nl_DRSLL, vdw_nl_LMKLL)
310 0 : CALL get_q0_on_grid_eo_vdw(rho, drho, q0, dispersion_env)
311 : CASE (vdw_nl_RVV10)
312 6 : CALL get_q0_on_grid_eo_rvv10(rho, drho, q0, dispersion_env)
313 : END SELECT
314 : ELSE
315 1528 : ALLOCATE (q0(np), dq0_drho(np), dq0_dgradrho(np))
316 0 : SELECT CASE (nl_type)
317 : CASE DEFAULT
318 0 : CPABORT("Unknown vdW-DF functional")
319 : CASE (vdw_nl_DRSLL, vdw_nl_LMKLL)
320 228 : CALL get_q0_on_grid_vdw(rho, drho, q0, dq0_drho, dq0_dgradrho, dispersion_env)
321 : CASE (vdw_nl_RVV10)
322 382 : 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 1552 : 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 388 : 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 0 : SELECT CASE (nl_type)
349 : CASE DEFAULT
350 0 : CPABORT("Unknown vdW-DF functional")
351 : CASE (vdw_nl_DRSLL, vdw_nl_LMKLL)
352 : !$OMP PARALLEL DO DEFAULT( NONE ) &
353 228 : !$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 388 : !$OMP SCHEDULE(DYNAMIC) ! use dynamic to allow for possibility of cases having (rho(i_grid) <= 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 1552 : DO i = 1, 3
377 1164 : lo(i) = LBOUND(tmp_r%array, i)
378 1164 : hi(i) = UBOUND(tmp_r%array, i)
379 1552 : n(i) = hi(i) - lo(i) + 1
380 : END DO
381 8924 : ALLOCATE (thetas_g(dispersion_env%nqs))
382 8148 : DO i = 1, dispersion_env%nqs
383 : !$OMP PARALLEL DO DEFAULT(NONE) &
384 : !$OMP SHARED(i,n,lo,thetas,tmp_r) &
385 7760 : !$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 7760 : CALL pw_pool%create_pw(thetas_g(i))
395 8148 : CALL pw_transfer(tmp_r, thetas_g(i))
396 : END DO
397 388 : 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 388 : sumnp = np
404 388 : CALL para_env%sum(sumnp)
405 388 : IF (use_virial) THEN
406 : ! calculates kernel contribution to stress
407 108 : CALL vdW_energy(thetas_g, dispersion_env, Ec_nl, energy_only, virial)
408 40 : SELECT CASE (nl_type)
409 : CASE (vdw_nl_RVV10)
410 276588 : 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 432 : DO idir = 1, 3
415 432 : virial%pv_xc(idir, idir) = virial%pv_xc(idir, idir) + Ec_nl
416 : END DO
417 : ELSE
418 280 : CALL vdW_energy(thetas_g, dispersion_env, Ec_nl, energy_only)
419 120 : SELECT CASE (nl_type)
420 : CASE (vdw_nl_RVV10)
421 565368 : Ec_nl = 0.5_dp*Ec_nl + beta*SUM(rho(:))*grid%vol/sumnp
422 : END SELECT
423 : END IF
424 388 : CALL para_env%sum(Ec_nl)
425 388 : IF (nl_type == vdw_nl_RVV10) Ec_nl = Ec_nl*dispersion_env%scale_rvv10
426 388 : edispersion = Ec_nl
427 :
428 388 : IF (energy_only) THEN
429 6 : 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 1528 : DO i = 1, 3
435 1146 : lo(i) = LBOUND(tmp_r%array, i)
436 1146 : hi(i) = UBOUND(tmp_r%array, i)
437 1528 : n(i) = hi(i) - lo(i) + 1
438 : END DO
439 8022 : DO i = 1, dispersion_env%nqs
440 7640 : CALL pw_transfer(thetas_g(i), tmp_r)
441 : !$OMP PARALLEL DO DEFAULT(NONE) &
442 : !$OMP SHARED(i,n,lo,thetas,tmp_r) &
443 8022 : !$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 1146 : ALLOCATE (potential(np), hpot(np))
461 382 : IF (use_virial) THEN
462 : ! calculates gradient contribution to stress
463 108 : grid => tmp_g%pw_grid
464 : CALL get_potential(q0, dq0_drho, dq0_dgradrho, rho, thetas, potential, hpot, &
465 108 : dispersion_env, drho, grid%dvol, virial)
466 : ELSE
467 : CALL get_potential(q0, dq0_drho, dq0_dgradrho, rho, thetas, potential, hpot, &
468 274 : dispersion_env)
469 : END IF
470 154 : SELECT CASE (nl_type)
471 : CASE (vdw_nl_RVV10)
472 743418 : potential(:) = (0.5_dp*potential(:) + beta)*dispersion_env%scale_rvv10
473 743800 : hpot(:) = 0.5_dp*dispersion_env%scale_rvv10*hpot(:)
474 : END SELECT
475 382 : CALL pw_pool%create_pw(vxc_r)
476 1528 : DO i = 1, 3
477 1146 : lo(i) = LBOUND(vxc_r%array, i)
478 1146 : hi(i) = UBOUND(vxc_r%array, i)
479 1528 : 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 382 : !$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 1528 : DO i = 1, 3
493 1146 : lo(i) = LBOUND(tmp_r%array, i)
494 1146 : hi(i) = UBOUND(tmp_r%array, i)
495 1528 : n(i) = hi(i) - lo(i) + 1
496 : END DO
497 1528 : DO idir = 1, 3
498 : !$OMP PARALLEL DO DEFAULT(NONE) &
499 : !$OMP SHARED(n,lo,tmp_r) &
500 1146 : !$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 2352 : DO ispin = 1, nspin
510 : !$OMP PARALLEL DO DEFAULT(NONE) &
511 : !$OMP SHARED(n,lo,tmp_r,hpot,drho_r,idir,ispin) &
512 2352 : !$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 1146 : CALL pw_transfer(tmp_r, tmp_g)
525 1146 : CALL pw_derive(tmp_g, nd(:, idir))
526 1146 : CALL pw_transfer(tmp_g, tmp_r)
527 1528 : CALL pw_axpy(tmp_r, vxc_r, -1._dp)
528 : END DO
529 382 : CALL pw_transfer(vxc_r, tmp_g)
530 382 : CALL pw_pool%give_back_pw(vxc_r)
531 382 : CALL xc_pw_pool%create_pw(vxc_r)
532 382 : CALL xc_pw_pool%create_pw(vxc_g)
533 382 : CALL pw_transfer(tmp_g, vxc_g)
534 382 : CALL pw_transfer(vxc_g, vxc_r)
535 784 : DO ispin = 1, nspin
536 784 : CALL pw_axpy(vxc_r, vxc_rho(ispin), 1._dp)
537 : END DO
538 382 : CALL xc_pw_pool%give_back_pw(vxc_r)
539 382 : CALL xc_pw_pool%give_back_pw(vxc_g)
540 382 : DEALLOCATE (q0, dq0_drho, dq0_dgradrho)
541 : END IF
542 :
543 388 : DEALLOCATE (thetas)
544 :
545 8148 : DO i = 1, dispersion_env%nqs
546 8148 : CALL pw_pool%give_back_pw(thetas_g(i))
547 : END DO
548 796 : DO ispin = 1, nspin
549 2020 : DO idir = 1, 3
550 1632 : CALL pw_pool%give_back_pw(drho_r(idir, ispin))
551 : END DO
552 : END DO
553 388 : CALL pw_pool%give_back_pw(tmp_r)
554 388 : CALL pw_pool%give_back_pw(tmp_g)
555 :
556 388 : DEALLOCATE (rho, drho, drho_r, thetas_g)
557 :
558 388 : CALL timestop(handle)
559 :
560 776 : 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 388 : 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 388 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: theta
586 388 : 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 388 : 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 388 : CALL timeset(routineN, handle)
596 388 : nqs = dispersion_env%nqs
597 :
598 388 : use_virial = PRESENT(virial)
599 388 : virial_thread(:, :) = 0.0_dp ! always initialize to avoid floating point exceptions in OMP REDUCTION
600 :
601 388 : vdW_xc_energy = 0._dp
602 388 : grid => thetas_g(1)%pw_grid
603 :
604 388 : IF (grid%grid_span == HALFSPACE) THEN
605 : g_multiplier = 2._dp
606 : ELSE
607 388 : g_multiplier = 1._dp
608 : END IF
609 :
610 388 : nl_type = dispersion_env%nl_type
611 :
612 388 : IF (.NOT. energy_only) THEN
613 1528 : ALLOCATE (u_vdW(grid%ngpts_cut_local, nqs))
614 86161402 : 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 388 : !$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 388 : IF (.NOT. energy_only) THEN
692 382 : DEALLOCATE (u_vdW)
693 : END IF
694 :
695 388 : vdW_xc_energy = vdW_xc_energy*grid%vol*0.5_dp
696 :
697 388 : IF (use_virial) THEN
698 432 : DO l = 1, 3
699 648 : DO m = 1, (l - 1)
700 324 : virial%pv_xc(l, m) = virial%pv_xc(l, m) + virial_thread(l, m)
701 648 : virial%pv_xc(m, l) = virial%pv_xc(l, m)
702 : END DO
703 324 : m = l
704 432 : virial%pv_xc(l, m) = virial%pv_xc(l, m) + virial_thread(l, m)
705 : END DO
706 : END IF
707 :
708 388 : CALL timestop(handle)
709 :
710 388 : 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 382 : SUBROUTINE get_potential(q0, dq0_drho, dq0_dgradrho, total_rho, u_vdW, potential, h_prefactor, &
736 382 : 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 382 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: y
755 : REAL(dp), DIMENSION(3, 3) :: virial_thread
756 382 : REAL(dp), DIMENSION(:), POINTER :: q_mesh
757 382 : REAL(dp), DIMENSION(:, :), POINTER :: d2y_dx2
758 :
759 382 : CALL timeset(routineN, handle)
760 :
761 382 : use_virial = PRESENT(virial)
762 382 : CPASSERT(.NOT. use_virial .OR. PRESENT(drho))
763 382 : CPASSERT(.NOT. use_virial .OR. PRESENT(dvol))
764 :
765 382 : virial_thread(:, :) = 0.0_dp ! always initialize to avoid floating point exceptions in OMP REDUCTION
766 382 : b_value = dispersion_env%b_value
767 382 : const = 1.0_dp/(3.0_dp*b_value**(3.0_dp/2.0_dp)*pi**(5.0_dp/4.0_dp))
768 4308051 : potential = 0.0_dp
769 4308051 : h_prefactor = 0.0_dp
770 :
771 382 : d2y_dx2 => dispersion_env%d2y_dx2
772 382 : q_mesh => dispersion_env%q_mesh
773 382 : nqs = dispersion_env%nqs
774 382 : 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 382 : !$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) /= 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 382 : IF (use_virial) THEN
869 432 : DO l = 1, 3
870 648 : DO m = 1, (l - 1)
871 324 : virial%pv_xc(l, m) = virial%pv_xc(l, m) + virial_thread(l, m)
872 648 : virial%pv_xc(m, l) = virial%pv_xc(l, m)
873 : END DO
874 324 : m = l
875 432 : virial%pv_xc(l, m) = virial%pv_xc(l, m) + virial_thread(l, m)
876 : END DO
877 : END IF
878 :
879 382 : CALL timestop(handle)
880 764 : 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 72326 : 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 72326 : multiplier = alpha
899 72326 : exponent = alpha
900 :
901 867912 : DO i = 2, hi
902 795586 : multiplier = multiplier*alpha
903 867912 : exponent = exponent + (multiplier/i)
904 : END DO
905 72326 : 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 4145733 : 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 4145733 : derivative = 0.0d0
926 4145733 : multiplier = 1.0d0
927 4145733 : exponent = 0.0d0
928 :
929 53894529 : DO i = 1, hi
930 49748796 : derivative = derivative + multiplier
931 49748796 : multiplier = multiplier*alpha
932 53894529 : exponent = exponent + (multiplier/i)
933 : END DO
934 4145733 : 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 228 : 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 228 : q_cut = dispersion_env%q_cut
970 228 : q_min = dispersion_env%q_min
971 228 : SELECT CASE (dispersion_env%nl_type)
972 : CASE DEFAULT
973 0 : CPABORT("Unknown vdW-DF functional")
974 : CASE (vdw_nl_DRSLL)
975 48 : Z_ab = -0.8491_dp
976 : CASE (vdw_nl_LMKLL)
977 228 : Z_ab = -1.887_dp
978 : END SELECT
979 :
980 : ! initialize q0-related arrays ...
981 3564633 : q0(:) = q_cut
982 3564633 : dq0_drho(:) = 0.0_dp
983 3564633 : dq0_dgradrho(:) = 0.0_dp
984 :
985 3564633 : 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 3564405 : IF (total_rho(i_grid) < epsr) CYCLE
994 : !! ------------------------------------------------------------------------------------
995 : !! Calculate some intermediate values needed to find q
996 : !! ------------------------------------------------------------------------------------
997 3492040 : kF = (3.0_dp*pi*pi*total_rho(i_grid))**(1.0_dp/3.0_dp)
998 3492040 : r_s = (3.0_dp/(4.0_dp*pi*total_rho(i_grid)))**(1.0_dp/3.0_dp)
999 3492040 : sqrt_r_s = SQRT(r_s)
1000 :
1001 : gradient_correction = -Z_ab/(36.0_dp*kF*total_rho(i_grid)**2) &
1002 3492040 : *(gradient_rho(i_grid, 1)**2 + gradient_rho(i_grid, 2)**2 + gradient_rho(i_grid, 3)**2)
1003 :
1004 3492040 : LDA_1 = 8.0_dp*pi/3.0_dp*(LDA_A*(1.0_dp + LDA_a1*r_s))
1005 3492040 : 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 3492040 : 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 3492040 : q__q_cut = q/q_cut
1015 3492040 : CALL calculate_exponent_derivative(m_cut, q__q_cut, exponent, dq0_dq)
1016 3492040 : q0(i_grid) = q_cut*(1.0_dp - EXP(-exponent))
1017 3492040 : 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 3492040 : IF (q0(i_grid) < q_min) THEN
1023 0 : 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 3492040 : + 2.0_dp*LDA_b4/3.0_dp*r_s**2)))
1042 :
1043 3564633 : 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 228 : 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 154 : 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 154 : q_cut = dispersion_env%q_cut
1078 154 : q_min = dispersion_env%q_min
1079 154 : b_value = dispersion_env%b_value
1080 154 : C_value = dispersion_env%c_value
1081 :
1082 : ! initialize q0-related arrays ...
1083 743418 : q0(:) = q_cut
1084 743418 : dq0_drho(:) = 0.0_dp
1085 743418 : dq0_dgradrho(:) = 0.0_dp
1086 :
1087 743418 : DO i_grid = 1, SIZE(total_rho)
1088 :
1089 743264 : 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 743418 : IF (total_rho(i_grid) > epsr) THEN
1093 :
1094 : !! Calculate some intermediate values needed to find q
1095 : !! ------------------------------------------------------------------------------------
1096 653693 : mod_grad = SQRT(gmod2)
1097 :
1098 653693 : wp2 = 16.0_dp*pi*total_rho(i_grid)
1099 653693 : wg2 = 4_dp*C_value*(mod_grad/total_rho(i_grid))**4
1100 :
1101 653693 : k = b_value*3.0_dp*pi*((total_rho(i_grid)/(9.0_dp*pi))**(1.0_dp/6.0_dp))
1102 653693 : w0 = SQRT(wg2 + wp2/3.0_dp)
1103 :
1104 653693 : q = w0/k
1105 :
1106 : !! Here, we calculate q0 by saturating q according
1107 : !! ---------------------------------------------------------------------------------------
1108 653693 : q__q_cut = q/q_cut
1109 653693 : CALL calculate_exponent_derivative(m_cut, q__q_cut, exponent, dq0_dq)
1110 653693 : q0(i_grid) = q_cut*(1.0_dp - EXP(-exponent))
1111 653693 : dq0_dq = dq0_dq*EXP(-exponent)
1112 :
1113 : !! ---------------------------------------------------------------------------------------
1114 653693 : IF (q0(i_grid) < q_min) THEN
1115 3 : q0(i_grid) = q_min
1116 : END IF
1117 :
1118 : !!---------------------------------Final values---------------------------------
1119 653693 : 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 653693 : dk_dn = k/(6.0_dp*total_rho(i_grid))
1121 :
1122 653693 : dq0_drho(i_grid) = dq0_dq*1.0_dp/(k**2.0)*(dw0_dn*k - dk_dn*w0)
1123 653693 : 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 154 : 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 0 : 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 0 : q_cut = dispersion_env%q_cut
1153 0 : q_min = dispersion_env%q_min
1154 0 : SELECT CASE (dispersion_env%nl_type)
1155 : CASE DEFAULT
1156 0 : CPABORT("Unknown vdW-DF functional")
1157 : CASE (vdw_nl_DRSLL)
1158 0 : Z_ab = -0.8491_dp
1159 : CASE (vdw_nl_LMKLL)
1160 0 : Z_ab = -1.887_dp
1161 : END SELECT
1162 :
1163 : ! initialize q0-related arrays ...
1164 0 : q0(:) = q_cut
1165 0 : 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 0 : IF (total_rho(i_grid) < epsr) CYCLE
1173 : !! ------------------------------------------------------------------------------------
1174 : !! Calculate some intermediate values needed to find q
1175 : !! ------------------------------------------------------------------------------------
1176 0 : kF = (3.0_dp*pi*pi*total_rho(i_grid))**(1.0_dp/3.0_dp)
1177 0 : r_s = (3.0_dp/(4.0_dp*pi*total_rho(i_grid)))**(1.0_dp/3.0_dp)
1178 0 : sqrt_r_s = SQRT(r_s)
1179 :
1180 : gradient_correction = -Z_ab/(36.0_dp*kF*total_rho(i_grid)**2) &
1181 0 : *(gradient_rho(i_grid, 1)**2 + gradient_rho(i_grid, 2)**2 + gradient_rho(i_grid, 3)**2)
1182 :
1183 0 : LDA_1 = 8.0_dp*pi/3.0_dp*(LDA_A*(1.0_dp + LDA_a1*r_s))
1184 0 : 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 0 : 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 0 : q__q_cut = q/q_cut
1195 0 : CALL calculate_exponent(m_cut, q__q_cut, exponent)
1196 0 : 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 0 : IF (q0(i_grid) < q_min) THEN
1203 0 : q0(i_grid) = q_min
1204 : END IF
1205 : END DO
1206 :
1207 0 : 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 6 : 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 6 : q_cut = dispersion_env%q_cut
1229 6 : q_min = dispersion_env%q_min
1230 6 : b_value = dispersion_env%b_value
1231 6 : C_value = dispersion_env%c_value
1232 :
1233 : ! initialize q0-related arrays ...
1234 98310 : q0(:) = q_cut
1235 :
1236 98310 : DO i_grid = 1, SIZE(total_rho)
1237 :
1238 98304 : 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 98310 : IF (total_rho(i_grid) > epsr) THEN
1242 :
1243 : !! Calculate some intermediate values needed to find q
1244 : !! ------------------------------------------------------------------------------------
1245 72326 : wp2 = 16.0_dp*pi*total_rho(i_grid)
1246 72326 : wg2 = 4_dp*C_value*(gmod2*gmod2)/(total_rho(i_grid)**4)
1247 :
1248 72326 : k = b_value*3.0_dp*pi*((total_rho(i_grid)/(9.0_dp*pi))**(1.0_dp/6.0_dp))
1249 72326 : w0 = SQRT(wg2 + wp2/3.0_dp)
1250 :
1251 72326 : q = w0/k
1252 :
1253 : !! Here, we calculate q0 by saturating q according
1254 : !! ---------------------------------------------------------------------------------------
1255 72326 : q__q_cut = q/q_cut
1256 72326 : CALL calculate_exponent(m_cut, q__q_cut, exponent)
1257 72326 : q0(i_grid) = q_cut*(1.0_dp - EXP(-exponent))
1258 :
1259 72326 : IF (q0(i_grid) < q_min) THEN
1260 1 : q0(i_grid) = q_min
1261 : END IF
1262 :
1263 : END IF
1264 :
1265 : END DO
1266 :
1267 6 : 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 388 : 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 388 : REAL(dp), ALLOCATABLE :: y(:)
1289 :
1290 388 : Nx = SIZE(x)
1291 388 : 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 388 : !$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 388 : 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 46 : 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 46 : REAL(dp), ALLOCATABLE :: temp_array(:), y(:)
1353 :
1354 46 : 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 46 : !$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 46 : 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 386626 : 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 386626 : REAL(dp), DIMENSION(:, :, :), POINTER :: d2phi_dk2, kernel
1417 :
1418 386626 : Nr_points = dispersion_env%nr_points
1419 386626 : r_max = dispersion_env%r_max
1420 386626 : dk = dispersion_env%dk
1421 386626 : kernel => dispersion_env%kernel
1422 386626 : 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 0 : CPASSERT(k < Nr_points*dk)
1430 : !! -------------------------------------------------------------------------------------
1431 162769546 : 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 386626 : 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 386626 : IF (MOD(k, dk) == 0) THEN
1441 4074 : DO q1_i = 1, dispersion_env%Nqs
1442 44814 : DO q2_i = 1, q1_i
1443 40740 : kernel_of_k(q1_i, q2_i) = kernel(k_i, q1_i, q2_i)
1444 44620 : kernel_of_k(q2_i, q1_i) = kernel(k_i, q2_i, q1_i)
1445 : END DO
1446 : END DO
1447 386626 : RETURN
1448 : END IF
1449 : !! ----------------------------------------------------------------------------------------
1450 : !! If we are not on a function point then we carry out the interpolation
1451 : !! ----------------------------------------------------------------------------------------
1452 386432 : const = dk*dk/6.0_dp
1453 386432 : A = (dk*(k_i + 1.0_dp) - k)/dk
1454 386432 : B = (k - dk*k_i)/dk
1455 386432 : C = (A**3 - A)*const
1456 386432 : D = (B**3 - B)*const
1457 8115072 : DO q1_i = 1, dispersion_env%Nqs
1458 89265792 : 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 81150720 : + (C*d2phi_dk2(k_i, q1_i, q2_i) + D*d2phi_dk2(k_i + 1, q1_i, q2_i))
1461 88879360 : kernel_of_k(q2_i, q1_i) = kernel_of_k(q1_i, q2_i)
1462 : END DO
1463 : END DO
1464 :
1465 386626 : 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 297321 : 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 297321 : REAL(dp), DIMENSION(:, :, :), POINTER :: d2phi_dk2, kernel
1483 :
1484 297321 : Nr_points = dispersion_env%nr_points
1485 297321 : dk = dispersion_env%dk
1486 297321 : kernel => dispersion_env%kernel
1487 297321 : d2phi_dk2 => dispersion_env%d2phi_dk2
1488 297321 : dk_6 = dk/6.0_dp
1489 :
1490 0 : CPASSERT(k < Nr_points*dk)
1491 :
1492 125172141 : dkernel_of_dk = 0.0_dp
1493 297321 : k_i = INT(k/dk)
1494 :
1495 297321 : A = (dk*(k_i + 1.0_dp) - k)/dk
1496 297321 : B = (k - dk*k_i)/dk
1497 297321 : dAdk = -1.0_dp/dk
1498 297321 : dBdk = 1.0_dp/dk
1499 297321 : dCdk = -(3*A**2 - 1.0_dp)*dk_6
1500 297321 : dDdk = (3*B**2 - 1.0_dp)*dk_6
1501 6243741 : DO q1_i = 1, dispersion_env%Nqs
1502 68681151 : 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 62437410 : + dCdk*d2phi_dk2(k_i, q1_i, q2_i) + dDdk*d2phi_dk2(k_i + 1, q1_i, q2_i)
1505 68383830 : dkernel_of_dk(q2_i, q1_i) = dkernel_of_dk(q1_i, q2_i)
1506 : END DO
1507 : END DO
1508 :
1509 297321 : END SUBROUTINE interpolate_dkernel_dk
1510 :
1511 : ! **************************************************************************************************
1512 :
1513 : END MODULE qs_dispersion_nonloc
|
