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