Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Exchange and Correlation functional calculations
10 : !> \par History
11 : !> (13-Feb-2001) JGH, based on earlier version of apsi
12 : !> 02.2003 Many many changes [fawzi]
13 : !> 03.2004 new xc interface [fawzi]
14 : !> 04.2004 kinetic functionals [fawzi]
15 : !> \author fawzi
16 : ! **************************************************************************************************
17 : MODULE xc
18 : #:include 'xc.fypp'
19 : USE cp_array_utils, ONLY: cp_3d_r_cp_type
20 : USE cp_linked_list_xc_deriv, ONLY: cp_sll_xc_deriv_next, &
21 : cp_sll_xc_deriv_type
22 : USE cp_log_handling, ONLY: cp_get_default_logger, &
23 : cp_logger_get_default_unit_nr, &
24 : cp_logger_type, &
25 : cp_to_string
26 : USE input_section_types, ONLY: section_get_ival, &
27 : section_get_lval, &
28 : section_get_rval, &
29 : section_vals_get_subs_vals, &
30 : section_vals_type, &
31 : section_vals_val_get
32 : USE kahan_sum, ONLY: accurate_dot_product, &
33 : accurate_sum
34 : USE kinds, ONLY: default_path_length, &
35 : dp
36 : USE pw_grid_types, ONLY: PW_MODE_DISTRIBUTED, &
37 : pw_grid_type
38 : USE pw_methods, ONLY: pw_axpy, &
39 : pw_copy, &
40 : pw_derive, &
41 : pw_scale, &
42 : pw_transfer, &
43 : pw_zero, pw_integrate_function
44 : USE pw_pool_types, ONLY: &
45 : pw_pool_type
46 : USE pw_types, ONLY: &
47 : pw_c1d_gs_type, pw_r3d_rs_type
48 : USE xc_derivative_desc, ONLY: &
49 : deriv_rho, deriv_rhoa, deriv_rhob, &
50 : deriv_norm_drhoa, deriv_norm_drhob, deriv_norm_drho, deriv_tau_a, deriv_tau_b, deriv_tau, &
51 : deriv_laplace_rho, deriv_laplace_rhoa, deriv_laplace_rhob, id_to_desc
52 : USE xc_derivative_set_types, ONLY: xc_derivative_set_type, &
53 : xc_dset_create, &
54 : xc_dset_get_derivative, &
55 : xc_dset_release, &
56 : xc_dset_zero_all, xc_dset_recover_pw
57 : USE xc_derivative_types, ONLY: xc_derivative_get, &
58 : xc_derivative_type
59 : USE xc_derivatives, ONLY: xc_functionals_eval, &
60 : xc_functionals_get_needs
61 : USE xc_input_constants, ONLY: &
62 : xc_debug_new_routine, xc_default_f_routine, xc_test_lsd_f_routine
63 : USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
64 : USE xc_rho_set_types, ONLY: xc_rho_set_create, &
65 : xc_rho_set_get, &
66 : xc_rho_set_release, &
67 : xc_rho_set_type, &
68 : xc_rho_set_update, xc_rho_set_recover_pw
69 : USE xc_util, ONLY: xc_pw_smooth, xc_pw_laplace, xc_pw_divergence, xc_requires_tmp_g
70 : #include "../base/base_uses.f90"
71 :
72 : IMPLICIT NONE
73 : PRIVATE
74 : PUBLIC :: xc_vxc_pw_create1, xc_vxc_pw_create, &
75 : xc_exc_calc, xc_calc_2nd_deriv_analytical, xc_calc_2nd_deriv_numerical, &
76 : xc_calc_2nd_deriv, xc_prep_2nd_deriv, divide_by_norm_drho, smooth_cutoff, &
77 : xc_uses_kinetic_energy_density, xc_uses_norm_drho
78 : PUBLIC :: calc_xc_density
79 :
80 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
81 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc'
82 :
83 : CONTAINS
84 :
85 : ! **************************************************************************************************
86 : !> \brief ...
87 : !> \param xc_fun_section ...
88 : !> \param lsd ...
89 : !> \return ...
90 : ! **************************************************************************************************
91 6906 : FUNCTION xc_uses_kinetic_energy_density(xc_fun_section, lsd) RESULT(res)
92 : TYPE(section_vals_type), POINTER, INTENT(IN) :: xc_fun_section
93 : LOGICAL, INTENT(IN) :: lsd
94 : LOGICAL :: res
95 :
96 : TYPE(xc_rho_cflags_type) :: needs
97 :
98 : needs = xc_functionals_get_needs(xc_fun_section, &
99 : lsd=lsd, &
100 13812 : calc_potential=.FALSE.)
101 6906 : res = (needs%tau_spin .OR. needs%tau)
102 :
103 6906 : END FUNCTION
104 :
105 : ! **************************************************************************************************
106 : !> \brief ...
107 : !> \param xc_fun_section ...
108 : !> \param lsd ...
109 : !> \return ...
110 : ! **************************************************************************************************
111 6922 : FUNCTION xc_uses_norm_drho(xc_fun_section, lsd) RESULT(res)
112 : TYPE(section_vals_type), POINTER, INTENT(IN) :: xc_fun_section
113 : LOGICAL, INTENT(IN) :: lsd
114 : LOGICAL :: res
115 :
116 : TYPE(xc_rho_cflags_type) :: needs
117 :
118 : needs = xc_functionals_get_needs(xc_fun_section, &
119 : lsd=lsd, &
120 13844 : calc_potential=.FALSE.)
121 6922 : res = (needs%norm_drho .OR. needs%norm_drho_spin)
122 :
123 6922 : END FUNCTION
124 :
125 : ! **************************************************************************************************
126 : !> \brief Exchange and Correlation functional calculations.
127 : !> depending on the selected functional_routine calls
128 : !> the correct routine
129 : !> \param vxc_rho will contain the v_xc part that depend on rho
130 : !> (if one of the chosen xc functionals has it it is allocated and you
131 : !> are responsible for it)
132 : !> \param vxc_tau will contain the kinetic tau part of v_xcthe functional
133 : !> (if one of the chosen xc functionals has it it is allocated and you
134 : !> are responsible for it)
135 : !> \param rho_r the value of the density in the real space
136 : !> \param rho_g value of the density in the g space (needs to be associated
137 : !> only for gradient corrections)
138 : !> \param tau value of the kinetic density tau on the grid (can be null,
139 : !> used only with meta functionals)
140 : !> \param exc the xc energy
141 : !> \param xc_section parameters selecting the xc and the method used to
142 : !> calculate it
143 : !> \param pw_pool the pool for the grids
144 : !> \param compute_virial should the virial be computed... if so virial_xc must be present
145 : !> \param virial_xc for calculating the GGA part of the stress
146 : !> \param exc_r ...
147 : !> \author fawzi
148 : ! **************************************************************************************************
149 210592 : SUBROUTINE xc_vxc_pw_create1(vxc_rho, vxc_tau, rho_r, rho_g, tau, exc, xc_section, &
150 : pw_pool, compute_virial, virial_xc, exc_r)
151 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: vxc_rho, vxc_tau
152 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, tau
153 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
154 : REAL(KIND=dp), INTENT(out) :: exc
155 : TYPE(section_vals_type), POINTER :: xc_section
156 : TYPE(pw_pool_type), POINTER :: pw_pool
157 : LOGICAL, INTENT(IN) :: compute_virial
158 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT) :: virial_xc
159 : TYPE(pw_r3d_rs_type), INTENT(INOUT), OPTIONAL :: exc_r
160 :
161 : INTEGER :: f_routine
162 :
163 105296 : CPASSERT(ASSOCIATED(rho_r))
164 105296 : CPASSERT(ASSOCIATED(xc_section))
165 105296 : CPASSERT(.NOT. ASSOCIATED(vxc_rho))
166 105296 : CPASSERT(.NOT. ASSOCIATED(vxc_tau))
167 :
168 105296 : virial_xc = 0.0_dp
169 :
170 : CALL section_vals_val_get(xc_section, "FUNCTIONAL_ROUTINE", &
171 105296 : i_val=f_routine)
172 105296 : SELECT CASE (f_routine)
173 : CASE (xc_default_f_routine)
174 : CALL xc_vxc_pw_create(vxc_rho=vxc_rho, vxc_tau=vxc_tau, tau=tau, exc_r=exc_r, &
175 : rho_r=rho_r, rho_g=rho_g, exc=exc, xc_section=xc_section, &
176 105296 : pw_pool=pw_pool, compute_virial=compute_virial, virial_xc=virial_xc)
177 : CASE (xc_debug_new_routine)
178 0 : CPASSERT(.NOT. PRESENT(exc_r))
179 : CALL xc_vxc_pw_create_debug(vxc_rho=vxc_rho, vxc_tau=vxc_tau, tau=tau, &
180 : rho_r=rho_r, rho_g=rho_g, exc=exc, xc_section=xc_section, &
181 0 : pw_pool=pw_pool)
182 : CASE (xc_test_lsd_f_routine)
183 0 : CPASSERT(.NOT. PRESENT(exc_r))
184 : CALL xc_vxc_pw_create_test_lsd(vxc_rho=vxc_rho, vxc_tau=vxc_tau, &
185 : tau=tau, rho_r=rho_r, rho_g=rho_g, exc=exc, &
186 105296 : xc_section=xc_section, pw_pool=pw_pool)
187 : CASE default
188 : END SELECT
189 :
190 105296 : END SUBROUTINE xc_vxc_pw_create1
191 :
192 : ! **************************************************************************************************
193 : !> \brief calculates vxc using lsd with rhoa=rhob=0.5*rho and compares
194 : !> with the lda result
195 : !> \param vxc_rho will contain the v_xc part that depend on rho
196 : !> (if one of the chosen xc functionals has it it is allocated and you
197 : !> are responsible for it)
198 : !> \param vxc_tau will contain the kinetic tau part of v_xc
199 : !> (if one of the chosen xc functionals has it it is allocated and you
200 : !> are responsible for it)
201 : !> \param rho_r the value of the density in the real space
202 : !> \param rho_g value of the density in the g space (needs to be associated
203 : !> only for gradient corrections)
204 : !> \param tau value of the kinetic density tau on the grid (can be null,
205 : !> used only with meta functionals)
206 : !> \param exc the xc energy
207 : !> \param xc_section which functional to calculate, and how
208 : !> \param pw_pool the pool for the grids
209 : !> \author Fawzi Mohamed
210 : !> \note
211 : !> for debugging only: leaks, and non parallel
212 : ! **************************************************************************************************
213 0 : SUBROUTINE xc_vxc_pw_create_test_lsd(vxc_rho, vxc_tau, rho_r, rho_g, tau, &
214 : exc, xc_section, pw_pool)
215 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: vxc_rho, vxc_tau
216 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, tau
217 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
218 : REAL(KIND=dp), INTENT(out) :: exc
219 : TYPE(section_vals_type), POINTER :: xc_section
220 : TYPE(pw_pool_type), POINTER :: pw_pool
221 :
222 : LOGICAL, PARAMETER :: compute_virial = .FALSE.
223 :
224 : CHARACTER(len=default_path_length) :: filename
225 0 : INTEGER, DIMENSION(:), POINTER :: split_desc
226 : INTEGER :: i, ii, ispin, j, k, order
227 : INTEGER, DIMENSION(2, 3) :: bo
228 : REAL(kind=dp) :: diff, exc2, maxdiff, tmp
229 : REAL(KIND=dp), DIMENSION(3, 3) :: virial_xc
230 0 : REAL(kind=dp), DIMENSION(:, :, :), POINTER :: pot, pot2, pot3
231 : TYPE(cp_sll_xc_deriv_type), POINTER :: deriv_iter
232 0 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: vxc_rho2, vxc_tau2
233 0 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho2_g
234 0 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho2_r, tau2
235 : TYPE(xc_derivative_set_type) :: dSet1, dSet2
236 : TYPE(xc_derivative_type), POINTER :: deriv, deriv2, deriv3
237 : TYPE(xc_rho_set_type) :: rho_set1, rho_set2
238 :
239 0 : NULLIFY (vxc_rho2, vxc_tau2, tau2, pot, pot3, pot3, &
240 0 : deriv, deriv2, deriv3, rho2_g)
241 :
242 0 : bo = rho_r(1)%pw_grid%bounds_local
243 :
244 0 : ALLOCATE (rho2_r(2))
245 0 : DO ispin = 1, 2
246 0 : CALL pw_pool%create_pw(rho2_r(ispin))
247 : END DO
248 0 : DO k = bo(1, 3), bo(2, 3)
249 0 : DO j = bo(1, 2), bo(2, 2)
250 0 : DO i = bo(1, 1), bo(2, 1)
251 0 : tmp = rho_r(1)%array(i, j, k)*0.5
252 0 : rho2_r(1)%array(i, j, k) = tmp
253 0 : rho2_r(2)%array(i, j, k) = tmp
254 : END DO
255 : END DO
256 : END DO
257 :
258 0 : IF (ASSOCIATED(tau)) THEN
259 0 : ALLOCATE (tau2(2))
260 0 : DO ispin = 1, 2
261 0 : CALL pw_pool%create_pw(tau2(ispin))
262 : END DO
263 :
264 0 : DO k = bo(1, 3), bo(2, 3)
265 0 : DO j = bo(1, 2), bo(2, 2)
266 0 : DO i = bo(1, 1), bo(2, 1)
267 0 : tmp = tau(1)%array(i, j, k)*0.5
268 0 : tau2(1)%array(i, j, k) = tmp
269 0 : tau2(2)%array(i, j, k) = tmp
270 : END DO
271 : END DO
272 : END DO
273 : END IF
274 :
275 0 : PRINT *, "about to calculate xc (lda)"
276 : CALL xc_rho_set_and_dset_create(rho_r=rho_r, rho_g=rho_g, &
277 : tau=tau, xc_section=xc_section, &
278 : pw_pool=pw_pool, rho_set=rho_set1, &
279 : deriv_set=dSet1, deriv_order=1, &
280 0 : calc_potential=.FALSE.)
281 : CALL xc_vxc_pw_create(rho_r=rho_r, rho_g=rho_g, tau=tau, &
282 : vxc_rho=vxc_rho, vxc_tau=vxc_tau, exc=exc, xc_section=xc_section, &
283 : pw_pool=pw_pool, &
284 0 : compute_virial=compute_virial, virial_xc=virial_xc)
285 0 : PRINT *, "did calculate xc (lda)"
286 0 : PRINT *, "about to calculate xc (lsd)"
287 : CALL xc_rho_set_and_dset_create(rho_set=rho_set2, deriv_set=dSet2, &
288 : rho_r=rho2_r, rho_g=rho2_g, tau=tau2, xc_section=xc_section, &
289 : pw_pool=pw_pool, deriv_order=1, &
290 0 : calc_potential=.FALSE.)
291 : CALL xc_vxc_pw_create(rho_r=rho2_r, rho_g=rho2_g, tau=tau2, &
292 : vxc_rho=vxc_rho2, vxc_tau=vxc_tau2, exc=exc2, xc_section=xc_section, &
293 0 : pw_pool=pw_pool, compute_virial=compute_virial, virial_xc=virial_xc)
294 0 : PRINT *, "did calculate xc (new)"
295 0 : PRINT *, "at (0,0,0) rho_r=", rho_r(1)%array(0, 0, 0), &
296 0 : "rho2_r(1)=", rho2_r(1)%array(0, 0, 0), &
297 0 : "rho2_r(2)=", rho2_r(2)%array(0, 0, 0), &
298 0 : "rho_r_sm=", rho_set1%rho(0, 0, 0), "rhoa2_r_sm=", rho_set2%rhoa(0, 0, 0), &
299 0 : "rhob2_r_sm=", rho_set2%rhob(0, 0, 0)
300 : OPEN (unit=120, file="rho.bindata", status="unknown", access='sequential', &
301 0 : form="unformatted", action="write")
302 0 : pot => rho_set1%rho
303 0 : WRITE (unit=120) pot(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(2, 3))
304 0 : CLOSE (unit=120)
305 : OPEN (unit=120, file="rhoa.bindata", status="unknown", access='sequential', &
306 0 : form="unformatted", action="write")
307 0 : pot => rho_set2%rhoa
308 0 : WRITE (unit=120) pot(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(2, 3))
309 0 : CLOSE (unit=120)
310 : OPEN (unit=120, file="rhob.bindata", status="unknown", access='sequential', &
311 0 : form="unformatted", action="write")
312 0 : pot => rho_set2%rhob
313 0 : WRITE (unit=120) pot(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(2, 3))
314 0 : CLOSE (unit=120)
315 : OPEN (unit=120, file="ndrho.bindata", status="unknown", access='sequential', &
316 0 : form="unformatted", action="write")
317 0 : pot => rho_set1%norm_drho
318 0 : WRITE (unit=120) pot(:, :, bo(2, 3))
319 0 : CLOSE (unit=120)
320 : OPEN (unit=120, file="ndrhoa.bindata", status="unknown", access='sequential', &
321 0 : form="unformatted", action="write")
322 0 : pot => rho_set2%norm_drhoa
323 0 : WRITE (unit=120) pot(:, :, bo(2, 3))
324 0 : CLOSE (unit=120)
325 : OPEN (unit=120, file="ndrhob.bindata", status="unknown", access='sequential', &
326 0 : form="unformatted", action="write")
327 0 : pot => rho_set2%norm_drhob
328 0 : WRITE (unit=120) pot(:, :, bo(2, 3))
329 0 : CLOSE (unit=120)
330 0 : IF (rho_set1%has%tau) THEN
331 : OPEN (unit=120, file="tau.bindata", status="unknown", access='sequential', &
332 0 : form="unformatted", action="write")
333 0 : pot => rho_set1%tau
334 0 : WRITE (unit=120) pot(:, :, bo(2, 3))
335 0 : CLOSE (unit=120)
336 : END IF
337 0 : IF (rho_set2%has%tau_spin) THEN
338 : OPEN (unit=120, file="tau_a.bindata", status="unknown", access='sequential', &
339 0 : form="unformatted", action="write")
340 0 : pot => rho_set2%tau_a
341 0 : WRITE (unit=120) pot(:, :, bo(2, 3))
342 0 : CLOSE (unit=120)
343 : OPEN (unit=120, file="tau_v.bindata", status="unknown", access='sequential', &
344 0 : form="unformatted", action="write")
345 0 : pot => rho_set2%tau_b
346 0 : WRITE (unit=120) pot(:, :, bo(2, 3))
347 0 : CLOSE (unit=120)
348 : END IF
349 : OPEN (unit=120, file="vxc.bindata", status="unknown", access='sequential', &
350 0 : form="unformatted", action="write")
351 0 : pot => vxc_rho(1)%array
352 0 : WRITE (unit=120) pot(:, :, bo(2, 3))
353 0 : CLOSE (unit=120)
354 : OPEN (unit=120, file="vxc2.bindata", status="unknown", access='sequential', &
355 0 : form="unformatted", action="write")
356 0 : pot => vxc_rho2(1)%array
357 0 : WRITE (unit=120) pot(:, :, bo(2, 3))
358 0 : CLOSE (unit=120)
359 0 : IF (ASSOCIATED(vxc_tau)) THEN
360 : OPEN (unit=120, file="vxc_tau.bindata", status="unknown", access='sequential', &
361 0 : form="unformatted", action="write")
362 0 : pot => vxc_tau(1)%array
363 0 : WRITE (unit=120) pot(:, :, bo(2, 3))
364 0 : CLOSE (unit=120)
365 : END IF
366 0 : IF (ASSOCIATED(vxc_tau2)) THEN
367 : OPEN (unit=120, file="vxc_tau2_a.bindata", status="unknown", access='sequential', &
368 0 : form="unformatted", action="write")
369 0 : pot => vxc_tau2(1)%array
370 0 : WRITE (unit=120) pot(:, :, bo(2, 3))
371 0 : CLOSE (unit=120)
372 : OPEN (unit=120, file="vxc_tau2_b.bindata", status="unknown", access='sequential', &
373 0 : form="unformatted", action="write")
374 0 : pot => vxc_tau2(2)%array
375 0 : WRITE (unit=120) pot(:, :, bo(2, 3))
376 0 : CLOSE (unit=120)
377 : END IF
378 :
379 0 : PRINT *, "calc diff on vxc"
380 0 : maxDiff = 0.0_dp
381 0 : DO ispin = 1, 1
382 0 : ii = 0
383 0 : DO k = bo(1, 3), bo(2, 3)
384 0 : DO j = bo(1, 2), bo(2, 2)
385 0 : DO i = bo(1, 1), bo(2, 1)
386 0 : ii = ii + 1
387 : diff = ABS(vxc_rho(ispin)%array(i, j, k) - &
388 0 : vxc_rho2(ispin)%array(i, j, k))
389 0 : IF (ii == 1) THEN
390 0 : PRINT *, "vxc", ispin, "=", vxc_rho(ispin)%array(i, j, k), &
391 0 : "vs", vxc_rho2(ispin)%array(i, j, k), "diff=", diff
392 : END IF
393 0 : IF (maxDiff < diff) THEN
394 0 : maxDiff = diff
395 0 : PRINT *, "diff=", diff, " at ", i, ",", j, ",", k, &
396 0 : " spin=", ispin, "rho=", rho_set1%rho(i, j, k), &
397 0 : " ndrho=", rho_set1%norm_drho(i, j, k)
398 : END IF
399 : END DO
400 : END DO
401 : END DO
402 : END DO
403 0 : PRINT *, "diff exc=", ABS(exc - exc2), "diff vxc=", maxdiff
404 : ! CPASSERT(maxdiff<5.e-11)
405 : ! CPASSERT(ABS(exc-exc2)<1.e-14)
406 :
407 0 : IF (ASSOCIATED(vxc_tau)) THEN
408 0 : PRINT *, "calc diff on vxc_tau"
409 0 : maxDiff = 0.0_dp
410 0 : DO ispin = 1, 1
411 0 : ii = 0
412 0 : DO k = bo(1, 3), bo(2, 3)
413 0 : DO j = bo(1, 2), bo(2, 2)
414 0 : DO i = bo(1, 1), bo(2, 1)
415 0 : ii = ii + 1
416 : diff = ABS(vxc_tau(ispin)%array(i, j, k) - &
417 0 : vxc_tau2(ispin)%array(i, j, k))
418 0 : IF (ii == 1) THEN
419 0 : PRINT *, "vxc_tau", ispin, "=", vxc_tau(ispin)%array(i, j, k), &
420 0 : "vs", vxc_tau2(ispin)%array(i, j, k), "diff=", diff
421 : END IF
422 0 : IF (maxDiff < diff) THEN
423 0 : maxDiff = diff
424 0 : PRINT *, "diff=", diff, " at ", i, ",", j, ",", k, &
425 0 : " spin=", ispin, "rho=", rho_set1%rho(i, j, k), &
426 0 : " ndrho=", rho_set1%norm_drho(i, j, k)
427 : END IF
428 : END DO
429 : END DO
430 : END DO
431 : END DO
432 0 : PRINT *, "diff exc=", ABS(exc - exc2), "diff vxc_tau=", maxdiff
433 : END IF
434 0 : deriv_iter => dSet1%derivs
435 0 : DO WHILE (cp_sll_xc_deriv_next(deriv_iter, el_att=deriv))
436 : CALL xc_derivative_get(deriv, order=order, &
437 0 : split_desc=split_desc, deriv_data=pot)
438 0 : SELECT CASE (order)
439 : CASE (0)
440 0 : filename = "e_0.bindata"
441 0 : deriv2 => xc_dset_get_derivative(dSet2, [INTEGER::])
442 : CASE (1)
443 0 : filename = "e_"//TRIM(id_to_desc(split_desc(1)))//".bindata"
444 0 : IF (split_desc(1) == deriv_rho) THEN
445 0 : deriv2 => xc_dset_get_derivative(dSet2, [deriv_rhoa])
446 0 : ELSEIF (split_desc(1) == deriv_tau) THEN
447 0 : deriv2 => xc_dset_get_derivative(dSet2, [deriv_tau_a])
448 0 : ELSEIF (split_desc(1) == deriv_norm_drho) THEN
449 0 : deriv2 => xc_dset_get_derivative(dSet2, [deriv_norm_drhoa])
450 0 : deriv3 => xc_dset_get_derivative(dSet2, [deriv_norm_drho])
451 0 : IF (ASSOCIATED(deriv3)) THEN
452 0 : IF (ASSOCIATED(deriv2)) THEN
453 : CALL xc_derivative_get(deriv2, &
454 0 : deriv_data=pot2)
455 : CALL xc_derivative_get(deriv3, &
456 0 : deriv_data=pot3)
457 0 : pot2 = pot2 + pot3
458 : ELSE
459 : deriv2 => deriv3
460 : END IF
461 0 : NULLIFY (deriv3, pot2, pot3)
462 : END IF
463 : ELSE
464 0 : CPABORT("Unknown derivative variable")
465 : END IF
466 : CASE default
467 0 : CPABORT("Unsupported derivative order")
468 : END SELECT
469 : CALL xc_derivative_get(deriv2, &
470 0 : deriv_data=pot2)
471 0 : PRINT *, "checking ", filename
472 0 : maxDiff = 0.0_dp
473 0 : DO k = bo(1, 3), bo(2, 3)
474 0 : DO j = bo(1, 2), bo(2, 2)
475 0 : DO i = bo(1, 1), bo(2, 1)
476 0 : diff = ABS(pot(i, j, k) - pot2(i, j, k))
477 0 : IF (maxDiff < diff) THEN
478 0 : maxDiff = diff
479 0 : PRINT *, "ediff(", i, j, k, ")=", maxDiff, &
480 0 : "rho=", rho_set1%rho(i, j, k), &
481 0 : "ndrho=", rho_set1%norm_drho(i, j, k)
482 : END IF
483 : END DO
484 : END DO
485 : END DO
486 0 : PRINT *, "maxdiff ", filename, "=", maxDiff
487 : OPEN (unit=120, file=TRIM(filename), status="unknown", &
488 : access='sequential', &
489 0 : form="unformatted")
490 0 : WRITE (unit=120) pot(:, :, bo(2, 3))
491 0 : CLOSE (unit=120)
492 : END DO
493 0 : deriv_iter => dSet2%derivs
494 0 : DO WHILE (cp_sll_xc_deriv_next(deriv_iter, el_att=deriv))
495 : CALL xc_derivative_get(deriv, order=order, &
496 0 : split_desc=split_desc, deriv_data=pot)
497 0 : SELECT CASE (order)
498 : CASE (0)
499 0 : filename = "e_0-2.bindata"
500 : CASE (1)
501 0 : filename = "e_"//TRIM(id_to_desc(split_desc(1)))//"-2.bindata"
502 : CASE default
503 0 : CPABORT("Unsupported derivative order")
504 : END SELECT
505 : OPEN (unit=120, file=TRIM(filename), status="unknown", &
506 : access='sequential', &
507 0 : form="unformatted")
508 0 : WRITE (unit=120) pot(:, :, bo(2, 3))
509 0 : CLOSE (unit=120)
510 : END DO
511 0 : CALL xc_rho_set_release(rho_set1)
512 0 : CALL xc_rho_set_release(rho_set2)
513 0 : CALL xc_dset_release(dSet2)
514 0 : CALL xc_dset_release(dSet1)
515 0 : DO ispin = 1, 2
516 0 : CALL pw_pool%give_back_pw(rho2_r(ispin))
517 0 : CALL pw_pool%give_back_pw(vxc_rho2(ispin))
518 0 : IF (ASSOCIATED(vxc_tau2)) THEN
519 0 : CALL pw_pool%give_back_pw(vxc_tau2(ispin))
520 : END IF
521 : END DO
522 0 : DEALLOCATE (vxc_rho2, rho2_r, rho2_g)
523 0 : IF (ASSOCIATED(vxc_tau2)) THEN
524 0 : DEALLOCATE (vxc_tau2)
525 : END IF
526 :
527 0 : END SUBROUTINE xc_vxc_pw_create_test_lsd
528 :
529 : ! **************************************************************************************************
530 : !> \brief calculates vxc outputting the yz plane of rho, and of the various components
531 : !> of the the derivatives and of vxc
532 : !> \param vxc_rho will contain the v_xc part that depend on rho
533 : !> (if one of the chosen xc functionals has it it is allocated and you
534 : !> are responsible for it)
535 : !> \param vxc_tau will contain the kinetic tau part of v_xc
536 : !> (if one of the chosen xc functionals has it it is allocated and you
537 : !> are responsible for it)
538 : !> \param rho_r the value of the density in the real space
539 : !> \param rho_g value of the density in the g space (needs to be associated
540 : !> only for gradient corrections)
541 : !> \param tau value of the kinetic density tau on the grid (can be null,
542 : !> used only with meta functionals)
543 : !> \param exc the xc energy
544 : !> \param xc_section which functional should be used, and how to do it
545 : !> \param pw_pool the pool for the grids
546 : !> \author Fawzi Mohamed
547 : !> \note
548 : !> for debugging only.
549 : ! **************************************************************************************************
550 0 : SUBROUTINE xc_vxc_pw_create_debug(vxc_rho, vxc_tau, rho_r, rho_g, tau, exc, &
551 : xc_section, pw_pool)
552 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: vxc_rho, vxc_tau
553 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, tau
554 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
555 : REAL(KIND=dp), INTENT(out) :: exc
556 : TYPE(section_vals_type), POINTER :: xc_section
557 : TYPE(pw_pool_type), POINTER :: pw_pool
558 :
559 : LOGICAL, PARAMETER :: compute_virial = .FALSE.
560 :
561 : CHARACTER(len=default_path_length) :: filename
562 0 : INTEGER, DIMENSION(:), POINTER :: split_desc
563 : INTEGER :: i, ispin, j, k, order
564 : INTEGER, DIMENSION(2, 3) :: bo
565 : REAL(KIND=dp), DIMENSION(3, 3) :: virial_xc
566 0 : REAL(kind=dp), DIMENSION(:, :, :), POINTER :: pot
567 : TYPE(cp_logger_type), POINTER :: logger
568 : TYPE(cp_sll_xc_deriv_type), POINTER :: deriv_iter
569 : TYPE(xc_derivative_set_type) :: dSet1
570 : TYPE(xc_derivative_type), POINTER :: deriv
571 : TYPE(xc_rho_set_type) :: rho_set1
572 :
573 0 : NULLIFY (pot, deriv)
574 0 : logger => cp_get_default_logger()
575 :
576 0 : bo = rho_r(1)%pw_grid%bounds_local
577 :
578 : CALL xc_rho_set_and_dset_create(rho_r=rho_r, rho_g=rho_g, &
579 : tau=tau, xc_section=xc_section, &
580 : pw_pool=pw_pool, rho_set=rho_set1, &
581 : deriv_set=dSet1, deriv_order=1, &
582 0 : calc_potential=.FALSE.)
583 :
584 : ! outputs 0,:,: plane
585 0 : IF (bo(1, 1) <= 0 .AND. 0 <= bo(2, 1)) THEN
586 0 : IF (rho_set1%has%rho_spin) THEN
587 : OPEN (unit=120, file="rhoa.bindata", status="unknown", access='sequential', &
588 0 : form="unformatted", action="write")
589 0 : pot => rho_set1%rhoa
590 0 : WRITE (unit=120) pot(0, :, :)
591 0 : CLOSE (unit=120)
592 : OPEN (unit=120, file="rhob.bindata", status="unknown", access='sequential', &
593 0 : form="unformatted", action="write")
594 0 : pot => rho_set1%rhob
595 0 : WRITE (unit=120) pot(0, :, :)
596 0 : CLOSE (unit=120)
597 : END IF
598 0 : IF (rho_set1%has%norm_drho) THEN
599 : OPEN (unit=120, file="ndrho.bindata", status="unknown", access='sequential', &
600 0 : form="unformatted", action="write")
601 0 : pot => rho_set1%norm_drho
602 0 : WRITE (unit=120) pot(0, :, :)
603 0 : CLOSE (unit=120)
604 : END IF
605 0 : IF (rho_set1%has%norm_drho_spin) THEN
606 : OPEN (unit=120, file="ndrhoa.bindata", status="unknown", access='sequential', &
607 0 : form="unformatted", action="write")
608 0 : pot => rho_set1%norm_drhoa
609 0 : WRITE (unit=120) pot(0, :, :)
610 0 : CLOSE (unit=120)
611 : OPEN (unit=120, file="ndrhob.bindata", status="unknown", access='sequential', &
612 0 : form="unformatted", action="write")
613 0 : pot => rho_set1%norm_drhob
614 0 : WRITE (unit=120) pot(0, :, :)
615 0 : CLOSE (unit=120)
616 : END IF
617 0 : IF (rho_set1%has%rho) THEN
618 : OPEN (unit=120, file="rho.bindata", status="unknown", access='sequential', &
619 0 : form="unformatted", action="write")
620 0 : pot => rho_set1%rho
621 0 : WRITE (unit=120) pot(0, :, :)
622 0 : CLOSE (unit=120)
623 : END IF
624 0 : IF (rho_set1%has%tau) THEN
625 : OPEN (unit=120, file="tau.bindata", status="unknown", access='sequential', &
626 0 : form="unformatted", action="write")
627 0 : pot => rho_set1%tau
628 0 : WRITE (unit=120) pot(0, :, :)
629 0 : CLOSE (unit=120)
630 : END IF
631 0 : IF (rho_set1%has%tau_spin) THEN
632 : OPEN (unit=120, file="tau_a.bindata", status="unknown", access='sequential', &
633 0 : form="unformatted", action="write")
634 0 : pot => rho_set1%tau_a
635 0 : WRITE (unit=120) pot(0, :, :)
636 0 : CLOSE (unit=120)
637 : OPEN (unit=120, file="tau_b.bindata", status="unknown", access='sequential', &
638 0 : form="unformatted", action="write")
639 0 : pot => rho_set1%tau_b
640 0 : WRITE (unit=120) pot(0, :, :)
641 0 : CLOSE (unit=120)
642 : END IF
643 :
644 0 : deriv_iter => dSet1%derivs
645 0 : DO WHILE (cp_sll_xc_deriv_next(deriv_iter, el_att=deriv))
646 : CALL xc_derivative_get(deriv, order=order, &
647 0 : split_desc=split_desc, deriv_data=pot)
648 0 : SELECT CASE (order)
649 : CASE (0)
650 0 : filename = "e_0.bindata"
651 : CASE (1)
652 0 : filename = "e_"//TRIM(id_to_desc(split_desc(1)))//".bindata"
653 : CASE default
654 0 : CPABORT("Unsupported derivative order")
655 : END SELECT
656 : OPEN (unit=120, file=TRIM(filename), status="unknown", &
657 : access='sequential', &
658 0 : form="unformatted")
659 0 : WRITE (unit=120) pot(0, :, :)
660 0 : CLOSE (unit=120)
661 : END DO
662 : END IF
663 :
664 : CALL xc_vxc_pw_create(vxc_rho=vxc_rho, vxc_tau=vxc_tau, &
665 : rho_r=rho_r, rho_g=rho_g, tau=tau, &
666 : exc=exc, xc_section=xc_section, &
667 : pw_pool=pw_pool, &
668 0 : compute_virial=compute_virial, virial_xc=virial_xc)
669 :
670 : ! outputs 0,:,: plane
671 0 : IF (bo(1, 1) <= 0 .AND. 0 <= bo(2, 1)) THEN
672 0 : IF (ASSOCIATED(vxc_rho)) THEN
673 0 : DO ispin = 1, SIZE(vxc_rho)
674 0 : WRITE (filename, "('vxc-',i1,'.bindata')") ispin
675 : OPEN (unit=120, file=filename, status="unknown", access='sequential', &
676 0 : form="unformatted", action="write")
677 0 : pot => vxc_rho(ispin)%array
678 0 : WRITE (unit=120) pot(0, :, :)
679 0 : CLOSE (unit=120)
680 :
681 0 : pot => vxc_rho(ispin)%array
682 0 : DO k = bo(1, 3), bo(2, 3)
683 0 : DO j = bo(1, 2), bo(2, 2)
684 0 : DO i = bo(1, 1), bo(2, 1)
685 0 : IF (ABS(pot(i, j, k)) > 10.0_dp) THEN
686 : WRITE (cp_logger_get_default_unit_nr(logger), &
687 : "(' vxc_',i1,'(',i6,',',i6,',',i6,')=',e11.4)", &
688 0 : advance="no") ispin, i, j, k, pot(i, j, k)
689 0 : IF (rho_set1%has%rho_spin) THEN
690 : WRITE (cp_logger_get_default_unit_nr(logger), &
691 : "(' rho=(',e11.4,',',e11.4,')')", advance="no") &
692 0 : rho_set1%rhoa(i, j, k), rho_set1%rhob(i, j, k)
693 0 : ELSE IF (rho_set1%has%rho) THEN
694 : WRITE (cp_logger_get_default_unit_nr(logger), &
695 0 : "(' rho=',e11.4)", advance="no") rho_set1%rho(i, j, k)
696 : END IF
697 0 : IF (rho_set1%has%norm_drho_spin) THEN
698 : WRITE (cp_logger_get_default_unit_nr(logger), &
699 : "(' ndrho=(',e11.4,',',e11.4,')')", advance="no") &
700 0 : rho_set1%norm_drhoa(i, j, k), &
701 0 : rho_set1%norm_drhob(i, j, k)
702 0 : ELSE IF (rho_set1%has%norm_drho) THEN
703 : WRITE (cp_logger_get_default_unit_nr(logger), &
704 0 : "(' ndrho=',e11.4)", advance="no") rho_set1%norm_drho(i, j, k)
705 : END IF
706 0 : IF (rho_set1%has%tau_spin) THEN
707 : WRITE (cp_logger_get_default_unit_nr(logger), &
708 : "(' tau=(',e11.4,',',e11.4,')')", advance="no") &
709 0 : rho_set1%tau_a(i, j, k), &
710 0 : rho_set1%tau_b(i, j, k)
711 0 : ELSE IF (rho_set1%has%tau) THEN
712 : WRITE (cp_logger_get_default_unit_nr(logger), &
713 0 : "(' tau=',e11.4)", advance="no") rho_set1%tau(i, j, k)
714 : END IF
715 :
716 0 : deriv_iter => dSet1%derivs
717 0 : DO WHILE (cp_sll_xc_deriv_next(deriv_iter, el_att=deriv))
718 : CALL xc_derivative_get(deriv, order=order, &
719 0 : split_desc=split_desc, deriv_data=pot)
720 0 : SELECT CASE (order)
721 : CASE (0)
722 0 : filename = " e_0"
723 : CASE (1)
724 0 : filename = " e_"//TRIM(id_to_desc(split_desc(1)))
725 : CASE default
726 0 : CPABORT("Unsupported derivative order")
727 : END SELECT
728 : WRITE (unit=cp_logger_get_default_unit_nr(logger), &
729 : fmt="(a,'=',e11.4)", advance="no") &
730 0 : TRIM(filename), pot(i, j, k)
731 : END DO
732 :
733 : WRITE (cp_logger_get_default_unit_nr(logger), &
734 0 : "()")
735 : END IF
736 : END DO
737 : END DO
738 : END DO
739 : END DO
740 : END IF
741 0 : IF (ASSOCIATED(vxc_tau)) THEN
742 0 : DO ispin = 1, SIZE(vxc_tau)
743 0 : WRITE (filename, "('vxc_tau_',i1,'.bindata')") ispin
744 : OPEN (unit=120, file=filename, status="unknown", access='sequential', &
745 0 : form="unformatted", action="write")
746 0 : pot => vxc_tau(ispin)%array
747 0 : WRITE (unit=120) pot(0, :, :)
748 0 : CLOSE (unit=120)
749 :
750 0 : pot => vxc_tau(ispin)%array
751 0 : DO k = bo(1, 3), bo(2, 3)
752 0 : DO j = bo(1, 2), bo(2, 2)
753 0 : DO i = bo(1, 1), bo(2, 1)
754 0 : IF (ABS(pot(i, j, k)) > 10.0_dp) THEN
755 : WRITE (cp_logger_get_default_unit_nr(logger), &
756 : "('vxc_tau_',i1,'(',i6,',',i6,',',i6,')=',e11.4)", &
757 0 : advance="no") ispin, i, j, k, pot(i, j, k)
758 0 : IF (rho_set1%has%rho_spin) THEN
759 : WRITE (cp_logger_get_default_unit_nr(logger), &
760 : "(' rho=(',e11.4,',',e11.4,')')", advance="no") &
761 0 : rho_set1%rhoa(i, j, k), rho_set1%rhob(i, j, k)
762 0 : ELSE IF (rho_set1%has%rho) THEN
763 : WRITE (cp_logger_get_default_unit_nr(logger), &
764 0 : "(' rho=',e11.4)", advance="no") rho_set1%rho(i, j, k)
765 : END IF
766 0 : IF (rho_set1%has%norm_drho_spin) THEN
767 : WRITE (cp_logger_get_default_unit_nr(logger), &
768 : "(' ndrho=(',e11.4,',',e11.4,')')", advance="no") &
769 0 : rho_set1%norm_drhoa(i, j, k), &
770 0 : rho_set1%norm_drhob(i, j, k)
771 0 : ELSE IF (rho_set1%has%norm_drho) THEN
772 : WRITE (cp_logger_get_default_unit_nr(logger), &
773 0 : "(' ndrho=',e11.4)", advance="no") rho_set1%norm_drho(i, j, k)
774 : END IF
775 0 : IF (rho_set1%has%tau_spin) THEN
776 : WRITE (cp_logger_get_default_unit_nr(logger), &
777 : "(' tau=(',e11.4,',',e11.4,')')", advance="no") &
778 0 : rho_set1%tau_a(i, j, k), &
779 0 : rho_set1%tau_b(i, j, k)
780 0 : ELSE IF (rho_set1%has%tau) THEN
781 : WRITE (cp_logger_get_default_unit_nr(logger), &
782 0 : "(' tau=',e11.4)", advance="no") rho_set1%tau(i, j, k)
783 : END IF
784 :
785 0 : deriv_iter => dSet1%derivs
786 0 : DO WHILE (cp_sll_xc_deriv_next(deriv_iter, el_att=deriv))
787 : CALL xc_derivative_get(deriv, order=order, &
788 0 : split_desc=split_desc, deriv_data=pot)
789 0 : SELECT CASE (order)
790 : CASE (0)
791 0 : filename = " e_0"
792 : CASE (1)
793 0 : filename = " e_"//TRIM(id_to_desc(split_desc(1)))
794 : CASE default
795 0 : CPABORT("Unsupported derivative order")
796 : END SELECT
797 : WRITE (unit=cp_logger_get_default_unit_nr(logger), &
798 : fmt="(a,'=',e11.4)", advance="no") &
799 0 : TRIM(filename), pot(i, j, k)
800 : END DO
801 :
802 0 : WRITE (cp_logger_get_default_unit_nr(logger), "()")
803 : END IF
804 : END DO
805 : END DO
806 : END DO
807 : END DO
808 : END IF
809 : END IF
810 :
811 0 : CALL xc_dset_release(dSet1)
812 0 : CALL xc_rho_set_release(rho_set1)
813 0 : END SUBROUTINE xc_vxc_pw_create_debug
814 :
815 : ! **************************************************************************************************
816 : !> \brief creates a xc_rho_set and a derivative set containing the derivatives
817 : !> of the functionals with the given deriv_order.
818 : !> \param rho_set will contain the rho set
819 : !> \param deriv_set will contain the derivatives
820 : !> \param deriv_order the order of the requested derivatives. If positive
821 : !> 0:deriv_order are calculated, if negative only -deriv_order is
822 : !> guaranteed to be valid. Orders not requested might be present,
823 : !> but might contain garbage.
824 : !> \param rho_r the value of the density in the real space
825 : !> \param rho_g value of the density in the g space (can be null, used only
826 : !> without smoothing of rho or deriv)
827 : !> \param tau value of the kinetic density tau on the grid (can be null,
828 : !> used only with meta functionals)
829 : !> \param xc_section the section describing the functional to use
830 : !> \param pw_pool the pool for the grids
831 : !> \param calc_potential if the basic components of the arguments
832 : !> should be kept in rho set (a basic component is for example drho
833 : !> when with lda a functional needs norm_drho)
834 : !> \author fawzi
835 : !> \note
836 : !> if any of the functionals is gradient corrected the full gradient is
837 : !> added to the rho set
838 : ! **************************************************************************************************
839 260302 : SUBROUTINE xc_rho_set_and_dset_create(rho_set, deriv_set, deriv_order, &
840 : rho_r, rho_g, tau, xc_section, pw_pool, &
841 : calc_potential)
842 :
843 : TYPE(xc_rho_set_type) :: rho_set
844 : TYPE(xc_derivative_set_type) :: deriv_set
845 : INTEGER, INTENT(in) :: deriv_order
846 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, tau
847 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
848 : TYPE(section_vals_type), POINTER :: xc_section
849 : TYPE(pw_pool_type), POINTER :: pw_pool
850 : LOGICAL, INTENT(in) :: calc_potential
851 :
852 : CHARACTER(len=*), PARAMETER :: routineN = 'xc_rho_set_and_dset_create'
853 :
854 : INTEGER :: handle, nspins
855 : LOGICAL :: lsd
856 : TYPE(section_vals_type), POINTER :: xc_fun_sections
857 :
858 130151 : CALL timeset(routineN, handle)
859 :
860 130151 : CPASSERT(ASSOCIATED(pw_pool))
861 :
862 130151 : nspins = SIZE(rho_r)
863 130151 : lsd = (nspins /= 1)
864 :
865 130151 : xc_fun_sections => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
866 :
867 130151 : CALL xc_dset_create(deriv_set, pw_pool)
868 :
869 : CALL xc_rho_set_create(rho_set, &
870 : rho_r(1)%pw_grid%bounds_local, &
871 : rho_cutoff=section_get_rval(xc_section, "density_cutoff"), &
872 : drho_cutoff=section_get_rval(xc_section, "gradient_cutoff"), &
873 130151 : tau_cutoff=section_get_rval(xc_section, "tau_cutoff"))
874 :
875 : CALL xc_rho_set_update(rho_set, rho_r, rho_g, tau, &
876 : xc_functionals_get_needs(xc_fun_sections, lsd, calc_potential), &
877 : section_get_ival(xc_section, "XC_GRID%XC_DERIV"), &
878 : section_get_ival(xc_section, "XC_GRID%XC_SMOOTH_RHO"), &
879 130151 : pw_pool)
880 :
881 : CALL xc_functionals_eval(xc_fun_sections, &
882 : lsd=lsd, &
883 : rho_set=rho_set, &
884 : deriv_set=deriv_set, &
885 130151 : deriv_order=deriv_order)
886 :
887 130151 : CALL divide_by_norm_drho(deriv_set, rho_set, lsd)
888 :
889 130151 : CALL timestop(handle)
890 :
891 130151 : END SUBROUTINE xc_rho_set_and_dset_create
892 :
893 : ! **************************************************************************************************
894 : !> \brief smooths the cutoff on rho with a function smoothderiv_rho that is 0
895 : !> for rho<rho_cutoff and 1 for rho>rho_cutoff*rho_smooth_cutoff_range:
896 : !> E= integral e_0*smoothderiv_rho => dE/d...= de/d... * smooth,
897 : !> dE/drho = de/drho * smooth + e_0 * dsmooth/drho
898 : !> \param pot the potential to smooth
899 : !> \param rho , rhoa,rhob: the value of the density (used to apply the cutoff)
900 : !> \param rhoa ...
901 : !> \param rhob ...
902 : !> \param rho_cutoff the value at whch the cutoff function must go to 0
903 : !> \param rho_smooth_cutoff_range range of the smoothing
904 : !> \param e_0 value of e_0, if given it is assumed that pot is the derivative
905 : !> wrt. to rho, and needs the dsmooth*e_0 contribution
906 : !> \param e_0_scale_factor ...
907 : !> \author Fawzi Mohamed
908 : ! **************************************************************************************************
909 250384 : SUBROUTINE smooth_cutoff(pot, rho, rhoa, rhob, rho_cutoff, &
910 : rho_smooth_cutoff_range, e_0, e_0_scale_factor)
911 : REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN), &
912 : POINTER :: pot, rho, rhoa, rhob
913 : REAL(kind=dp), INTENT(in) :: rho_cutoff, rho_smooth_cutoff_range
914 : REAL(kind=dp), DIMENSION(:, :, :), OPTIONAL, &
915 : POINTER :: e_0
916 : REAL(kind=dp), INTENT(in), OPTIONAL :: e_0_scale_factor
917 :
918 : INTEGER :: i, j, k
919 : INTEGER, DIMENSION(2, 3) :: bo
920 : REAL(kind=dp) :: my_e_0_scale_factor, my_rho, my_rho_n, my_rho_n2, rho_smooth_cutoff, &
921 : rho_smooth_cutoff_2, rho_smooth_cutoff_range_2
922 :
923 250384 : CPASSERT(ASSOCIATED(pot))
924 1001536 : bo(1, :) = LBOUND(pot)
925 1001536 : bo(2, :) = UBOUND(pot)
926 250384 : my_e_0_scale_factor = 1.0_dp
927 250384 : IF (PRESENT(e_0_scale_factor)) my_e_0_scale_factor = e_0_scale_factor
928 250384 : rho_smooth_cutoff = rho_cutoff*rho_smooth_cutoff_range
929 250384 : rho_smooth_cutoff_2 = (rho_cutoff + rho_smooth_cutoff)/2
930 250384 : rho_smooth_cutoff_range_2 = rho_smooth_cutoff_2 - rho_cutoff
931 :
932 250384 : IF (rho_smooth_cutoff_range > 0.0_dp) THEN
933 2 : IF (PRESENT(e_0)) THEN
934 0 : CPASSERT(ASSOCIATED(e_0))
935 0 : IF (ASSOCIATED(rho)) THEN
936 : !$OMP PARALLEL DO DEFAULT(NONE) &
937 : !$OMP SHARED(bo,e_0,pot,rho,rho_cutoff,rho_smooth_cutoff,rho_smooth_cutoff_2, &
938 : !$OMP rho_smooth_cutoff_range_2,my_e_0_scale_factor) &
939 : !$OMP PRIVATE(k,j,i,my_rho,my_rho_n,my_rho_n2) &
940 0 : !$OMP COLLAPSE(3)
941 : DO k = bo(1, 3), bo(2, 3)
942 : DO j = bo(1, 2), bo(2, 2)
943 : DO i = bo(1, 1), bo(2, 1)
944 : my_rho = rho(i, j, k)
945 : IF (my_rho < rho_smooth_cutoff) THEN
946 : IF (my_rho < rho_cutoff) THEN
947 : pot(i, j, k) = 0.0_dp
948 : ELSEIF (my_rho < rho_smooth_cutoff_2) THEN
949 : my_rho_n = (my_rho - rho_cutoff)/rho_smooth_cutoff_range_2
950 : my_rho_n2 = my_rho_n*my_rho_n
951 : pot(i, j, k) = pot(i, j, k)* &
952 : my_rho_n2*(my_rho_n - 0.5_dp*my_rho_n2) + &
953 : my_e_0_scale_factor*e_0(i, j, k)* &
954 : my_rho_n2*(3.0_dp - 2.0_dp*my_rho_n) &
955 : /rho_smooth_cutoff_range_2
956 : ELSE
957 : my_rho_n = 2.0_dp - (my_rho - rho_cutoff)/rho_smooth_cutoff_range_2
958 : my_rho_n2 = my_rho_n*my_rho_n
959 : pot(i, j, k) = pot(i, j, k)* &
960 : (1.0_dp - my_rho_n2*(my_rho_n - 0.5_dp*my_rho_n2)) &
961 : + my_e_0_scale_factor*e_0(i, j, k)* &
962 : my_rho_n2*(3.0_dp - 2.0_dp*my_rho_n) &
963 : /rho_smooth_cutoff_range_2
964 : END IF
965 : END IF
966 : END DO
967 : END DO
968 : END DO
969 : !$OMP END PARALLEL DO
970 : ELSE
971 : !$OMP PARALLEL DO DEFAULT(NONE) &
972 : !$OMP SHARED(bo,pot,e_0,rhoa,rhob,rho_cutoff,rho_smooth_cutoff,rho_smooth_cutoff_2, &
973 : !$OMP rho_smooth_cutoff_range_2,my_e_0_scale_factor) &
974 : !$OMP PRIVATE(k,j,i,my_rho,my_rho_n,my_rho_n2) &
975 0 : !$OMP COLLAPSE(3)
976 : DO k = bo(1, 3), bo(2, 3)
977 : DO j = bo(1, 2), bo(2, 2)
978 : DO i = bo(1, 1), bo(2, 1)
979 : my_rho = rhoa(i, j, k) + rhob(i, j, k)
980 : IF (my_rho < rho_smooth_cutoff) THEN
981 : IF (my_rho < rho_cutoff) THEN
982 : pot(i, j, k) = 0.0_dp
983 : ELSEIF (my_rho < rho_smooth_cutoff_2) THEN
984 : my_rho_n = (my_rho - rho_cutoff)/rho_smooth_cutoff_range_2
985 : my_rho_n2 = my_rho_n*my_rho_n
986 : pot(i, j, k) = pot(i, j, k)* &
987 : my_rho_n2*(my_rho_n - 0.5_dp*my_rho_n2) + &
988 : my_e_0_scale_factor*e_0(i, j, k)* &
989 : my_rho_n2*(3.0_dp - 2.0_dp*my_rho_n) &
990 : /rho_smooth_cutoff_range_2
991 : ELSE
992 : my_rho_n = 2.0_dp - (my_rho - rho_cutoff)/rho_smooth_cutoff_range_2
993 : my_rho_n2 = my_rho_n*my_rho_n
994 : pot(i, j, k) = pot(i, j, k)* &
995 : (1.0_dp - my_rho_n2*(my_rho_n - 0.5_dp*my_rho_n2)) &
996 : + my_e_0_scale_factor*e_0(i, j, k)* &
997 : my_rho_n2*(3.0_dp - 2.0_dp*my_rho_n) &
998 : /rho_smooth_cutoff_range_2
999 : END IF
1000 : END IF
1001 : END DO
1002 : END DO
1003 : END DO
1004 : !$OMP END PARALLEL DO
1005 : END IF
1006 : ELSE
1007 2 : IF (ASSOCIATED(rho)) THEN
1008 : !$OMP PARALLEL DO DEFAULT(NONE) &
1009 : !$OMP SHARED(bo,pot,rho_cutoff,rho_smooth_cutoff,rho_smooth_cutoff_2, &
1010 : !$OMP rho_smooth_cutoff_range_2,rho) &
1011 : !$OMP PRIVATE(k,j,i,my_rho,my_rho_n,my_rho_n2) &
1012 2 : !$OMP COLLAPSE(3)
1013 : DO k = bo(1, 3), bo(2, 3)
1014 : DO j = bo(1, 2), bo(2, 2)
1015 : DO i = bo(1, 1), bo(2, 1)
1016 : my_rho = rho(i, j, k)
1017 : IF (my_rho < rho_smooth_cutoff) THEN
1018 : IF (my_rho < rho_cutoff) THEN
1019 : pot(i, j, k) = 0.0_dp
1020 : ELSEIF (my_rho < rho_smooth_cutoff_2) THEN
1021 : my_rho_n = (my_rho - rho_cutoff)/rho_smooth_cutoff_range_2
1022 : my_rho_n2 = my_rho_n*my_rho_n
1023 : pot(i, j, k) = pot(i, j, k)* &
1024 : my_rho_n2*(my_rho_n - 0.5_dp*my_rho_n2)
1025 : ELSE
1026 : my_rho_n = 2.0_dp - (my_rho - rho_cutoff)/rho_smooth_cutoff_range_2
1027 : my_rho_n2 = my_rho_n*my_rho_n
1028 : pot(i, j, k) = pot(i, j, k)* &
1029 : (1.0_dp - my_rho_n2*(my_rho_n - 0.5_dp*my_rho_n2))
1030 : END IF
1031 : END IF
1032 : END DO
1033 : END DO
1034 : END DO
1035 : !$OMP END PARALLEL DO
1036 : ELSE
1037 : !$OMP PARALLEL DO DEFAULT(NONE) &
1038 : !$OMP SHARED(bo,pot,rho_cutoff,rho_smooth_cutoff,rho_smooth_cutoff_2, &
1039 : !$OMP rho_smooth_cutoff_range_2,rhoa,rhob) &
1040 : !$OMP PRIVATE(k,j,i,my_rho,my_rho_n,my_rho_n2) &
1041 0 : !$OMP COLLAPSE(3)
1042 : DO k = bo(1, 3), bo(2, 3)
1043 : DO j = bo(1, 2), bo(2, 2)
1044 : DO i = bo(1, 1), bo(2, 1)
1045 : my_rho = rhoa(i, j, k) + rhob(i, j, k)
1046 : IF (my_rho < rho_smooth_cutoff) THEN
1047 : IF (my_rho < rho_cutoff) THEN
1048 : pot(i, j, k) = 0.0_dp
1049 : ELSEIF (my_rho < rho_smooth_cutoff_2) THEN
1050 : my_rho_n = (my_rho - rho_cutoff)/rho_smooth_cutoff_range_2
1051 : my_rho_n2 = my_rho_n*my_rho_n
1052 : pot(i, j, k) = pot(i, j, k)* &
1053 : my_rho_n2*(my_rho_n - 0.5_dp*my_rho_n2)
1054 : ELSE
1055 : my_rho_n = 2.0_dp - (my_rho - rho_cutoff)/rho_smooth_cutoff_range_2
1056 : my_rho_n2 = my_rho_n*my_rho_n
1057 : pot(i, j, k) = pot(i, j, k)* &
1058 : (1.0_dp - my_rho_n2*(my_rho_n - 0.5_dp*my_rho_n2))
1059 : END IF
1060 : END IF
1061 : END DO
1062 : END DO
1063 : END DO
1064 : !$OMP END PARALLEL DO
1065 : END IF
1066 : END IF
1067 : END IF
1068 250384 : END SUBROUTINE smooth_cutoff
1069 :
1070 64 : SUBROUTINE calc_xc_density(pot, rho, rho_cutoff)
1071 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: pot
1072 : TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(INOUT) :: rho
1073 : REAL(kind=dp), INTENT(in) :: rho_cutoff
1074 :
1075 : INTEGER :: i, j, k, nspins
1076 : INTEGER, DIMENSION(2, 3) :: bo
1077 : REAL(kind=dp) :: eps1, eps2, my_rho, my_pot
1078 :
1079 256 : bo(1, :) = LBOUND(pot%array)
1080 256 : bo(2, :) = UBOUND(pot%array)
1081 64 : nspins = SIZE(rho)
1082 :
1083 64 : eps1 = rho_cutoff*1.E-4_dp
1084 64 : eps2 = rho_cutoff
1085 :
1086 3160 : DO k = bo(1, 3), bo(2, 3)
1087 161272 : DO j = bo(1, 2), bo(2, 2)
1088 4529376 : DO i = bo(1, 1), bo(2, 1)
1089 4368168 : my_pot = pot%array(i, j, k)
1090 4368168 : IF (nspins == 2) THEN
1091 339714 : my_rho = rho(1)%array(i, j, k) + rho(2)%array(i, j, k)
1092 : ELSE
1093 4028454 : my_rho = rho(1)%array(i, j, k)
1094 : END IF
1095 4526280 : IF (my_rho > eps1) THEN
1096 4139616 : pot%array(i, j, k) = my_pot/my_rho
1097 228552 : ELSE IF (my_rho < eps2) THEN
1098 228552 : pot%array(i, j, k) = 0.0_dp
1099 : ELSE
1100 0 : pot%array(i, j, k) = MIN(my_pot/my_rho, my_rho**(1._dp/3._dp))
1101 : END IF
1102 : END DO
1103 : END DO
1104 : END DO
1105 :
1106 64 : END SUBROUTINE calc_xc_density
1107 :
1108 : ! **************************************************************************************************
1109 : !> \brief Exchange and Correlation functional calculations
1110 : !> \param vxc_rho will contain the v_xc part that depend on rho
1111 : !> (if one of the chosen xc functionals has it it is allocated and you
1112 : !> are responsible for it)
1113 : !> \param vxc_tau will contain the kinetic tau part of v_xc
1114 : !> (if one of the chosen xc functionals has it it is allocated and you
1115 : !> are responsible for it)
1116 : !> \param exc the xc energy
1117 : !> \param rho_r the value of the density in the real space
1118 : !> \param rho_g value of the density in the g space (needs to be associated
1119 : !> only for gradient corrections)
1120 : !> \param tau value of the kinetic density tau on the grid (can be null,
1121 : !> used only with meta functionals)
1122 : !> \param xc_section which functional to calculate, and how to do it
1123 : !> \param pw_pool the pool for the grids
1124 : !> \param compute_virial ...
1125 : !> \param virial_xc ...
1126 : !> \param exc_r the value of the xc functional in the real space
1127 : !> \par History
1128 : !> JGH (13-Jun-2002): adaptation to new functionals
1129 : !> Fawzi (11.2002): drho_g(1:3)->drho_g
1130 : !> Fawzi (1.2003). lsd version
1131 : !> Fawzi (11.2003): version using the new xc interface
1132 : !> Fawzi (03.2004): fft free for smoothed density and derivs, gga lsd
1133 : !> Fawzi (04.2004): metafunctionals
1134 : !> mguidon (12.2008) : laplace functionals
1135 : !> \author fawzi; based LDA version of JGH, based on earlier version of apsi
1136 : !> \note
1137 : !> Beware: some really dirty pointer handling!
1138 : !> energy should be kept consistent with xc_exc_calc
1139 : ! **************************************************************************************************
1140 110240 : SUBROUTINE xc_vxc_pw_create(vxc_rho, vxc_tau, exc, rho_r, rho_g, tau, xc_section, &
1141 : pw_pool, compute_virial, virial_xc, exc_r)
1142 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: vxc_rho, vxc_tau
1143 : REAL(KIND=dp), INTENT(out) :: exc
1144 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, tau
1145 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
1146 : TYPE(section_vals_type), POINTER :: xc_section
1147 : TYPE(pw_pool_type), POINTER :: pw_pool
1148 : LOGICAL :: compute_virial
1149 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT) :: virial_xc
1150 : TYPE(pw_r3d_rs_type), INTENT(INOUT), OPTIONAL :: exc_r
1151 :
1152 : CHARACTER(len=*), PARAMETER :: routineN = 'xc_vxc_pw_create'
1153 : INTEGER, DIMENSION(2), PARAMETER :: norm_drho_spin_name = [deriv_norm_drhoa, deriv_norm_drhob]
1154 :
1155 : INTEGER :: handle, idir, ispin, jdir, &
1156 : npoints, nspins, &
1157 : xc_deriv_method_id, xc_rho_smooth_id, deriv_id
1158 : INTEGER, DIMENSION(2, 3) :: bo
1159 : LOGICAL :: dealloc_pw_to_deriv, has_laplace, &
1160 : has_tau, lsd, use_virial, has_gradient, &
1161 : has_derivs, has_rho, dealloc_pw_to_deriv_rho
1162 : REAL(KIND=dp) :: density_smooth_cut_range, drho_cutoff, &
1163 : rho_cutoff
1164 110240 : REAL(kind=dp), DIMENSION(:, :, :), POINTER :: deriv_data, norm_drho, norm_drho_spin, &
1165 220480 : rho, rhoa, rhob
1166 : TYPE(cp_sll_xc_deriv_type), POINTER :: pos
1167 : TYPE(pw_grid_type), POINTER :: pw_grid
1168 771680 : TYPE(pw_r3d_rs_type), DIMENSION(3) :: pw_to_deriv, pw_to_deriv_rho
1169 : TYPE(pw_c1d_gs_type) :: tmp_g, vxc_g
1170 : TYPE(pw_r3d_rs_type) :: v_drho_r, virial_pw
1171 : TYPE(xc_derivative_set_type) :: deriv_set
1172 : TYPE(xc_derivative_type), POINTER :: deriv_att
1173 : TYPE(xc_rho_set_type) :: rho_set
1174 :
1175 110240 : CALL timeset(routineN, handle)
1176 110240 : NULLIFY (norm_drho_spin, norm_drho, pos)
1177 :
1178 110240 : pw_grid => rho_r(1)%pw_grid
1179 :
1180 110240 : CPASSERT(ASSOCIATED(xc_section))
1181 110240 : CPASSERT(ASSOCIATED(pw_pool))
1182 110240 : CPASSERT(.NOT. ASSOCIATED(vxc_rho))
1183 110240 : CPASSERT(.NOT. ASSOCIATED(vxc_tau))
1184 110240 : nspins = SIZE(rho_r)
1185 110240 : lsd = (nspins /= 1)
1186 110240 : IF (lsd) THEN
1187 22403 : CPASSERT(nspins == 2)
1188 : END IF
1189 :
1190 110240 : use_virial = compute_virial
1191 110240 : virial_xc = 0.0_dp
1192 :
1193 1102400 : bo = rho_r(1)%pw_grid%bounds_local
1194 110240 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
1195 :
1196 : ! calculate the potential derivatives
1197 : CALL xc_rho_set_and_dset_create(rho_set=rho_set, deriv_set=deriv_set, &
1198 : deriv_order=1, rho_r=rho_r, rho_g=rho_g, tau=tau, &
1199 : xc_section=xc_section, &
1200 : pw_pool=pw_pool, &
1201 110240 : calc_potential=.TRUE.)
1202 :
1203 : CALL section_vals_val_get(xc_section, "XC_GRID%XC_DERIV", &
1204 110240 : i_val=xc_deriv_method_id)
1205 : CALL section_vals_val_get(xc_section, "XC_GRID%XC_SMOOTH_RHO", &
1206 110240 : i_val=xc_rho_smooth_id)
1207 : CALL section_vals_val_get(xc_section, "DENSITY_SMOOTH_CUTOFF_RANGE", &
1208 110240 : r_val=density_smooth_cut_range)
1209 :
1210 : CALL xc_rho_set_get(rho_set, rho_cutoff=rho_cutoff, &
1211 110240 : drho_cutoff=drho_cutoff)
1212 :
1213 110240 : CALL check_for_derivatives(deriv_set, lsd, has_rho, has_gradient, has_tau, has_laplace)
1214 : ! check for unknown derivatives
1215 110240 : has_derivs = has_rho .OR. has_gradient .OR. has_tau .OR. has_laplace
1216 :
1217 463363 : ALLOCATE (vxc_rho(nspins))
1218 :
1219 : CALL xc_rho_set_get(rho_set, rho=rho, rhoa=rhoa, rhob=rhob, &
1220 110240 : can_return_null=.TRUE.)
1221 :
1222 : ! recover the vxc arrays
1223 110240 : IF (lsd) THEN
1224 22403 : CALL xc_dset_recover_pw(deriv_set, [deriv_rhoa], vxc_rho(1), pw_grid, pw_pool)
1225 22403 : CALL xc_dset_recover_pw(deriv_set, [deriv_rhob], vxc_rho(2), pw_grid, pw_pool)
1226 :
1227 : ELSE
1228 87837 : CALL xc_dset_recover_pw(deriv_set, [deriv_rho], vxc_rho(1), pw_grid, pw_pool)
1229 : END IF
1230 :
1231 110240 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_norm_drho])
1232 110240 : IF (ASSOCIATED(deriv_att)) THEN
1233 58723 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
1234 :
1235 : CALL xc_rho_set_get(rho_set, norm_drho=norm_drho, &
1236 : rho_cutoff=rho_cutoff, &
1237 : drho_cutoff=drho_cutoff, &
1238 58723 : can_return_null=.TRUE.)
1239 58723 : CALL xc_rho_set_recover_pw(rho_set, pw_grid, pw_pool, dealloc_pw_to_deriv_rho, drho=pw_to_deriv_rho)
1240 :
1241 58723 : CPASSERT(ASSOCIATED(deriv_data))
1242 58723 : IF (use_virial) THEN
1243 1526 : CALL pw_pool%create_pw(virial_pw)
1244 1526 : CALL pw_zero(virial_pw)
1245 6104 : DO idir = 1, 3
1246 4578 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(virial_pw,pw_to_deriv_rho,deriv_data,idir)
1247 : virial_pw%array(:, :, :) = pw_to_deriv_rho(idir)%array(:, :, :)*deriv_data(:, :, :)
1248 : !$OMP END PARALLEL WORKSHARE
1249 15260 : DO jdir = 1, idir
1250 : virial_xc(idir, jdir) = -pw_grid%dvol* &
1251 : accurate_dot_product(virial_pw%array(:, :, :), &
1252 9156 : pw_to_deriv_rho(jdir)%array(:, :, :))
1253 13734 : virial_xc(jdir, idir) = virial_xc(idir, jdir)
1254 : END DO
1255 : END DO
1256 1526 : CALL pw_pool%give_back_pw(virial_pw)
1257 : END IF ! use_virial
1258 234892 : DO idir = 1, 3
1259 176169 : CPASSERT(ASSOCIATED(pw_to_deriv_rho(idir)%array))
1260 234892 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(deriv_data,pw_to_deriv_rho,idir)
1261 : pw_to_deriv_rho(idir)%array(:, :, :) = pw_to_deriv_rho(idir)%array(:, :, :)*deriv_data(:, :, :)
1262 : !$OMP END PARALLEL WORKSHARE
1263 : END DO
1264 :
1265 : ! Deallocate pw to save memory
1266 58723 : CALL pw_pool%give_back_cr3d(deriv_att%deriv_data)
1267 :
1268 : END IF
1269 :
1270 110240 : IF ((has_gradient .AND. xc_requires_tmp_g(xc_deriv_method_id)) .OR. pw_grid%spherical) THEN
1271 57997 : CALL pw_pool%create_pw(vxc_g)
1272 57997 : IF (.NOT. pw_grid%spherical) THEN
1273 57997 : CALL pw_pool%create_pw(tmp_g)
1274 : END IF
1275 : END IF
1276 :
1277 242883 : DO ispin = 1, nspins
1278 :
1279 132643 : IF (lsd) THEN
1280 44806 : IF (ispin == 1) THEN
1281 : CALL xc_rho_set_get(rho_set, norm_drhoa=norm_drho_spin, &
1282 22403 : can_return_null=.TRUE.)
1283 22403 : IF (ASSOCIATED(norm_drho_spin)) CALL xc_rho_set_recover_pw( &
1284 13968 : rho_set, pw_grid, pw_pool, dealloc_pw_to_deriv, drhoa=pw_to_deriv)
1285 : ELSE
1286 : CALL xc_rho_set_get(rho_set, norm_drhob=norm_drho_spin, &
1287 22403 : can_return_null=.TRUE.)
1288 22403 : IF (ASSOCIATED(norm_drho_spin)) CALL xc_rho_set_recover_pw( &
1289 13968 : rho_set, pw_grid, pw_pool, dealloc_pw_to_deriv, drhob=pw_to_deriv)
1290 : END IF
1291 :
1292 89612 : deriv_att => xc_dset_get_derivative(deriv_set, [norm_drho_spin_name(ispin)])
1293 44806 : IF (ASSOCIATED(deriv_att)) THEN
1294 : CPASSERT(lsd)
1295 27936 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
1296 :
1297 27936 : IF (use_virial) THEN
1298 100 : CALL pw_pool%create_pw(virial_pw)
1299 100 : CALL pw_zero(virial_pw)
1300 400 : DO idir = 1, 3
1301 300 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(deriv_data,pw_to_deriv,virial_pw,idir)
1302 : virial_pw%array(:, :, :) = pw_to_deriv(idir)%array(:, :, :)*deriv_data(:, :, :)
1303 : !$OMP END PARALLEL WORKSHARE
1304 1000 : DO jdir = 1, idir
1305 : virial_xc(idir, jdir) = virial_xc(idir, jdir) - pw_grid%dvol* &
1306 : accurate_dot_product(virial_pw%array(:, :, :), &
1307 600 : pw_to_deriv(jdir)%array(:, :, :))
1308 900 : virial_xc(jdir, idir) = virial_xc(idir, jdir)
1309 : END DO
1310 : END DO
1311 100 : CALL pw_pool%give_back_pw(virial_pw)
1312 : END IF ! use_virial
1313 :
1314 111744 : DO idir = 1, 3
1315 111744 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(deriv_data,idir,pw_to_deriv)
1316 : pw_to_deriv(idir)%array(:, :, :) = deriv_data(:, :, :)*pw_to_deriv(idir)%array(:, :, :)
1317 : !$OMP END PARALLEL WORKSHARE
1318 : END DO
1319 : END IF ! deriv_att
1320 :
1321 : END IF ! LSD
1322 :
1323 132643 : IF (ASSOCIATED(pw_to_deriv_rho(1)%array)) THEN
1324 71997 : IF (.NOT. ASSOCIATED(pw_to_deriv(1)%array)) THEN
1325 45449 : pw_to_deriv = pw_to_deriv_rho
1326 45449 : dealloc_pw_to_deriv = ((.NOT. lsd) .OR. (ispin == 2))
1327 45449 : dealloc_pw_to_deriv = dealloc_pw_to_deriv .AND. dealloc_pw_to_deriv_rho
1328 : ELSE
1329 : ! This branch is called in case of open-shell systems
1330 : ! Add the contributions from norm_drho and norm_drho_spin
1331 106192 : DO idir = 1, 3
1332 79644 : CALL pw_axpy(pw_to_deriv_rho(idir), pw_to_deriv(idir))
1333 106192 : IF (ispin == 2) THEN
1334 39822 : IF (dealloc_pw_to_deriv_rho) THEN
1335 39822 : CALL pw_pool%give_back_pw(pw_to_deriv_rho(idir))
1336 : END IF
1337 : END IF
1338 : END DO
1339 : END IF
1340 : END IF
1341 :
1342 132643 : IF (ASSOCIATED(pw_to_deriv(1)%array)) THEN
1343 293540 : DO idir = 1, 3
1344 293540 : CALL pw_scale(pw_to_deriv(idir), -1.0_dp)
1345 : END DO
1346 :
1347 73385 : CALL xc_pw_divergence(xc_deriv_method_id, pw_to_deriv, tmp_g, vxc_g, vxc_rho(ispin))
1348 :
1349 73385 : IF (dealloc_pw_to_deriv) THEN
1350 293540 : DO idir = 1, 3
1351 293540 : CALL pw_pool%give_back_pw(pw_to_deriv(idir))
1352 : END DO
1353 : END IF
1354 : END IF
1355 :
1356 : ! Add laplace part to vxc_rho
1357 132643 : IF (has_laplace) THEN
1358 986 : IF (lsd) THEN
1359 416 : IF (ispin == 1) THEN
1360 : deriv_id = deriv_laplace_rhoa
1361 : ELSE
1362 208 : deriv_id = deriv_laplace_rhob
1363 : END IF
1364 : ELSE
1365 : deriv_id = deriv_laplace_rho
1366 : END IF
1367 :
1368 1972 : CALL xc_dset_recover_pw(deriv_set, [deriv_id], pw_to_deriv(1), pw_grid)
1369 :
1370 986 : IF (use_virial) CALL virial_laplace(rho_r(ispin), pw_pool, virial_xc, pw_to_deriv(1)%array)
1371 :
1372 986 : CALL xc_pw_laplace(pw_to_deriv(1), pw_pool, xc_deriv_method_id)
1373 :
1374 986 : CALL pw_axpy(pw_to_deriv(1), vxc_rho(ispin))
1375 :
1376 986 : CALL pw_pool%give_back_pw(pw_to_deriv(1))
1377 : END IF
1378 :
1379 132643 : IF (pw_grid%spherical) THEN
1380 : ! filter vxc
1381 0 : CALL pw_transfer(vxc_rho(ispin), vxc_g)
1382 0 : CALL pw_transfer(vxc_g, vxc_rho(ispin))
1383 : END IF
1384 : CALL smooth_cutoff(pot=vxc_rho(ispin)%array, rho=rho, rhoa=rhoa, rhob=rhob, &
1385 : rho_cutoff=rho_cutoff*density_smooth_cut_range, &
1386 132643 : rho_smooth_cutoff_range=density_smooth_cut_range)
1387 :
1388 132643 : v_drho_r = vxc_rho(ispin)
1389 132643 : CALL pw_pool%create_pw(vxc_rho(ispin))
1390 132643 : CALL xc_pw_smooth(v_drho_r, vxc_rho(ispin), xc_rho_smooth_id)
1391 242883 : CALL pw_pool%give_back_pw(v_drho_r)
1392 : END DO
1393 :
1394 110240 : CALL pw_pool%give_back_pw(vxc_g)
1395 110240 : CALL pw_pool%give_back_pw(tmp_g)
1396 :
1397 : ! 0-deriv -> value of exc
1398 : ! this has to be kept consistent with xc_exc_calc
1399 110240 : IF (has_derivs) THEN
1400 109920 : CALL xc_dset_recover_pw(deriv_set, [INTEGER::], v_drho_r, pw_grid)
1401 :
1402 : CALL smooth_cutoff(pot=v_drho_r%array, rho=rho, rhoa=rhoa, rhob=rhob, &
1403 : rho_cutoff=rho_cutoff, &
1404 109920 : rho_smooth_cutoff_range=density_smooth_cut_range)
1405 :
1406 109920 : exc = pw_integrate_function(v_drho_r)
1407 : !
1408 : ! return the xc functional value at the grid points
1409 : !
1410 109920 : IF (PRESENT(exc_r)) THEN
1411 96 : exc_r = v_drho_r
1412 : ELSE
1413 109824 : CALL v_drho_r%release()
1414 : END IF
1415 : ELSE
1416 320 : exc = 0.0_dp
1417 : END IF
1418 :
1419 110240 : CALL xc_rho_set_release(rho_set, pw_pool=pw_pool)
1420 :
1421 : ! tau part
1422 110240 : IF (has_tau) THEN
1423 7968 : ALLOCATE (vxc_tau(nspins))
1424 2454 : IF (lsd) THEN
1425 606 : CALL xc_dset_recover_pw(deriv_set, [deriv_tau_a], vxc_tau(1), pw_grid)
1426 606 : CALL xc_dset_recover_pw(deriv_set, [deriv_tau_b], vxc_tau(2), pw_grid)
1427 : ELSE
1428 1848 : CALL xc_dset_recover_pw(deriv_set, [deriv_tau], vxc_tau(1), pw_grid)
1429 : END IF
1430 5514 : DO ispin = 1, nspins
1431 5514 : CPASSERT(ASSOCIATED(vxc_tau(ispin)%array))
1432 : END DO
1433 : END IF
1434 110240 : CALL xc_dset_release(deriv_set)
1435 :
1436 110240 : CALL timestop(handle)
1437 :
1438 2315040 : END SUBROUTINE xc_vxc_pw_create
1439 :
1440 : ! **************************************************************************************************
1441 : !> \brief calculates just the exchange and correlation energy
1442 : !> (no vxc)
1443 : !> \param rho_r realspace density on the grid
1444 : !> \param rho_g g-space density on the grid
1445 : !> \param tau kinetic energy density on the grid
1446 : !> \param xc_section XC parameters
1447 : !> \param pw_pool pool of plain-wave grids
1448 : !> \return the XC energy
1449 : !> \par History
1450 : !> 11.2003 created [fawzi]
1451 : !> \author fawzi
1452 : !> \note
1453 : !> has to be kept consistent with xc_vxc_pw_create
1454 : ! **************************************************************************************************
1455 15638 : FUNCTION xc_exc_calc(rho_r, rho_g, tau, xc_section, pw_pool) &
1456 : RESULT(exc)
1457 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, tau
1458 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
1459 : TYPE(section_vals_type), POINTER :: xc_section
1460 : TYPE(pw_pool_type), POINTER :: pw_pool
1461 : REAL(kind=dp) :: exc
1462 :
1463 : CHARACTER(len=*), PARAMETER :: routineN = 'xc_exc_calc'
1464 :
1465 : INTEGER :: handle
1466 : REAL(dp) :: density_smooth_cut_range, rho_cutoff
1467 7819 : REAL(dp), DIMENSION(:, :, :), POINTER :: e_0
1468 : TYPE(xc_derivative_set_type) :: deriv_set
1469 : TYPE(xc_derivative_type), POINTER :: deriv
1470 : TYPE(xc_rho_set_type) :: rho_set
1471 :
1472 7819 : CALL timeset(routineN, handle)
1473 7819 : NULLIFY (deriv, e_0)
1474 7819 : exc = 0.0_dp
1475 :
1476 : ! this has to be consistent with what is done in xc_vxc_pw_create
1477 : CALL xc_rho_set_and_dset_create(rho_set=rho_set, &
1478 : deriv_set=deriv_set, deriv_order=0, &
1479 : rho_r=rho_r, rho_g=rho_g, tau=tau, xc_section=xc_section, &
1480 : pw_pool=pw_pool, &
1481 7819 : calc_potential=.FALSE.)
1482 7819 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::])
1483 :
1484 7819 : IF (ASSOCIATED(deriv)) THEN
1485 7819 : CALL xc_derivative_get(deriv, deriv_data=e_0)
1486 :
1487 : CALL section_vals_val_get(xc_section, "DENSITY_CUTOFF", &
1488 7819 : r_val=rho_cutoff)
1489 : CALL section_vals_val_get(xc_section, "DENSITY_SMOOTH_CUTOFF_RANGE", &
1490 7819 : r_val=density_smooth_cut_range)
1491 : CALL smooth_cutoff(pot=e_0, rho=rho_set%rho, &
1492 : rhoa=rho_set%rhoa, rhob=rho_set%rhob, &
1493 : rho_cutoff=rho_cutoff, &
1494 7819 : rho_smooth_cutoff_range=density_smooth_cut_range)
1495 :
1496 7819 : exc = accurate_sum(e_0)*rho_r(1)%pw_grid%dvol
1497 7819 : IF (rho_r(1)%pw_grid%para%mode == PW_MODE_DISTRIBUTED) THEN
1498 7694 : CALL rho_r(1)%pw_grid%para%group%sum(exc)
1499 : END IF
1500 :
1501 7819 : CALL xc_rho_set_release(rho_set, pw_pool=pw_pool)
1502 7819 : CALL xc_dset_release(deriv_set)
1503 : END IF
1504 7819 : CALL timestop(handle)
1505 148561 : END FUNCTION xc_exc_calc
1506 :
1507 : ! **************************************************************************************************
1508 : !> \brief Caller routine to calculate the second order potential in the direction of rho1_r
1509 : !> \param v_xc XC potential, will be allocated, to be integrated with the KS density
1510 : !> \param v_xc_tau ...
1511 : !> \param deriv_set XC derivatives from xc_prep_2nd_deriv
1512 : !> \param rho_set XC rho set from KS rho from xc_prep_2nd_deriv
1513 : !> \param rho1_r first-order density in r space
1514 : !> \param rho1_g first-order density in g space
1515 : !> \param tau1_r ...
1516 : !> \param pw_pool pw pool to create new grids
1517 : !> \param xc_section XC section to calculate the derivatives from
1518 : !> \param gapw whether to carry out GAPW (not possible with numerical derivatives)
1519 : !> \param vxg GAPW potential
1520 : !> \param lsd_singlets ...
1521 : !> \param do_excitations ...
1522 : !> \param do_triplet ...
1523 : !> \param do_tddft ...
1524 : !> \param compute_virial ...
1525 : !> \param virial_xc virial terms will be collected here
1526 : ! **************************************************************************************************
1527 14310 : SUBROUTINE xc_calc_2nd_deriv(v_xc, v_xc_tau, deriv_set, rho_set, rho1_r, rho1_g, tau1_r, &
1528 : pw_pool, xc_section, gapw, vxg, &
1529 : lsd_singlets, do_excitations, do_triplet, do_tddft, &
1530 : compute_virial, virial_xc)
1531 :
1532 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_xc, v_xc_tau
1533 : TYPE(xc_derivative_set_type) :: deriv_set
1534 : TYPE(xc_rho_set_type) :: rho_set
1535 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho1_r, tau1_r
1536 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho1_g
1537 : TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool
1538 : TYPE(section_vals_type), INTENT(IN), POINTER :: xc_section
1539 : LOGICAL, INTENT(IN) :: gapw
1540 : REAL(KIND=dp), DIMENSION(:, :, :, :), OPTIONAL, &
1541 : POINTER :: vxg
1542 : LOGICAL, INTENT(IN), OPTIONAL :: lsd_singlets, do_excitations, &
1543 : do_triplet, do_tddft, compute_virial
1544 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT), &
1545 : OPTIONAL :: virial_xc
1546 :
1547 : CHARACTER(len=*), PARAMETER :: routineN = 'xc_calc_2nd_deriv'
1548 :
1549 : INTEGER :: handle, ispin, nspins
1550 : INTEGER, DIMENSION(2, 3) :: bo
1551 : LOGICAL :: lsd, my_compute_virial, &
1552 : my_do_excitations, my_do_tddft, &
1553 : my_do_triplet, my_lsd_singlets
1554 : REAL(KIND=dp) :: fac
1555 : TYPE(section_vals_type), POINTER :: xc_fun_section
1556 : TYPE(xc_rho_cflags_type) :: needs
1557 : TYPE(xc_rho_set_type) :: rho1_set
1558 :
1559 14310 : CALL timeset(routineN, handle)
1560 :
1561 14310 : my_compute_virial = .FALSE.
1562 14310 : IF (PRESENT(compute_virial)) my_compute_virial = compute_virial
1563 :
1564 14310 : my_do_tddft = .FALSE.
1565 14310 : IF (PRESENT(do_tddft)) my_do_tddft = do_tddft
1566 :
1567 14310 : my_do_excitations = .FALSE.
1568 14310 : IF (PRESENT(do_excitations)) my_do_excitations = do_excitations
1569 :
1570 14310 : my_lsd_singlets = .FALSE.
1571 14310 : IF (PRESENT(lsd_singlets)) my_lsd_singlets = lsd_singlets
1572 :
1573 14310 : my_do_triplet = .FALSE.
1574 14310 : IF (PRESENT(do_triplet)) my_do_triplet = do_triplet
1575 :
1576 14310 : nspins = SIZE(rho1_r)
1577 14310 : lsd = (nspins == 2)
1578 14310 : IF (nspins == 1 .AND. my_do_excitations .AND. (my_lsd_singlets .OR. my_do_triplet)) THEN
1579 0 : nspins = 2
1580 0 : lsd = .TRUE.
1581 : END IF
1582 :
1583 14310 : NULLIFY (v_xc, v_xc_tau)
1584 59264 : ALLOCATE (v_xc(nspins))
1585 30644 : DO ispin = 1, nspins
1586 16334 : CALL pw_pool%create_pw(v_xc(ispin))
1587 30644 : CALL pw_zero(v_xc(ispin))
1588 : END DO
1589 :
1590 14310 : xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
1591 14310 : needs = xc_functionals_get_needs(xc_fun_section, lsd, .TRUE.)
1592 :
1593 14310 : IF (needs%tau .OR. needs%tau_spin) THEN
1594 536 : IF (.NOT. ASSOCIATED(tau1_r)) &
1595 0 : CPABORT("Tau-dependent functionals requires allocated kinetic energy density grid")
1596 1736 : ALLOCATE (v_xc_tau(nspins))
1597 1200 : DO ispin = 1, nspins
1598 664 : CALL pw_pool%create_pw(v_xc_tau(ispin))
1599 14974 : CALL pw_zero(v_xc_tau(ispin))
1600 : END DO
1601 : END IF
1602 :
1603 14310 : IF (section_get_lval(xc_section, "2ND_DERIV_ANALYTICAL") .AND. .NOT. my_do_tddft) THEN
1604 : !------!
1605 : ! rho1 !
1606 : !------!
1607 135300 : bo = rho1_r(1)%pw_grid%bounds_local
1608 : ! create the place where to store the argument for the functionals
1609 : CALL xc_rho_set_create(rho1_set, bo, &
1610 : rho_cutoff=section_get_rval(xc_section, "DENSITY_CUTOFF"), &
1611 : drho_cutoff=section_get_rval(xc_section, "GRADIENT_CUTOFF"), &
1612 13530 : tau_cutoff=section_get_rval(xc_section, "TAU_CUTOFF"))
1613 :
1614 : ! calculate the arguments needed by the functionals
1615 : CALL xc_rho_set_update(rho1_set, rho1_r, rho1_g, tau1_r, needs, &
1616 : section_get_ival(xc_section, "XC_GRID%XC_DERIV"), &
1617 : section_get_ival(xc_section, "XC_GRID%XC_SMOOTH_RHO"), &
1618 13530 : pw_pool)
1619 :
1620 13530 : fac = 0._dp
1621 13530 : IF (nspins == 1 .AND. my_do_excitations) THEN
1622 1074 : IF (my_lsd_singlets) fac = 1.0_dp
1623 1074 : IF (my_do_triplet) fac = -1.0_dp
1624 : END IF
1625 :
1626 : CALL xc_calc_2nd_deriv_analytical(v_xc, v_xc_tau, deriv_set, rho_set, &
1627 : rho1_set, pw_pool, xc_section, gapw, vxg=vxg, &
1628 13530 : tddfpt_fac=fac, compute_virial=compute_virial, virial_xc=virial_xc)
1629 :
1630 13530 : CALL xc_rho_set_release(rho1_set)
1631 :
1632 : ELSE
1633 780 : IF (gapw) CPABORT("Numerical 2nd derivatives not implemented with GAPW")
1634 :
1635 : CALL xc_calc_2nd_deriv_numerical(v_xc, v_xc_tau, rho_set, rho1_r, rho1_g, tau1_r, &
1636 : pw_pool, xc_section, &
1637 : my_do_excitations .AND. my_do_triplet, &
1638 780 : compute_virial, virial_xc, deriv_set)
1639 : END IF
1640 :
1641 14310 : CALL timestop(handle)
1642 :
1643 314820 : END SUBROUTINE xc_calc_2nd_deriv
1644 :
1645 : ! **************************************************************************************************
1646 : !> \brief calculates 2nd derivative numerically
1647 : !> \param v_xc potential to be calculated (has to be allocated already)
1648 : !> \param v_tau tau-part of the potential to be calculated (has to be allocated already)
1649 : !> \param rho_set KS density from xc_prep_2nd_deriv
1650 : !> \param rho1_r first-order density in r-space
1651 : !> \param rho1_g first-order density in g-space
1652 : !> \param tau1_r first-order kinetic-energy density in r-space
1653 : !> \param pw_pool pw pool for new grids
1654 : !> \param xc_section XC section to calculate the derivatives from
1655 : !> \param do_triplet ...
1656 : !> \param calc_virial whether to calculate virial terms
1657 : !> \param virial_xc collects stress tensor components (no metaGGAs!)
1658 : !> \param deriv_set deriv set from xc_prep_2nd_deriv (only for virials)
1659 : ! **************************************************************************************************
1660 832 : SUBROUTINE xc_calc_2nd_deriv_numerical(v_xc, v_tau, rho_set, rho1_r, rho1_g, tau1_r, &
1661 : pw_pool, xc_section, &
1662 : do_triplet, calc_virial, virial_xc, deriv_set)
1663 :
1664 : TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(IN), POINTER :: v_xc, v_tau
1665 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
1666 : TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(IN), POINTER :: rho1_r, tau1_r
1667 : TYPE(pw_c1d_gs_type), DIMENSION(:), INTENT(IN), POINTER :: rho1_g
1668 : TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool
1669 : TYPE(section_vals_type), INTENT(IN), POINTER :: xc_section
1670 : LOGICAL, INTENT(IN) :: do_triplet
1671 : LOGICAL, INTENT(IN), OPTIONAL :: calc_virial
1672 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT), &
1673 : OPTIONAL :: virial_xc
1674 : TYPE(xc_derivative_set_type), OPTIONAL :: deriv_set
1675 :
1676 : CHARACTER(len=*), PARAMETER :: routineN = 'xc_calc_2nd_deriv_numerical'
1677 : REAL(KIND=dp), DIMENSION(-4:4, 4), PARAMETER :: &
1678 : weights = RESHAPE([0.0_dp, 0.0_dp, 0.0_dp, -0.5_dp, 0.0_dp, 0.5_dp, 0.0_dp, 0.0_dp, 0.0_dp, &
1679 : 0.0_dp, 0.0_dp, 1.0_dp/12.0_dp, -2.0_dp/3.0_dp, 0.0_dp, 2.0_dp/3.0_dp, -1.0_dp/12.0_dp, 0.0_dp, 0.0_dp, &
1680 : 0.0_dp, -1.0_dp/60.0_dp, 0.15_dp, -0.75_dp, 0.0_dp, 0.75_dp, -0.15_dp, 1.0_dp/60.0_dp, 0.0_dp, &
1681 : 1.0_dp/280.0_dp, -4.0_dp/105.0_dp, 0.2_dp, -0.8_dp, 0.0_dp, 0.8_dp, -0.2_dp, 4.0_dp/105.0_dp, -1.0_dp/280.0_dp], [9, 4])
1682 :
1683 : INTEGER :: handle, idir, ispin, nspins, istep, nsteps
1684 : INTEGER, DIMENSION(2, 3) :: bo
1685 : LOGICAL :: gradient_f, lsd, my_calc_virial, tau_f, laplace_f, rho_f
1686 : REAL(KIND=dp) :: exc, gradient_cut, h, weight, step, rho_cutoff
1687 832 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: dr1dr, dra1dra, drb1drb
1688 : REAL(KIND=dp), DIMENSION(3, 3) :: virial_dummy
1689 832 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: norm_drho, norm_drho2, norm_drho2a, &
1690 832 : norm_drho2b, norm_drhoa, norm_drhob, &
1691 2496 : rho, rho1, rho1a, rho1b, rhoa, rhob, &
1692 1664 : tau_a, tau_b, tau, tau1, tau1a, tau1b, laplace, laplace1, &
1693 832 : laplacea, laplaceb, laplace1a, laplace1b, &
1694 1664 : laplace2, laplace2a, laplace2b, deriv_data
1695 19968 : TYPE(cp_3d_r_cp_type), DIMENSION(3) :: drho, drho1, drho1a, drho1b, drhoa, drhob
1696 : TYPE(pw_r3d_rs_type) :: v_drho, v_drhoa, v_drhob
1697 832 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: vxc_rho, vxc_tau
1698 832 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
1699 832 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, tau_r
1700 : TYPE(pw_r3d_rs_type) :: virial_pw, v_laplace, v_laplacea, v_laplaceb
1701 : TYPE(section_vals_type), POINTER :: xc_fun_section
1702 : TYPE(xc_derivative_set_type) :: deriv_set1
1703 : TYPE(xc_rho_cflags_type) :: needs
1704 : TYPE(xc_rho_set_type) :: rho1_set, rho2_set
1705 :
1706 832 : CALL timeset(routineN, handle)
1707 :
1708 832 : my_calc_virial = .FALSE.
1709 832 : IF (PRESENT(calc_virial) .AND. PRESENT(virial_xc)) my_calc_virial = calc_virial
1710 :
1711 832 : nspins = SIZE(v_xc)
1712 :
1713 832 : NULLIFY (tau, tau_r, tau_a, tau_b)
1714 :
1715 832 : h = section_get_rval(xc_section, "STEP_SIZE")
1716 832 : nsteps = section_get_ival(xc_section, "NSTEPS")
1717 832 : IF (nsteps < LBOUND(weights, 2) .OR. nspins > UBOUND(weights, 2)) THEN
1718 0 : CPABORT("The number of steps must be a value from 1 to 4.")
1719 : END IF
1720 :
1721 832 : IF (nspins == 2) THEN
1722 324 : NULLIFY (vxc_rho, rho_g, vxc_tau)
1723 972 : ALLOCATE (rho_r(2))
1724 972 : DO ispin = 1, nspins
1725 972 : CALL pw_pool%create_pw(rho_r(ispin))
1726 : END DO
1727 324 : IF (ASSOCIATED(tau1_r) .AND. ASSOCIATED(v_tau)) THEN
1728 162 : ALLOCATE (tau_r(2))
1729 162 : DO ispin = 1, nspins
1730 162 : CALL pw_pool%create_pw(tau_r(ispin))
1731 : END DO
1732 : END IF
1733 324 : CALL xc_rho_set_get(rho_set, can_return_null=.TRUE., rhoa=rhoa, rhob=rhob, tau_a=tau_a, tau_b=tau_b)
1734 2496 : DO istep = -nsteps, nsteps
1735 2172 : IF (istep == 0) CYCLE
1736 1848 : weight = weights(istep, nsteps)/h
1737 1848 : step = REAL(istep, dp)*h
1738 : CALL calc_resp_potential_numer_ab(rho_r, rho_g, rho1_r, rhoa, rhob, vxc_rho, &
1739 1848 : tau_r, tau1_r, tau_a, tau_b, vxc_tau, xc_section, pw_pool, step)
1740 5544 : DO ispin = 1, nspins
1741 3696 : CALL pw_axpy(vxc_rho(ispin), v_xc(ispin), weight)
1742 5544 : IF (ASSOCIATED(vxc_tau) .AND. ASSOCIATED(v_tau)) THEN
1743 456 : CALL pw_axpy(vxc_tau(ispin), v_tau(ispin), weight)
1744 : END IF
1745 : END DO
1746 5544 : DO ispin = 1, nspins
1747 5544 : CALL vxc_rho(ispin)%release()
1748 : END DO
1749 1848 : DEALLOCATE (vxc_rho)
1750 2172 : IF (ASSOCIATED(vxc_tau)) THEN
1751 684 : DO ispin = 1, nspins
1752 684 : CALL vxc_tau(ispin)%release()
1753 : END DO
1754 228 : DEALLOCATE (vxc_tau)
1755 : END IF
1756 : END DO
1757 508 : ELSE IF (nspins == 1 .AND. do_triplet) THEN
1758 24 : NULLIFY (vxc_rho, vxc_tau, rho_g)
1759 72 : ALLOCATE (rho_r(2))
1760 72 : DO ispin = 1, 2
1761 72 : CALL pw_pool%create_pw(rho_r(ispin))
1762 : END DO
1763 24 : IF (ASSOCIATED(tau1_r) .AND. ASSOCIATED(v_tau)) THEN
1764 0 : ALLOCATE (tau_r(2))
1765 0 : DO ispin = 1, nspins
1766 0 : CALL pw_pool%create_pw(tau_r(ispin))
1767 : END DO
1768 : END IF
1769 24 : CALL xc_rho_set_get(rho_set, can_return_null=.TRUE., rhoa=rhoa, rhob=rhob, tau_a=tau_a, tau_b=tau_b)
1770 192 : DO istep = -nsteps, nsteps
1771 168 : IF (istep == 0) CYCLE
1772 144 : weight = weights(istep, nsteps)/h
1773 144 : step = REAL(istep, dp)*h
1774 : ! K(alpha,alpha)
1775 144 : !$OMP PARALLEL DEFAULT(NONE) SHARED(rho_r,rhoa,rhob,step,rho1_r,tau_r,tau_a,tau_b,tau1_r)
1776 : !$OMP WORKSHARE
1777 : rho_r(1)%array(:, :, :) = rhoa(:, :, :) + step*rho1_r(1)%array(:, :, :)
1778 : !$OMP END WORKSHARE NOWAIT
1779 : !$OMP WORKSHARE
1780 : rho_r(2)%array(:, :, :) = rhob(:, :, :)
1781 : !$OMP END WORKSHARE NOWAIT
1782 : IF (ASSOCIATED(tau1_r)) THEN
1783 : !$OMP WORKSHARE
1784 : tau_r(1)%array(:, :, :) = tau_a(:, :, :) + step*tau1_r(1)%array(:, :, :)
1785 : !$OMP END WORKSHARE NOWAIT
1786 : !$OMP WORKSHARE
1787 : tau_r(2)%array(:, :, :) = tau_b(:, :, :)
1788 : !$OMP END WORKSHARE NOWAIT
1789 : END IF
1790 : !$OMP END PARALLEL
1791 : CALL xc_vxc_pw_create(vxc_rho, vxc_tau, exc, rho_r, rho_g, tau_r, xc_section, &
1792 144 : pw_pool, .FALSE., virial_dummy)
1793 144 : CALL pw_axpy(vxc_rho(1), v_xc(1), weight)
1794 144 : IF (ASSOCIATED(vxc_tau) .AND. ASSOCIATED(v_tau)) THEN
1795 0 : CALL pw_axpy(vxc_tau(1), v_tau(1), weight)
1796 : END IF
1797 432 : DO ispin = 1, 2
1798 432 : CALL vxc_rho(ispin)%release()
1799 : END DO
1800 144 : DEALLOCATE (vxc_rho)
1801 144 : IF (ASSOCIATED(vxc_tau)) THEN
1802 0 : DO ispin = 1, 2
1803 0 : CALL vxc_tau(ispin)%release()
1804 : END DO
1805 0 : DEALLOCATE (vxc_tau)
1806 : END IF
1807 144 : !$OMP PARALLEL DEFAULT(NONE) SHARED(rho_r,rhoa,rhob,step,rho1_r,tau_r,tau_a,tau_b,tau1_r)
1808 : !$OMP WORKSHARE
1809 : ! K(alpha,beta)
1810 : rho_r(1)%array(:, :, :) = rhoa(:, :, :)
1811 : !$OMP END WORKSHARE NOWAIT
1812 : !$OMP WORKSHARE
1813 : rho_r(2)%array(:, :, :) = rhob(:, :, :) + step*rho1_r(1)%array(:, :, :)
1814 : !$OMP END WORKSHARE NOWAIT
1815 : IF (ASSOCIATED(tau1_r)) THEN
1816 : !$OMP WORKSHARE
1817 : tau_r(1)%array(:, :, :) = tau_a(:, :, :)
1818 : !$OMP END WORKSHARE NOWAIT
1819 : !$OMP WORKSHARE
1820 : tau_r(2)%array(:, :, :) = tau_b(:, :, :) + step*tau1_r(1)%array(:, :, :)
1821 : !$OMP END WORKSHARE NOWAIT
1822 : END IF
1823 : !$OMP END PARALLEL
1824 : CALL xc_vxc_pw_create(vxc_rho, vxc_tau, exc, rho_r, rho_g, tau_r, xc_section, &
1825 144 : pw_pool, .FALSE., virial_dummy)
1826 144 : CALL pw_axpy(vxc_rho(1), v_xc(1), weight)
1827 144 : IF (ASSOCIATED(vxc_tau) .AND. ASSOCIATED(v_tau)) THEN
1828 0 : CALL pw_axpy(vxc_tau(1), v_tau(1), weight)
1829 : END IF
1830 432 : DO ispin = 1, 2
1831 432 : CALL vxc_rho(ispin)%release()
1832 : END DO
1833 144 : DEALLOCATE (vxc_rho)
1834 168 : IF (ASSOCIATED(vxc_tau)) THEN
1835 0 : DO ispin = 1, 2
1836 0 : CALL vxc_tau(ispin)%release()
1837 : END DO
1838 0 : DEALLOCATE (vxc_tau)
1839 : END IF
1840 : END DO
1841 : ELSE
1842 484 : NULLIFY (vxc_rho, rho_r, rho_g, vxc_tau, tau_r, tau)
1843 968 : ALLOCATE (rho_r(1))
1844 484 : CALL pw_pool%create_pw(rho_r(1))
1845 484 : IF (ASSOCIATED(tau1_r) .AND. ASSOCIATED(v_tau)) THEN
1846 92 : ALLOCATE (tau_r(1))
1847 46 : CALL pw_pool%create_pw(tau_r(1))
1848 : END IF
1849 484 : CALL xc_rho_set_get(rho_set, can_return_null=.TRUE., rho=rho, tau=tau)
1850 3776 : DO istep = -nsteps, nsteps
1851 3292 : IF (istep == 0) CYCLE
1852 2808 : weight = weights(istep, nsteps)/h
1853 2808 : step = REAL(istep, dp)*h
1854 2808 : !$OMP PARALLEL DEFAULT(NONE) SHARED(rho_r,rho,step,rho1_r,tau1_r,tau,tau_r)
1855 : !$OMP WORKSHARE
1856 : rho_r(1)%array(:, :, :) = rho(:, :, :) + step*rho1_r(1)%array(:, :, :)
1857 : !$OMP END WORKSHARE NOWAIT
1858 : IF (ASSOCIATED(tau1_r) .AND. ASSOCIATED(tau) .AND. ASSOCIATED(tau1_r)) THEN
1859 : !$OMP WORKSHARE
1860 : tau_r(1)%array(:, :, :) = tau(:, :, :) + step*tau1_r(1)%array(:, :, :)
1861 : !$OMP END WORKSHARE NOWAIT
1862 : END IF
1863 : !$OMP END PARALLEL
1864 : CALL xc_vxc_pw_create(vxc_rho, vxc_tau, exc, rho_r, rho_g, tau_r, xc_section, &
1865 2808 : pw_pool, .FALSE., virial_dummy)
1866 2808 : CALL pw_axpy(vxc_rho(1), v_xc(1), weight)
1867 2808 : IF (ASSOCIATED(vxc_tau) .AND. ASSOCIATED(v_tau)) THEN
1868 276 : CALL pw_axpy(vxc_tau(1), v_tau(1), weight)
1869 : END IF
1870 2808 : CALL vxc_rho(1)%release()
1871 2808 : DEALLOCATE (vxc_rho)
1872 3292 : IF (ASSOCIATED(vxc_tau)) THEN
1873 276 : CALL vxc_tau(1)%release()
1874 276 : DEALLOCATE (vxc_tau)
1875 : END IF
1876 : END DO
1877 : END IF
1878 :
1879 832 : IF (my_calc_virial) THEN
1880 36 : lsd = (nspins == 2)
1881 36 : IF (nspins == 1 .AND. do_triplet) THEN
1882 0 : lsd = .TRUE.
1883 : END IF
1884 :
1885 36 : CALL check_for_derivatives(deriv_set, (nspins == 2), rho_f, gradient_f, tau_f, laplace_f)
1886 :
1887 : ! Calculate the virial terms
1888 : ! Those arising from the first derivatives are treated like in xc_calc_2nd_deriv_analytical
1889 : ! Those arising from the second derivatives are calculated numerically
1890 : ! We assume that all metaGGA functionals require the gradient
1891 36 : IF (gradient_f) THEN
1892 360 : bo = rho_set%local_bounds
1893 :
1894 : ! Create the work grid for the virial terms
1895 36 : CALL allocate_pw(virial_pw, pw_pool, bo)
1896 :
1897 36 : gradient_cut = section_get_rval(xc_section, "GRADIENT_CUTOFF")
1898 :
1899 : ! create the container to store the argument of the functionals
1900 : CALL xc_rho_set_create(rho1_set, bo, &
1901 : rho_cutoff=section_get_rval(xc_section, "DENSITY_CUTOFF"), &
1902 : drho_cutoff=gradient_cut, &
1903 36 : tau_cutoff=section_get_rval(xc_section, "TAU_CUTOFF"))
1904 :
1905 36 : xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
1906 36 : needs = xc_functionals_get_needs(xc_fun_section, lsd, .TRUE.)
1907 :
1908 : ! calculate the arguments needed by the functionals
1909 : CALL xc_rho_set_update(rho1_set, rho1_r, rho1_g, tau1_r, needs, &
1910 : section_get_ival(xc_section, "XC_GRID%XC_DERIV"), &
1911 : section_get_ival(xc_section, "XC_GRID%XC_SMOOTH_RHO"), &
1912 36 : pw_pool)
1913 :
1914 36 : IF (lsd) THEN
1915 : CALL xc_rho_set_get(rho_set, drhoa=drhoa, drhob=drhob, norm_drho=norm_drho, &
1916 : norm_drhoa=norm_drhoa, norm_drhob=norm_drhob, tau_a=tau_a, tau_b=tau_b, &
1917 10 : laplace_rhoa=laplacea, laplace_rhob=laplaceb, can_return_null=.TRUE.)
1918 : CALL xc_rho_set_get(rho1_set, rhoa=rho1a, rhob=rho1b, drhoa=drho1a, drhob=drho1b, laplace_rhoa=laplace1a, &
1919 10 : laplace_rhob=laplace1b, can_return_null=.TRUE.)
1920 :
1921 10 : CALL calc_drho_from_ab(drho, drhoa, drhob)
1922 10 : CALL calc_drho_from_ab(drho1, drho1a, drho1b)
1923 : ELSE
1924 26 : CALL xc_rho_set_get(rho_set, drho=drho, norm_drho=norm_drho, tau=tau, laplace_rho=laplace, can_return_null=.TRUE.)
1925 26 : CALL xc_rho_set_get(rho1_set, rho=rho1, drho=drho1, laplace_rho=laplace1, can_return_null=.TRUE.)
1926 : END IF
1927 :
1928 36 : CALL prepare_dr1dr(dr1dr, drho, drho1)
1929 :
1930 36 : IF (lsd) THEN
1931 10 : CALL prepare_dr1dr(dra1dra, drhoa, drho1a)
1932 10 : CALL prepare_dr1dr(drb1drb, drhob, drho1b)
1933 :
1934 10 : CALL allocate_pw(v_drho, pw_pool, bo)
1935 10 : CALL allocate_pw(v_drhoa, pw_pool, bo)
1936 10 : CALL allocate_pw(v_drhob, pw_pool, bo)
1937 :
1938 10 : IF (ASSOCIATED(norm_drhoa)) CALL apply_drho(deriv_set, [deriv_norm_drhoa], virial_pw, drhoa, drho1a, virial_xc, &
1939 10 : norm_drhoa, gradient_cut, dra1dra, v_drhoa%array)
1940 10 : IF (ASSOCIATED(norm_drhob)) CALL apply_drho(deriv_set, [deriv_norm_drhob], virial_pw, drhob, drho1b, virial_xc, &
1941 10 : norm_drhob, gradient_cut, drb1drb, v_drhob%array)
1942 10 : IF (ASSOCIATED(norm_drho)) CALL apply_drho(deriv_set, [deriv_norm_drho], virial_pw, drho, drho1, virial_xc, &
1943 6 : norm_drho, gradient_cut, dr1dr, v_drho%array)
1944 10 : IF (laplace_f) THEN
1945 2 : CALL xc_derivative_get(xc_dset_get_derivative(deriv_set, [deriv_laplace_rhoa]), deriv_data=deriv_data)
1946 2 : CPASSERT(ASSOCIATED(deriv_data))
1947 15026 : virial_pw%array(:, :, :) = -rho1a(:, :, :)
1948 2 : CALL virial_laplace(virial_pw, pw_pool, virial_xc, deriv_data)
1949 :
1950 2 : CALL allocate_pw(v_laplacea, pw_pool, bo)
1951 :
1952 2 : CALL xc_derivative_get(xc_dset_get_derivative(deriv_set, [deriv_laplace_rhob]), deriv_data=deriv_data)
1953 2 : CPASSERT(ASSOCIATED(deriv_data))
1954 15026 : virial_pw%array(:, :, :) = -rho1b(:, :, :)
1955 2 : CALL virial_laplace(virial_pw, pw_pool, virial_xc, deriv_data)
1956 :
1957 2 : CALL allocate_pw(v_laplaceb, pw_pool, bo)
1958 : END IF
1959 :
1960 : ELSE
1961 :
1962 : ! Create the work grid for the potential of the gradient part
1963 26 : CALL allocate_pw(v_drho, pw_pool, bo)
1964 :
1965 : CALL apply_drho(deriv_set, [deriv_norm_drho], virial_pw, drho, drho1, virial_xc, &
1966 26 : norm_drho, gradient_cut, dr1dr, v_drho%array)
1967 26 : IF (laplace_f) THEN
1968 2 : CALL xc_derivative_get(xc_dset_get_derivative(deriv_set, [deriv_laplace_rho]), deriv_data=deriv_data)
1969 2 : CPASSERT(ASSOCIATED(deriv_data))
1970 28862 : virial_pw%array(:, :, :) = -rho1(:, :, :)
1971 2 : CALL virial_laplace(virial_pw, pw_pool, virial_xc, deriv_data)
1972 :
1973 2 : CALL allocate_pw(v_laplace, pw_pool, bo)
1974 : END IF
1975 :
1976 : END IF
1977 :
1978 36 : IF (lsd) THEN
1979 150260 : rho_r(1)%array = rhoa
1980 150260 : rho_r(2)%array = rhob
1981 : ELSE
1982 701552 : rho_r(1)%array = rho
1983 : END IF
1984 36 : IF (ASSOCIATED(tau1_r)) THEN
1985 8 : IF (lsd) THEN
1986 60104 : tau_r(1)%array = tau_a
1987 60104 : tau_r(2)%array = tau_b
1988 : ELSE
1989 115448 : tau_r(1)%array = tau
1990 : END IF
1991 : END IF
1992 :
1993 : ! Create deriv sets with same densities but different gradients
1994 36 : CALL xc_dset_create(deriv_set1, pw_pool)
1995 :
1996 36 : rho_cutoff = section_get_rval(xc_section, "DENSITY_CUTOFF")
1997 :
1998 : ! create the place where to store the argument for the functionals
1999 : CALL xc_rho_set_create(rho2_set, bo, &
2000 : rho_cutoff=rho_cutoff, &
2001 : drho_cutoff=section_get_rval(xc_section, "GRADIENT_CUTOFF"), &
2002 36 : tau_cutoff=section_get_rval(xc_section, "TAU_CUTOFF"))
2003 :
2004 : ! calculate the arguments needed by the functionals
2005 : CALL xc_rho_set_update(rho2_set, rho_r, rho_g, tau_r, needs, &
2006 : section_get_ival(xc_section, "XC_GRID%XC_DERIV"), &
2007 : section_get_ival(xc_section, "XC_GRID%XC_SMOOTH_RHO"), &
2008 36 : pw_pool)
2009 :
2010 36 : IF (lsd) THEN
2011 : CALL xc_rho_set_get(rho1_set, rhoa=rho1a, rhob=rho1b, tau_a=tau1a, tau_b=tau1b, &
2012 10 : laplace_rhoa=laplace1a, laplace_rhob=laplace1b, can_return_null=.TRUE.)
2013 : CALL xc_rho_set_get(rho2_set, norm_drhoa=norm_drho2a, norm_drhob=norm_drho2b, &
2014 10 : norm_drho=norm_drho2, laplace_rhoa=laplace2a, laplace_rhob=laplace2b, can_return_null=.TRUE.)
2015 :
2016 64 : DO istep = -nsteps, nsteps
2017 54 : IF (istep == 0) CYCLE
2018 44 : weight = weights(istep, nsteps)/h
2019 44 : step = REAL(istep, dp)*h
2020 44 : IF (ASSOCIATED(norm_drhoa)) THEN
2021 44 : CALL get_derivs_rho(norm_drho2a, norm_drhoa, step, xc_fun_section, lsd, rho2_set, deriv_set1)
2022 : CALL update_deriv_rho(deriv_set1, [deriv_rhoa], bo, &
2023 44 : norm_drhoa, gradient_cut, weight, rho1a, v_drhoa%array)
2024 : CALL update_deriv_rho(deriv_set1, [deriv_rhob], bo, &
2025 44 : norm_drhoa, gradient_cut, weight, rho1b, v_drhoa%array)
2026 : CALL update_deriv_rho(deriv_set1, [deriv_norm_drhoa], bo, &
2027 44 : norm_drhoa, gradient_cut, weight, dra1dra, v_drhoa%array)
2028 : CALL update_deriv_drho_ab(deriv_set1, [deriv_norm_drhob], bo, &
2029 44 : norm_drhoa, gradient_cut, weight, dra1dra, drb1drb, v_drhoa%array, v_drhob%array)
2030 : CALL update_deriv_drho_ab(deriv_set1, [deriv_norm_drho], bo, &
2031 44 : norm_drhoa, gradient_cut, weight, dra1dra, dr1dr, v_drhoa%array, v_drho%array)
2032 44 : IF (tau_f) THEN
2033 : CALL update_deriv_rho(deriv_set1, [deriv_tau_a], bo, &
2034 8 : norm_drhoa, gradient_cut, weight, tau1a, v_drhoa%array)
2035 : CALL update_deriv_rho(deriv_set1, [deriv_tau_b], bo, &
2036 8 : norm_drhoa, gradient_cut, weight, tau1b, v_drhoa%array)
2037 : END IF
2038 44 : IF (laplace_f) THEN
2039 : CALL update_deriv_rho(deriv_set1, [deriv_laplace_rhoa], bo, &
2040 4 : norm_drhoa, gradient_cut, weight, laplace1a, v_drhoa%array)
2041 : CALL update_deriv_rho(deriv_set1, [deriv_laplace_rhob], bo, &
2042 4 : norm_drhoa, gradient_cut, weight, laplace1b, v_drhoa%array)
2043 : END IF
2044 : END IF
2045 :
2046 44 : IF (ASSOCIATED(norm_drhob)) THEN
2047 44 : CALL get_derivs_rho(norm_drho2b, norm_drhob, step, xc_fun_section, lsd, rho2_set, deriv_set1)
2048 : CALL update_deriv_rho(deriv_set1, [deriv_rhoa], bo, &
2049 44 : norm_drhob, gradient_cut, weight, rho1a, v_drhob%array)
2050 : CALL update_deriv_rho(deriv_set1, [deriv_rhob], bo, &
2051 44 : norm_drhob, gradient_cut, weight, rho1b, v_drhob%array)
2052 : CALL update_deriv_rho(deriv_set1, [deriv_norm_drhob], bo, &
2053 44 : norm_drhob, gradient_cut, weight, drb1drb, v_drhob%array)
2054 : CALL update_deriv_drho_ab(deriv_set1, [deriv_norm_drhoa], bo, &
2055 44 : norm_drhob, gradient_cut, weight, drb1drb, dra1dra, v_drhob%array, v_drhoa%array)
2056 : CALL update_deriv_drho_ab(deriv_set1, [deriv_norm_drho], bo, &
2057 44 : norm_drhob, gradient_cut, weight, drb1drb, dr1dr, v_drhob%array, v_drho%array)
2058 44 : IF (tau_f) THEN
2059 : CALL update_deriv_rho(deriv_set1, [deriv_tau_a], bo, &
2060 8 : norm_drhob, gradient_cut, weight, tau1a, v_drhob%array)
2061 : CALL update_deriv_rho(deriv_set1, [deriv_tau_b], bo, &
2062 8 : norm_drhob, gradient_cut, weight, tau1b, v_drhob%array)
2063 : END IF
2064 44 : IF (laplace_f) THEN
2065 : CALL update_deriv_rho(deriv_set1, [deriv_laplace_rhoa], bo, &
2066 4 : norm_drhob, gradient_cut, weight, laplace1a, v_drhob%array)
2067 : CALL update_deriv_rho(deriv_set1, [deriv_laplace_rhob], bo, &
2068 4 : norm_drhob, gradient_cut, weight, laplace1b, v_drhob%array)
2069 : END IF
2070 : END IF
2071 :
2072 44 : IF (ASSOCIATED(norm_drho)) THEN
2073 20 : CALL get_derivs_rho(norm_drho2, norm_drho, step, xc_fun_section, lsd, rho2_set, deriv_set1)
2074 : CALL update_deriv_rho(deriv_set1, [deriv_rhoa], bo, &
2075 20 : norm_drho, gradient_cut, weight, rho1a, v_drho%array)
2076 : CALL update_deriv_rho(deriv_set1, [deriv_rhob], bo, &
2077 20 : norm_drho, gradient_cut, weight, rho1b, v_drho%array)
2078 : CALL update_deriv_rho(deriv_set1, [deriv_norm_drho], bo, &
2079 20 : norm_drho, gradient_cut, weight, dr1dr, v_drho%array)
2080 : CALL update_deriv_drho_ab(deriv_set1, [deriv_norm_drhoa], bo, &
2081 20 : norm_drho, gradient_cut, weight, dr1dr, dra1dra, v_drho%array, v_drhoa%array)
2082 : CALL update_deriv_drho_ab(deriv_set1, [deriv_norm_drhob], bo, &
2083 20 : norm_drho, gradient_cut, weight, dr1dr, drb1drb, v_drho%array, v_drhob%array)
2084 20 : IF (tau_f) THEN
2085 : CALL update_deriv_rho(deriv_set1, [deriv_tau_a], bo, &
2086 8 : norm_drho, gradient_cut, weight, tau1a, v_drho%array)
2087 : CALL update_deriv_rho(deriv_set1, [deriv_tau_b], bo, &
2088 8 : norm_drho, gradient_cut, weight, tau1b, v_drho%array)
2089 : END IF
2090 20 : IF (laplace_f) THEN
2091 : CALL update_deriv_rho(deriv_set1, [deriv_laplace_rhoa], bo, &
2092 4 : norm_drho, gradient_cut, weight, laplace1a, v_drho%array)
2093 : CALL update_deriv_rho(deriv_set1, [deriv_laplace_rhob], bo, &
2094 4 : norm_drho, gradient_cut, weight, laplace1b, v_drho%array)
2095 : END IF
2096 : END IF
2097 :
2098 54 : IF (laplace_f) THEN
2099 :
2100 4 : CALL get_derivs_rho(laplace2a, laplacea, step, xc_fun_section, lsd, rho2_set, deriv_set1)
2101 :
2102 : ! Obtain the numerical 2nd derivatives w.r.t. to drho and collect the potential
2103 : CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [deriv_rhoa], bo, &
2104 4 : weight, rho1a, v_laplacea%array)
2105 : CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [deriv_rhob], bo, &
2106 4 : weight, rho1b, v_laplacea%array)
2107 4 : IF (ASSOCIATED(norm_drho)) THEN
2108 : CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [deriv_norm_drho], bo, &
2109 4 : weight, dr1dr, v_laplacea%array)
2110 : END IF
2111 4 : IF (ASSOCIATED(norm_drhoa)) THEN
2112 : CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [deriv_norm_drhoa], bo, &
2113 4 : weight, dra1dra, v_laplacea%array)
2114 : END IF
2115 4 : IF (ASSOCIATED(norm_drhob)) THEN
2116 : CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [deriv_norm_drhob], bo, &
2117 4 : weight, drb1drb, v_laplacea%array)
2118 : END IF
2119 :
2120 4 : IF (ASSOCIATED(tau1a)) THEN
2121 : CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [deriv_tau_a], bo, &
2122 4 : weight, tau1a, v_laplacea%array)
2123 : END IF
2124 4 : IF (ASSOCIATED(tau1b)) THEN
2125 : CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [deriv_tau_b], bo, &
2126 4 : weight, tau1b, v_laplacea%array)
2127 : END IF
2128 :
2129 : CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [deriv_laplace_rhoa], bo, &
2130 4 : weight, laplace1a, v_laplacea%array)
2131 :
2132 : CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [deriv_laplace_rhob], bo, &
2133 4 : weight, laplace1b, v_laplacea%array)
2134 :
2135 : ! The same for the beta spin
2136 4 : CALL get_derivs_rho(laplace2b, laplaceb, step, xc_fun_section, lsd, rho2_set, deriv_set1)
2137 :
2138 : ! Obtain the numerical 2nd derivatives w.r.t. to drho and collect the potential
2139 : CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [deriv_rhoa], bo, &
2140 4 : weight, rho1a, v_laplaceb%array)
2141 : CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [deriv_rhob], bo, &
2142 4 : weight, rho1b, v_laplaceb%array)
2143 4 : IF (ASSOCIATED(norm_drho)) THEN
2144 : CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [deriv_norm_drho], bo, &
2145 4 : weight, dr1dr, v_laplaceb%array)
2146 : END IF
2147 4 : IF (ASSOCIATED(norm_drhoa)) THEN
2148 : CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [deriv_norm_drhoa], bo, &
2149 4 : weight, dra1dra, v_laplaceb%array)
2150 : END IF
2151 4 : IF (ASSOCIATED(norm_drhob)) THEN
2152 : CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [deriv_norm_drhob], bo, &
2153 4 : weight, drb1drb, v_laplaceb%array)
2154 : END IF
2155 :
2156 4 : IF (tau_f) THEN
2157 : CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [deriv_tau_a], bo, &
2158 4 : weight, tau1a, v_laplaceb%array)
2159 : CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [deriv_tau_b], bo, &
2160 4 : weight, tau1b, v_laplaceb%array)
2161 : END IF
2162 :
2163 : CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [deriv_laplace_rhoa], bo, &
2164 4 : weight, laplace1a, v_laplaceb%array)
2165 :
2166 : CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [deriv_laplace_rhob], bo, &
2167 4 : weight, laplace1b, v_laplaceb%array)
2168 : END IF
2169 : END DO
2170 :
2171 10 : CALL virial_drho_drho(virial_pw, drhoa, v_drhoa, virial_xc)
2172 10 : CALL virial_drho_drho(virial_pw, drhob, v_drhob, virial_xc)
2173 10 : CALL virial_drho_drho(virial_pw, drho, v_drho, virial_xc)
2174 :
2175 10 : CALL deallocate_pw(v_drho, pw_pool)
2176 10 : CALL deallocate_pw(v_drhoa, pw_pool)
2177 10 : CALL deallocate_pw(v_drhob, pw_pool)
2178 :
2179 10 : IF (laplace_f) THEN
2180 15026 : virial_pw%array(:, :, :) = -rhoa(:, :, :)
2181 2 : CALL virial_laplace(virial_pw, pw_pool, virial_xc, v_laplacea%array)
2182 2 : CALL deallocate_pw(v_laplacea, pw_pool)
2183 :
2184 15026 : virial_pw%array(:, :, :) = -rhob(:, :, :)
2185 2 : CALL virial_laplace(virial_pw, pw_pool, virial_xc, v_laplaceb%array)
2186 2 : CALL deallocate_pw(v_laplaceb, pw_pool)
2187 : END IF
2188 :
2189 10 : CALL deallocate_pw(virial_pw, pw_pool)
2190 :
2191 40 : DO idir = 1, 3
2192 30 : DEALLOCATE (drho(idir)%array)
2193 40 : DEALLOCATE (drho1(idir)%array)
2194 : END DO
2195 10 : DEALLOCATE (dra1dra, drb1drb)
2196 :
2197 : ELSE
2198 26 : CALL xc_rho_set_get(rho1_set, rho=rho1, tau=tau1, laplace_rho=laplace1, can_return_null=.TRUE.)
2199 26 : CALL xc_rho_set_get(rho2_set, norm_drho=norm_drho2, laplace_rho=laplace2, can_return_null=.TRUE.)
2200 :
2201 200 : DO istep = -nsteps, nsteps
2202 174 : IF (istep == 0) CYCLE
2203 148 : weight = weights(istep, nsteps)/h
2204 148 : step = REAL(istep, dp)*h
2205 148 : CALL get_derivs_rho(norm_drho2, norm_drho, step, xc_fun_section, lsd, rho2_set, deriv_set1)
2206 :
2207 : ! Obtain the numerical 2nd derivatives w.r.t. to drho and collect the potential
2208 : CALL update_deriv_rho(deriv_set1, [deriv_rho], bo, &
2209 148 : norm_drho, gradient_cut, weight, rho1, v_drho%array)
2210 : CALL update_deriv_rho(deriv_set1, [deriv_norm_drho], bo, &
2211 148 : norm_drho, gradient_cut, weight, dr1dr, v_drho%array)
2212 :
2213 148 : IF (tau_f) THEN
2214 : CALL update_deriv_rho(deriv_set1, [deriv_tau], bo, &
2215 24 : norm_drho, gradient_cut, weight, tau1, v_drho%array)
2216 : END IF
2217 174 : IF (laplace_f) THEN
2218 : CALL update_deriv_rho(deriv_set1, [deriv_laplace_rho], bo, &
2219 12 : norm_drho, gradient_cut, weight, laplace1, v_drho%array)
2220 :
2221 12 : CALL get_derivs_rho(laplace2, laplace, step, xc_fun_section, lsd, rho2_set, deriv_set1)
2222 :
2223 : ! Obtain the numerical 2nd derivatives w.r.t. to drho and collect the potential
2224 : CALL update_deriv(deriv_set1, laplace, rho_cutoff, [deriv_rho], bo, &
2225 12 : weight, rho1, v_laplace%array)
2226 : CALL update_deriv(deriv_set1, laplace, rho_cutoff, [deriv_norm_drho], bo, &
2227 12 : weight, dr1dr, v_laplace%array)
2228 :
2229 12 : IF (tau_f) THEN
2230 : CALL update_deriv(deriv_set1, laplace, rho_cutoff, [deriv_tau], bo, &
2231 12 : weight, tau1, v_laplace%array)
2232 : END IF
2233 :
2234 : CALL update_deriv(deriv_set1, laplace, rho_cutoff, [deriv_laplace_rho], bo, &
2235 12 : weight, laplace1, v_laplace%array)
2236 : END IF
2237 : END DO
2238 :
2239 : ! Calculate the virial contribution from the potential
2240 26 : CALL virial_drho_drho(virial_pw, drho, v_drho, virial_xc)
2241 :
2242 26 : CALL deallocate_pw(v_drho, pw_pool)
2243 :
2244 26 : IF (laplace_f) THEN
2245 28862 : virial_pw%array(:, :, :) = -rho(:, :, :)
2246 2 : CALL virial_laplace(virial_pw, pw_pool, virial_xc, v_laplace%array)
2247 2 : CALL deallocate_pw(v_laplace, pw_pool)
2248 : END IF
2249 :
2250 26 : CALL deallocate_pw(virial_pw, pw_pool)
2251 : END IF
2252 :
2253 : END IF
2254 :
2255 36 : CALL xc_dset_release(deriv_set1)
2256 :
2257 36 : DEALLOCATE (dr1dr)
2258 :
2259 36 : CALL xc_rho_set_release(rho1_set)
2260 36 : CALL xc_rho_set_release(rho2_set)
2261 : END IF
2262 :
2263 2012 : DO ispin = 1, SIZE(rho_r)
2264 2012 : CALL pw_pool%give_back_pw(rho_r(ispin))
2265 : END DO
2266 832 : DEALLOCATE (rho_r)
2267 :
2268 832 : IF (ASSOCIATED(tau_r)) THEN
2269 254 : DO ispin = 1, SIZE(tau_r)
2270 254 : CALL pw_pool%give_back_pw(tau_r(ispin))
2271 : END DO
2272 100 : DEALLOCATE (tau_r)
2273 : END IF
2274 :
2275 832 : CALL timestop(handle)
2276 :
2277 37440 : END SUBROUTINE xc_calc_2nd_deriv_numerical
2278 :
2279 : ! **************************************************************************************************
2280 : !> \brief ...
2281 : !> \param rho_r ...
2282 : !> \param rho_g ...
2283 : !> \param rho1_r ...
2284 : !> \param rhoa ...
2285 : !> \param rhob ...
2286 : !> \param vxc_rho ...
2287 : !> \param tau_r ...
2288 : !> \param tau1_r ...
2289 : !> \param tau_a ...
2290 : !> \param tau_b ...
2291 : !> \param vxc_tau ...
2292 : !> \param xc_section ...
2293 : !> \param pw_pool ...
2294 : !> \param step ...
2295 : ! **************************************************************************************************
2296 1848 : SUBROUTINE calc_resp_potential_numer_ab(rho_r, rho_g, rho1_r, rhoa, rhob, vxc_rho, &
2297 : tau_r, tau1_r, tau_a, tau_b, vxc_tau, &
2298 : xc_section, pw_pool, step)
2299 :
2300 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER, INTENT(IN) :: vxc_rho, vxc_tau
2301 : TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(IN) :: rho1_r
2302 : TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(IN), POINTER :: tau1_r
2303 : TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool
2304 : TYPE(section_vals_type), INTENT(IN), POINTER :: xc_section
2305 : REAL(KIND=dp), INTENT(IN) :: step
2306 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER, INTENT(IN) :: rhoa, rhob, tau_a, tau_b
2307 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER, INTENT(IN) :: rho_r
2308 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
2309 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: tau_r
2310 :
2311 : CHARACTER(len=*), PARAMETER :: routineN = 'calc_resp_potential_numer_ab'
2312 :
2313 : INTEGER :: handle
2314 : REAL(KIND=dp) :: exc
2315 : REAL(KIND=dp), DIMENSION(3, 3) :: virial_dummy
2316 :
2317 1848 : CALL timeset(routineN, handle)
2318 :
2319 1848 : !$OMP PARALLEL DEFAULT(NONE) SHARED(rho_r,rhoa,rhob,step,rho1_r,tau_r,tau_a,tau_b,tau1_r)
2320 : !$OMP WORKSHARE
2321 : rho_r(1)%array(:, :, :) = rhoa(:, :, :) + step*rho1_r(1)%array(:, :, :)
2322 : !$OMP END WORKSHARE NOWAIT
2323 : !$OMP WORKSHARE
2324 : rho_r(2)%array(:, :, :) = rhob(:, :, :) + step*rho1_r(2)%array(:, :, :)
2325 : !$OMP END WORKSHARE NOWAIT
2326 : IF (ASSOCIATED(tau1_r) .AND. ASSOCIATED(tau_r) .AND. ASSOCIATED(tau_a) .AND. ASSOCIATED(tau_b)) THEN
2327 : !$OMP WORKSHARE
2328 : tau_r(1)%array(:, :, :) = tau_a(:, :, :) + step*tau1_r(1)%array(:, :, :)
2329 : !$OMP END WORKSHARE NOWAIT
2330 : !$OMP WORKSHARE
2331 : tau_r(2)%array(:, :, :) = tau_b(:, :, :) + step*tau1_r(2)%array(:, :, :)
2332 : !$OMP END WORKSHARE NOWAIT
2333 : END IF
2334 : !$OMP END PARALLEL
2335 : CALL xc_vxc_pw_create(vxc_rho, vxc_tau, exc, rho_r, rho_g, tau_r, xc_section, &
2336 1848 : pw_pool, .FALSE., virial_dummy)
2337 :
2338 1848 : CALL timestop(handle)
2339 :
2340 1848 : END SUBROUTINE calc_resp_potential_numer_ab
2341 :
2342 : ! **************************************************************************************************
2343 : !> \brief calculates stress tensor and potential contributions from the first derivative
2344 : !> \param deriv_set ...
2345 : !> \param description ...
2346 : !> \param virial_pw ...
2347 : !> \param drho ...
2348 : !> \param drho1 ...
2349 : !> \param virial_xc ...
2350 : !> \param norm_drho ...
2351 : !> \param gradient_cut ...
2352 : !> \param dr1dr ...
2353 : !> \param v_drho ...
2354 : ! **************************************************************************************************
2355 52 : SUBROUTINE apply_drho(deriv_set, description, virial_pw, drho, drho1, virial_xc, norm_drho, gradient_cut, dr1dr, v_drho)
2356 :
2357 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
2358 : INTEGER, DIMENSION(:), INTENT(in) :: description
2359 : TYPE(pw_r3d_rs_type), INTENT(IN) :: virial_pw
2360 : TYPE(cp_3d_r_cp_type), DIMENSION(3), INTENT(IN) :: drho, drho1
2361 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT) :: virial_xc
2362 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: norm_drho
2363 : REAL(KIND=dp), INTENT(IN) :: gradient_cut
2364 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: dr1dr
2365 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: v_drho
2366 :
2367 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_drho'
2368 :
2369 : INTEGER :: handle
2370 52 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: deriv_data
2371 : TYPE(xc_derivative_type), POINTER :: deriv_att
2372 :
2373 52 : CALL timeset(routineN, handle)
2374 :
2375 52 : deriv_att => xc_dset_get_derivative(deriv_set, description)
2376 52 : IF (ASSOCIATED(deriv_att)) THEN
2377 52 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
2378 52 : CALL virial_drho_drho1(virial_pw, drho, drho1, deriv_data, virial_xc)
2379 :
2380 52 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(dr1dr,gradient_cut,norm_drho,v_drho,deriv_data)
2381 : v_drho(:, :, :) = v_drho(:, :, :) + &
2382 : deriv_data(:, :, :)*dr1dr(:, :, :)/MAX(gradient_cut, norm_drho(:, :, :))**2
2383 : !$OMP END PARALLEL WORKSHARE
2384 : END IF
2385 :
2386 52 : CALL timestop(handle)
2387 :
2388 52 : END SUBROUTINE apply_drho
2389 :
2390 : ! **************************************************************************************************
2391 : !> \brief adds potential contributions from derivatives of rho or diagonal terms of norm_drho
2392 : !> \param deriv_set1 ...
2393 : !> \param description ...
2394 : !> \param bo ...
2395 : !> \param norm_drho norm_drho of which derivative is calculated
2396 : !> \param gradient_cut ...
2397 : !> \param h ...
2398 : !> \param rho1 function to contract the derivative with (rho1 for rho, dr1dr for norm_drho)
2399 : !> \param v_drho ...
2400 : ! **************************************************************************************************
2401 728 : SUBROUTINE update_deriv_rho(deriv_set1, description, bo, norm_drho, gradient_cut, weight, rho1, v_drho)
2402 :
2403 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set1
2404 : INTEGER, DIMENSION(:), INTENT(in) :: description
2405 : INTEGER, DIMENSION(2, 3), INTENT(IN) :: bo
2406 : REAL(KIND=dp), DIMENSION(bo(1, 1):bo(2, 1), bo(1, &
2407 : 2):bo(2, 2), bo(1, 3):bo(2, 3)), INTENT(IN) :: norm_drho
2408 : REAL(KIND=dp), INTENT(IN) :: gradient_cut, weight
2409 : REAL(KIND=dp), DIMENSION(bo(1, 1):bo(2, 1), bo(1, &
2410 : 2):bo(2, 2), bo(1, 3):bo(2, 3)), INTENT(IN) :: rho1
2411 : REAL(KIND=dp), DIMENSION(bo(1, 1):bo(2, 1), bo(1, &
2412 : 2):bo(2, 2), bo(1, 3):bo(2, 3)), INTENT(INOUT) :: v_drho
2413 :
2414 : CHARACTER(len=*), PARAMETER :: routineN = 'update_deriv_rho'
2415 :
2416 : INTEGER :: handle, i, j, k
2417 : REAL(KIND=dp) :: de
2418 728 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: deriv_data1
2419 : TYPE(xc_derivative_type), POINTER :: deriv_att1
2420 :
2421 728 : CALL timeset(routineN, handle)
2422 :
2423 : ! Obtain the numerical 2nd derivatives w.r.t. to drho and collect the potential
2424 728 : deriv_att1 => xc_dset_get_derivative(deriv_set1, description)
2425 728 : IF (ASSOCIATED(deriv_att1)) THEN
2426 728 : CALL xc_derivative_get(deriv_att1, deriv_data=deriv_data1)
2427 : !$OMP PARALLEL DO DEFAULT(NONE) &
2428 : !$OMP SHARED(bo,deriv_data1,weight,norm_drho,v_drho,rho1,gradient_cut) &
2429 : !$OMP PRIVATE(i,j,k,de) &
2430 728 : !$OMP COLLAPSE(3)
2431 : DO k = bo(1, 3), bo(2, 3)
2432 : DO j = bo(1, 2), bo(2, 2)
2433 : DO i = bo(1, 1), bo(2, 1)
2434 : de = weight*deriv_data1(i, j, k)/MAX(gradient_cut, norm_drho(i, j, k))**2
2435 : v_drho(i, j, k) = v_drho(i, j, k) - de*rho1(i, j, k)
2436 : END DO
2437 : END DO
2438 : END DO
2439 : !$OMP END PARALLEL DO
2440 : END IF
2441 :
2442 728 : CALL timestop(handle)
2443 :
2444 728 : END SUBROUTINE update_deriv_rho
2445 :
2446 : ! **************************************************************************************************
2447 : !> \brief adds potential contributions from derivatives of a component with positive and negative values
2448 : !> \param deriv_set1 ...
2449 : !> \param description ...
2450 : !> \param bo ...
2451 : !> \param h ...
2452 : !> \param rho1 function to contract the derivative with (rho1 for rho, dr1dr for norm_drho)
2453 : !> \param v ...
2454 : ! **************************************************************************************************
2455 120 : SUBROUTINE update_deriv(deriv_set1, rho, rho_cutoff, description, bo, weight, rho1, v)
2456 :
2457 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set1
2458 : INTEGER, DIMENSION(:), INTENT(in) :: description
2459 : INTEGER, DIMENSION(2, 3), INTENT(IN) :: bo
2460 : REAL(KIND=dp), INTENT(IN) :: weight, rho_cutoff
2461 : REAL(KIND=dp), DIMENSION(bo(1, 1):bo(2, 1), bo(1, &
2462 : 2):bo(2, 2), bo(1, 3):bo(2, 3)), INTENT(IN) :: rho, rho1
2463 : REAL(KIND=dp), DIMENSION(bo(1, 1):bo(2, 1), bo(1, &
2464 : 2):bo(2, 2), bo(1, 3):bo(2, 3)), INTENT(INOUT) :: v
2465 :
2466 : CHARACTER(len=*), PARAMETER :: routineN = 'update_deriv'
2467 :
2468 : INTEGER :: handle, i, j, k
2469 : REAL(KIND=dp) :: de
2470 120 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: deriv_data1
2471 : TYPE(xc_derivative_type), POINTER :: deriv_att1
2472 :
2473 120 : CALL timeset(routineN, handle)
2474 :
2475 : ! Obtain the numerical 2nd derivatives w.r.t. to drho and collect the potential
2476 120 : deriv_att1 => xc_dset_get_derivative(deriv_set1, description)
2477 120 : IF (ASSOCIATED(deriv_att1)) THEN
2478 120 : CALL xc_derivative_get(deriv_att1, deriv_data=deriv_data1)
2479 : !$OMP PARALLEL DO DEFAULT(NONE) &
2480 : !$OMP SHARED(bo,deriv_data1,weight,v,rho1,rho, rho_cutoff) &
2481 : !$OMP PRIVATE(i,j,k,de) &
2482 120 : !$OMP COLLAPSE(3)
2483 : DO k = bo(1, 3), bo(2, 3)
2484 : DO j = bo(1, 2), bo(2, 2)
2485 : DO i = bo(1, 1), bo(2, 1)
2486 : ! We have to consider that the given density (mostly the Laplacian) may have positive and negative values
2487 : de = weight*deriv_data1(i, j, k)/SIGN(MAX(ABS(rho(i, j, k)), rho_cutoff), rho(i, j, k))
2488 : v(i, j, k) = v(i, j, k) + de*rho1(i, j, k)
2489 : END DO
2490 : END DO
2491 : END DO
2492 : !$OMP END PARALLEL DO
2493 : END IF
2494 :
2495 120 : CALL timestop(handle)
2496 :
2497 120 : END SUBROUTINE update_deriv
2498 :
2499 : ! **************************************************************************************************
2500 : !> \brief adds mixed derivatives of norm_drho
2501 : !> \param deriv_set1 ...
2502 : !> \param description ...
2503 : !> \param bo ...
2504 : !> \param norm_drhoa norm_drho of which derivatives is calculated
2505 : !> \param gradient_cut ...
2506 : !> \param h ...
2507 : !> \param dra1dra dr1dr corresponding to norm_drho
2508 : !> \param drb1drb ...
2509 : !> \param v_drhoa potential corresponding to norm_drho
2510 : !> \param v_drhob ...
2511 : ! **************************************************************************************************
2512 216 : SUBROUTINE update_deriv_drho_ab(deriv_set1, description, bo, &
2513 216 : norm_drhoa, gradient_cut, weight, dra1dra, drb1drb, v_drhoa, v_drhob)
2514 :
2515 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set1
2516 : INTEGER, DIMENSION(:), INTENT(in) :: description
2517 : INTEGER, DIMENSION(2, 3), INTENT(IN) :: bo
2518 : REAL(KIND=dp), DIMENSION(bo(1, 1):bo(2, 1), bo(1, &
2519 : 2):bo(2, 2), bo(1, 3):bo(2, 3)), INTENT(IN) :: norm_drhoa
2520 : REAL(KIND=dp), INTENT(IN) :: gradient_cut, weight
2521 : REAL(KIND=dp), DIMENSION(bo(1, 1):bo(2, 1), bo(1, &
2522 : 2):bo(2, 2), bo(1, 3):bo(2, 3)), INTENT(IN) :: dra1dra, drb1drb
2523 : REAL(KIND=dp), DIMENSION(bo(1, 1):bo(2, 1), bo(1, &
2524 : 2):bo(2, 2), bo(1, 3):bo(2, 3)), INTENT(INOUT) :: v_drhoa, v_drhob
2525 :
2526 : CHARACTER(len=*), PARAMETER :: routineN = 'update_deriv_drho_ab'
2527 :
2528 : INTEGER :: handle, i, j, k
2529 : REAL(KIND=dp) :: de
2530 216 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: deriv_data1
2531 : TYPE(xc_derivative_type), POINTER :: deriv_att1
2532 :
2533 216 : CALL timeset(routineN, handle)
2534 :
2535 216 : deriv_att1 => xc_dset_get_derivative(deriv_set1, description)
2536 216 : IF (ASSOCIATED(deriv_att1)) THEN
2537 168 : CALL xc_derivative_get(deriv_att1, deriv_data=deriv_data1)
2538 : !$OMP PARALLEL DO DEFAULT(NONE) &
2539 : !$OMP PRIVATE(k,j,i,de) &
2540 : !$OMP SHARED(bo,drb1drb,dra1dra,deriv_data1,weight,gradient_cut,norm_drhoa,v_drhoa,v_drhob) &
2541 168 : !$OMP COLLAPSE(3)
2542 : DO k = bo(1, 3), bo(2, 3)
2543 : DO j = bo(1, 2), bo(2, 2)
2544 : DO i = bo(1, 1), bo(2, 1)
2545 : ! We introduce a factor of two because we will average between both numerical derivatives
2546 : de = 0.5_dp*weight*deriv_data1(i, j, k)/MAX(gradient_cut, norm_drhoa(i, j, k))**2
2547 : v_drhoa(i, j, k) = v_drhoa(i, j, k) - de*drb1drb(i, j, k)
2548 : v_drhob(i, j, k) = v_drhob(i, j, k) - de*dra1dra(i, j, k)
2549 : END DO
2550 : END DO
2551 : END DO
2552 : !$OMP END PARALLEL DO
2553 : END IF
2554 :
2555 216 : CALL timestop(handle)
2556 :
2557 216 : END SUBROUTINE update_deriv_drho_ab
2558 :
2559 : ! **************************************************************************************************
2560 : !> \brief calculate derivative sets for helper points
2561 : !> \param norm_drho2 norm_drho of new points
2562 : !> \param norm_drho norm_drho of KS density
2563 : !> \param h ...
2564 : !> \param xc_fun_section ...
2565 : !> \param lsd ...
2566 : !> \param rho2_set rho_set for new points
2567 : !> \param deriv_set1 will contain derivatives of the perturbed density
2568 : ! **************************************************************************************************
2569 276 : SUBROUTINE get_derivs_rho(norm_drho2, norm_drho, step, xc_fun_section, lsd, rho2_set, deriv_set1)
2570 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: norm_drho2
2571 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: norm_drho
2572 : REAL(KIND=dp), INTENT(IN) :: step
2573 : TYPE(section_vals_type), INTENT(IN), POINTER :: xc_fun_section
2574 : LOGICAL, INTENT(IN) :: lsd
2575 : TYPE(xc_rho_set_type), INTENT(INOUT) :: rho2_set
2576 : TYPE(xc_derivative_set_type) :: deriv_set1
2577 :
2578 : CHARACTER(len=*), PARAMETER :: routineN = 'get_derivs_rho'
2579 :
2580 : INTEGER :: handle
2581 :
2582 276 : CALL timeset(routineN, handle)
2583 :
2584 : ! Copy the densities, do one step into the direction of drho
2585 276 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(norm_drho,norm_drho2,step)
2586 : norm_drho2 = norm_drho*(1.0_dp + step)
2587 : !$OMP END PARALLEL WORKSHARE
2588 :
2589 276 : CALL xc_dset_zero_all(deriv_set1)
2590 :
2591 : ! Calculate the derivatives of the functional
2592 : CALL xc_functionals_eval(xc_fun_section, &
2593 : lsd=lsd, &
2594 : rho_set=rho2_set, &
2595 : deriv_set=deriv_set1, &
2596 276 : deriv_order=1)
2597 :
2598 : ! Return to the original values
2599 276 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(norm_drho,norm_drho2)
2600 : norm_drho2 = norm_drho
2601 : !$OMP END PARALLEL WORKSHARE
2602 :
2603 276 : CALL divide_by_norm_drho(deriv_set1, rho2_set, lsd)
2604 :
2605 276 : CALL timestop(handle)
2606 :
2607 276 : END SUBROUTINE get_derivs_rho
2608 :
2609 : ! **************************************************************************************************
2610 : !> \brief Calculates the second derivative of E_xc at rho in the direction
2611 : !> rho1 (if you see the second derivative as bilinear form)
2612 : !> partial_rho|_(rho=rho) partial_rho|_(rho=rho) E_xc drho(rho1)drho
2613 : !> The other direction is still undetermined, thus it returns
2614 : !> a potential (partial integration is performed to reduce it to
2615 : !> function of rho, removing the dependence from its partial derivs)
2616 : !> Has to be called after the setup by xc_prep_2nd_deriv.
2617 : !> \param v_xc exchange-correlation potential
2618 : !> \param v_xc_tau ...
2619 : !> \param deriv_set derivatives of the exchange-correlation potential
2620 : !> \param rho_set object containing the density at which the derivatives were calculated
2621 : !> \param rho1_set object containing the density with which to fold
2622 : !> \param pw_pool the pool for the grids
2623 : !> \param xc_section XC parameters
2624 : !> \param gapw Gaussian and augmented plane waves calculation
2625 : !> \param vxg ...
2626 : !> \param tddfpt_fac factor that multiplies the crossterms (tddfpt triplets
2627 : !> on a closed shell system it should be -1, defaults to 1)
2628 : !> \param compute_virial ...
2629 : !> \param virial_xc ...
2630 : !> \note
2631 : !> The old version of this routine was smarter: it handled split_desc(1)
2632 : !> and split_desc(2) separately, thus the code automatically handled all
2633 : !> possible cross terms (you only had to check if it was diagonal to avoid
2634 : !> double counting). I think that is the way to go if you want to add more
2635 : !> terms (tau,rho in LSD,...). The problem with the old code was that it
2636 : !> because of the old functional structure it sometime guessed wrongly
2637 : !> which derivative was where. There were probably still bugs with gradient
2638 : !> corrected functionals (never tested), and it didn't contain first
2639 : !> derivatives with respect to drho (that contribute also to the second
2640 : !> derivative wrt. rho).
2641 : !> The code was a little complex because it really tried to handle any
2642 : !> functional derivative in the most efficient way with the given contents of
2643 : !> rho_set.
2644 : !> Anyway I strongly encourage whoever wants to modify this code to give a
2645 : !> look to the old version. [fawzi]
2646 : ! **************************************************************************************************
2647 31770 : SUBROUTINE xc_calc_2nd_deriv_analytical(v_xc, v_xc_tau, deriv_set, rho_set, rho1_set, &
2648 : pw_pool, xc_section, gapw, vxg, tddfpt_fac, &
2649 : compute_virial, virial_xc)
2650 :
2651 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_xc, v_xc_tau
2652 : TYPE(xc_derivative_set_type) :: deriv_set
2653 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set, rho1_set
2654 : TYPE(pw_pool_type), POINTER :: pw_pool
2655 : TYPE(section_vals_type), POINTER :: xc_section
2656 : LOGICAL, INTENT(IN), OPTIONAL :: gapw
2657 : REAL(kind=dp), DIMENSION(:, :, :, :), OPTIONAL, &
2658 : POINTER :: vxg
2659 : REAL(kind=dp), INTENT(in), OPTIONAL :: tddfpt_fac
2660 : LOGICAL, INTENT(IN), OPTIONAL :: compute_virial
2661 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT), &
2662 : OPTIONAL :: virial_xc
2663 :
2664 : CHARACTER(len=*), PARAMETER :: routineN = 'xc_calc_2nd_deriv_analytical'
2665 :
2666 : INTEGER :: handle, i, ia, idir, ir, ispin, j, jdir, &
2667 : k, nspins, xc_deriv_method_id
2668 : INTEGER, DIMENSION(2, 3) :: bo
2669 : LOGICAL :: gradient_f, lsd, my_compute_virial, &
2670 : my_gapw, tau_f, laplace_f, rho_f
2671 : REAL(KIND=dp) :: fac, gradient_cut, tmp, factor2
2672 31770 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: dr1dr, dra1dra, drb1drb
2673 63540 : REAL(kind=dp), DIMENSION(:, :, :), POINTER :: deriv_data, e_drhoa, e_drhob, &
2674 31770 : e_drho, norm_drho, norm_drhoa, &
2675 31770 : norm_drhob, rho1, rho1a, rho1b, &
2676 31770 : tau1, tau1a, tau1b, laplace1, laplace1a, laplace1b, &
2677 31770 : rho, rhoa, rhob
2678 603630 : TYPE(cp_3d_r_cp_type), DIMENSION(3) :: drho, drho1, drho1a, drho1b, drhoa, drhob
2679 31770 : TYPE(pw_r3d_rs_type), DIMENSION(:), ALLOCATABLE :: v_drhoa, v_drhob, v_drho, v_laplace
2680 31770 : TYPE(pw_r3d_rs_type), DIMENSION(:, :), ALLOCATABLE :: v_drho_r
2681 : TYPE(pw_r3d_rs_type) :: virial_pw
2682 : TYPE(pw_c1d_gs_type) :: tmp_g, vxc_g
2683 : TYPE(xc_derivative_type), POINTER :: deriv_att
2684 :
2685 31770 : CALL timeset(routineN, handle)
2686 :
2687 31770 : NULLIFY (e_drhoa, e_drhob, e_drho)
2688 :
2689 31770 : my_gapw = .FALSE.
2690 31770 : IF (PRESENT(gapw)) my_gapw = gapw
2691 :
2692 31770 : my_compute_virial = .FALSE.
2693 31770 : IF (PRESENT(compute_virial)) my_compute_virial = compute_virial
2694 :
2695 31770 : CPASSERT(ASSOCIATED(v_xc))
2696 31770 : CPASSERT(ASSOCIATED(xc_section))
2697 31770 : IF (my_gapw) THEN
2698 10060 : CPASSERT(PRESENT(vxg))
2699 : END IF
2700 31770 : IF (my_compute_virial) THEN
2701 358 : CPASSERT(PRESENT(virial_xc))
2702 : END IF
2703 :
2704 : CALL section_vals_val_get(xc_section, "XC_GRID%XC_DERIV", &
2705 31770 : i_val=xc_deriv_method_id)
2706 31770 : CALL xc_rho_set_get(rho_set, drho_cutoff=gradient_cut)
2707 31770 : nspins = SIZE(v_xc)
2708 31770 : lsd = ASSOCIATED(rho_set%rhoa)
2709 31770 : fac = 0.0_dp
2710 31770 : factor2 = 1.0_dp
2711 31770 : IF (PRESENT(tddfpt_fac)) fac = tddfpt_fac
2712 31770 : IF (PRESENT(tddfpt_fac)) factor2 = tddfpt_fac
2713 :
2714 317700 : bo = rho_set%local_bounds
2715 :
2716 31770 : CALL check_for_derivatives(deriv_set, lsd, rho_f, gradient_f, tau_f, laplace_f)
2717 :
2718 31770 : IF (tau_f) THEN
2719 636 : CPASSERT(ASSOCIATED(v_xc_tau))
2720 : END IF
2721 :
2722 31770 : IF (gradient_f) THEN
2723 204130 : ALLOCATE (v_drho_r(3, nspins), v_drho(nspins))
2724 40826 : DO ispin = 1, nspins
2725 84288 : DO idir = 1, 3
2726 84288 : CALL allocate_pw(v_drho_r(idir, ispin), pw_pool, bo)
2727 : END DO
2728 40826 : CALL allocate_pw(v_drho(ispin), pw_pool, bo)
2729 : END DO
2730 :
2731 19754 : IF (xc_requires_tmp_g(xc_deriv_method_id) .AND. .NOT. my_gapw) THEN
2732 12060 : IF (ASSOCIATED(pw_pool)) THEN
2733 12060 : CALL pw_pool%create_pw(tmp_g)
2734 12060 : CALL pw_pool%create_pw(vxc_g)
2735 : ELSE
2736 : ! remember to refix for gapw
2737 0 : CPABORT("XC_DERIV method is not implemented in GAPW")
2738 : END IF
2739 : END IF
2740 : END IF
2741 :
2742 66840 : DO ispin = 1, nspins
2743 993206681 : v_xc(ispin)%array = 0.0_dp
2744 : END DO
2745 :
2746 31770 : IF (tau_f) THEN
2747 1346 : DO ispin = 1, nspins
2748 14404518 : v_xc_tau(ispin)%array = 0.0_dp
2749 : END DO
2750 : END IF
2751 :
2752 31770 : IF (laplace_f .AND. my_gapw) &
2753 0 : CPABORT("Laplace-dependent functional not implemented with GAPW!")
2754 :
2755 31770 : IF (my_compute_virial .AND. (gradient_f .OR. laplace_f)) CALL allocate_pw(virial_pw, pw_pool, bo)
2756 :
2757 31770 : IF (lsd) THEN
2758 :
2759 : !-------------------!
2760 : ! UNrestricted case !
2761 : !-------------------!
2762 :
2763 4594 : CALL xc_rho_set_get(rho1_set, rhoa=rho1a, rhob=rho1b)
2764 :
2765 4594 : IF (gradient_f) THEN
2766 : CALL xc_rho_set_get(rho_set, drhoa=drhoa, drhob=drhob, &
2767 2188 : norm_drho=norm_drho, norm_drhoa=norm_drhoa, norm_drhob=norm_drhob)
2768 2188 : CALL xc_rho_set_get(rho1_set, drhoa=drho1a, drhob=drho1b)
2769 :
2770 2188 : CALL calc_drho_from_ab(drho, drhoa, drhob)
2771 2188 : CALL calc_drho_from_ab(drho1, drho1a, drho1b)
2772 :
2773 2188 : CALL prepare_dr1dr(dra1dra, drhoa, drho1a)
2774 2188 : IF (nspins /= 1) THEN
2775 1318 : CALL prepare_dr1dr(drb1drb, drhob, drho1b)
2776 1318 : CALL prepare_dr1dr(dr1dr, drho, drho1)
2777 : ELSE
2778 870 : CALL prepare_dr1dr(drb1drb, drhob, drho1b)
2779 870 : CALL prepare_dr1dr_ab(dr1dr, drhoa, drhob, drho1a, drho1b, fac)
2780 : END IF
2781 :
2782 15764 : ALLOCATE (v_drhoa(nspins), v_drhob(nspins))
2783 5694 : DO ispin = 1, nspins
2784 3506 : CALL allocate_pw(v_drhoa(ispin), pw_pool, bo)
2785 5694 : CALL allocate_pw(v_drhob(ispin), pw_pool, bo)
2786 : END DO
2787 :
2788 : END IF
2789 :
2790 4594 : IF (laplace_f) THEN
2791 38 : CALL xc_rho_set_get(rho1_set, laplace_rhoa=laplace1a, laplace_rhob=laplace1b)
2792 :
2793 190 : ALLOCATE (v_laplace(nspins))
2794 114 : DO ispin = 1, nspins
2795 114 : CALL allocate_pw(v_laplace(ispin), pw_pool, bo)
2796 : END DO
2797 :
2798 38 : IF (my_compute_virial) CALL xc_rho_set_get(rho_set, rhoa=rhoa, rhob=rhob)
2799 : END IF
2800 :
2801 4594 : IF (tau_f) THEN
2802 74 : CALL xc_rho_set_get(rho1_set, tau_a=tau1a, tau_b=tau1b)
2803 : END IF
2804 :
2805 4594 : IF (nspins /= 1) THEN
2806 :
2807 33348 : $:add_2nd_derivative_terms(arguments_openshell)
2808 :
2809 : ELSE
2810 :
2811 1294 : $:add_2nd_derivative_terms(arguments_triplet_outer, arguments_triplet_inner)
2812 :
2813 : END IF
2814 :
2815 4594 : IF (gradient_f) THEN
2816 :
2817 2188 : IF (my_compute_virial) THEN
2818 10 : CALL virial_drho_drho(virial_pw, drhoa, v_drhoa(1), virial_xc)
2819 10 : CALL virial_drho_drho(virial_pw, drhob, v_drhob(2), virial_xc)
2820 40 : DO idir = 1, 3
2821 30 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(drho,idir,v_drho,virial_pw)
2822 : virial_pw%array(:, :, :) = drho(idir)%array(:, :, :)*(v_drho(1)%array(:, :, :) + v_drho(2)%array(:, :, :))
2823 : !$OMP END PARALLEL WORKSHARE
2824 100 : DO jdir = 1, idir
2825 : tmp = -0.5_dp*virial_pw%pw_grid%dvol*accurate_dot_product(virial_pw%array(:, :, :), &
2826 60 : drho(jdir)%array(:, :, :))
2827 60 : virial_xc(jdir, idir) = virial_xc(jdir, idir) + tmp
2828 90 : virial_xc(idir, jdir) = virial_xc(jdir, idir)
2829 : END DO
2830 : END DO
2831 : END IF ! my_compute_virial
2832 :
2833 2188 : IF (my_gapw) THEN
2834 : !$OMP PARALLEL DO DEFAULT(NONE) &
2835 : !$OMP PRIVATE(ia,idir,ispin,ir) &
2836 : !$OMP SHARED(bo,nspins,vxg,drhoa,drhob,v_drhoa,v_drhob,v_drho, &
2837 : !$OMP e_drhoa,e_drhob,e_drho,drho1a,drho1b,fac,drho,drho1) &
2838 602 : !$OMP COLLAPSE(3)
2839 : DO ir = bo(1, 2), bo(2, 2)
2840 : DO ia = bo(1, 1), bo(2, 1)
2841 : DO idir = 1, 3
2842 : DO ispin = 1, nspins
2843 : vxg(idir, ia, ir, ispin) = &
2844 : -(v_drhoa(ispin)%array(ia, ir, 1)*drhoa(idir)%array(ia, ir, 1) + &
2845 : v_drhob(ispin)%array(ia, ir, 1)*drhob(idir)%array(ia, ir, 1) + &
2846 : v_drho(ispin)%array(ia, ir, 1)*drho(idir)%array(ia, ir, 1))
2847 : END DO
2848 : IF (ASSOCIATED(e_drhoa)) THEN
2849 : vxg(idir, ia, ir, 1) = vxg(idir, ia, ir, 1) + &
2850 : e_drhoa(ia, ir, 1)*drho1a(idir)%array(ia, ir, 1)
2851 : END IF
2852 : IF (nspins /= 1 .AND. ASSOCIATED(e_drhob)) THEN
2853 : vxg(idir, ia, ir, 2) = vxg(idir, ia, ir, 2) + &
2854 : e_drhob(ia, ir, 1)*drho1b(idir)%array(ia, ir, 1)
2855 : END IF
2856 : IF (ASSOCIATED(e_drho)) THEN
2857 : IF (nspins /= 1) THEN
2858 : vxg(idir, ia, ir, 1) = vxg(idir, ia, ir, 1) + &
2859 : e_drho(ia, ir, 1)*drho1(idir)%array(ia, ir, 1)
2860 : vxg(idir, ia, ir, 2) = vxg(idir, ia, ir, 2) + &
2861 : e_drho(ia, ir, 1)*drho1(idir)%array(ia, ir, 1)
2862 : ELSE
2863 : vxg(idir, ia, ir, 1) = vxg(idir, ia, ir, 1) + &
2864 : e_drho(ia, ir, 1)*(drho1a(idir)%array(ia, ir, 1) + &
2865 : fac*drho1b(idir)%array(ia, ir, 1))
2866 : END IF
2867 : END IF
2868 : END DO
2869 : END DO
2870 : END DO
2871 : !$OMP END PARALLEL DO
2872 : ELSE
2873 :
2874 : ! partial integration
2875 6344 : DO idir = 1, 3
2876 :
2877 13086 : DO ispin = 1, nspins
2878 13086 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(v_drho_r,v_drhoa,v_drhob,v_drho,drhoa,drhob,drho,ispin,idir)
2879 : v_drho_r(idir, ispin)%array(:, :, :) = &
2880 : v_drhoa(ispin)%array(:, :, :)*drhoa(idir)%array(:, :, :) + &
2881 : v_drhob(ispin)%array(:, :, :)*drhob(idir)%array(:, :, :) + &
2882 : v_drho(ispin)%array(:, :, :)*drho(idir)%array(:, :, :)
2883 : !$OMP END PARALLEL WORKSHARE
2884 : END DO
2885 4758 : IF (ASSOCIATED(e_drhoa)) THEN
2886 4758 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(v_drho_r,e_drhoa,drho1a,idir)
2887 : v_drho_r(idir, 1)%array(:, :, :) = v_drho_r(idir, 1)%array(:, :, :) - &
2888 : e_drhoa(:, :, :)*drho1a(idir)%array(:, :, :)
2889 : !$OMP END PARALLEL WORKSHARE
2890 : END IF
2891 4758 : IF (nspins /= 1 .AND. ASSOCIATED(e_drhob)) THEN
2892 3570 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(v_drho_r,e_drhob,drho1b,idir)
2893 : v_drho_r(idir, 2)%array(:, :, :) = v_drho_r(idir, 2)%array(:, :, :) - &
2894 : e_drhob(:, :, :)*drho1b(idir)%array(:, :, :)
2895 : !$OMP END PARALLEL WORKSHARE
2896 : END IF
2897 6344 : IF (ASSOCIATED(e_drho)) THEN
2898 : !$OMP PARALLEL DO DEFAULT(NONE) &
2899 : !$OMP PRIVATE(k,j,i) &
2900 : !$OMP SHARED(bo,v_drho_r,e_drho,drho1a,drho1b,drho1,fac,idir,nspins) &
2901 4758 : !$OMP COLLAPSE(3)
2902 : DO k = bo(1, 3), bo(2, 3)
2903 : DO j = bo(1, 2), bo(2, 2)
2904 : DO i = bo(1, 1), bo(2, 1)
2905 : IF (nspins /= 1) THEN
2906 : v_drho_r(idir, 1)%array(i, j, k) = v_drho_r(idir, 1)%array(i, j, k) - &
2907 : e_drho(i, j, k)*drho1(idir)%array(i, j, k)
2908 : v_drho_r(idir, 2)%array(i, j, k) = v_drho_r(idir, 2)%array(i, j, k) - &
2909 : e_drho(i, j, k)*drho1(idir)%array(i, j, k)
2910 : ELSE
2911 : v_drho_r(idir, 1)%array(i, j, k) = v_drho_r(idir, 1)%array(i, j, k) - &
2912 : e_drho(i, j, k)*(drho1a(idir)%array(i, j, k) + &
2913 : fac*drho1b(idir)%array(i, j, k))
2914 : END IF
2915 : END DO
2916 : END DO
2917 : END DO
2918 : !$OMP END PARALLEL DO
2919 : END IF
2920 : END DO
2921 :
2922 4362 : DO ispin = 1, nspins
2923 : ! partial integration
2924 4362 : CALL xc_pw_divergence(xc_deriv_method_id, v_drho_r(:, ispin), tmp_g, vxc_g, v_xc(ispin))
2925 : END DO ! ispin
2926 :
2927 : END IF
2928 :
2929 8752 : DO idir = 1, 3
2930 6564 : DEALLOCATE (drho(idir)%array)
2931 8752 : DEALLOCATE (drho1(idir)%array)
2932 : END DO
2933 :
2934 5694 : DO ispin = 1, nspins
2935 3506 : CALL deallocate_pw(v_drhoa(ispin), pw_pool)
2936 5694 : CALL deallocate_pw(v_drhob(ispin), pw_pool)
2937 : END DO
2938 :
2939 2188 : DEALLOCATE (v_drhoa, v_drhob)
2940 :
2941 : END IF ! gradient_f
2942 :
2943 4594 : IF (laplace_f .AND. my_compute_virial) THEN
2944 15026 : virial_pw%array(:, :, :) = -rhoa(:, :, :)
2945 2 : CALL virial_laplace(virial_pw, pw_pool, virial_xc, v_laplace(1)%array)
2946 15026 : virial_pw%array(:, :, :) = -rhob(:, :, :)
2947 2 : CALL virial_laplace(virial_pw, pw_pool, virial_xc, v_laplace(2)%array)
2948 : END IF
2949 :
2950 : ELSE
2951 :
2952 : !-----------------!
2953 : ! restricted case !
2954 : !-----------------!
2955 :
2956 27176 : CALL xc_rho_set_get(rho1_set, rho=rho1)
2957 :
2958 27176 : IF (gradient_f) THEN
2959 17566 : CALL xc_rho_set_get(rho_set, drho=drho, norm_drho=norm_drho)
2960 17566 : CALL xc_rho_set_get(rho1_set, drho=drho1)
2961 17566 : CALL prepare_dr1dr(dr1dr, drho, drho1)
2962 : END IF
2963 :
2964 27176 : IF (laplace_f) THEN
2965 136 : CALL xc_rho_set_get(rho1_set, laplace_rho=laplace1)
2966 :
2967 544 : ALLOCATE (v_laplace(nspins))
2968 272 : DO ispin = 1, nspins
2969 272 : CALL allocate_pw(v_laplace(ispin), pw_pool, bo)
2970 : END DO
2971 :
2972 136 : IF (my_compute_virial) CALL xc_rho_set_get(rho_set, rho=rho)
2973 : END IF
2974 :
2975 27176 : IF (tau_f) THEN
2976 562 : CALL xc_rho_set_get(rho1_set, tau=tau1)
2977 : END IF
2978 :
2979 321692 : $:add_2nd_derivative_terms(arguments_closedshell)
2980 :
2981 27176 : IF (gradient_f) THEN
2982 :
2983 17566 : IF (my_compute_virial) THEN
2984 222 : CALL virial_drho_drho(virial_pw, drho, v_drho(1), virial_xc)
2985 : END IF ! my_compute_virial
2986 :
2987 17566 : IF (my_gapw) THEN
2988 :
2989 26640 : DO idir = 1, 3
2990 : !$OMP PARALLEL DO DEFAULT(NONE) &
2991 : !$OMP PRIVATE(ia,ir) &
2992 : !$OMP SHARED(bo,vxg,drho,v_drho,e_drho,drho1,idir,factor2) &
2993 26640 : !$OMP COLLAPSE(2)
2994 : DO ia = bo(1, 1), bo(2, 1)
2995 : DO ir = bo(1, 2), bo(2, 2)
2996 : vxg(idir, ia, ir, 1) = -drho(idir)%array(ia, ir, 1)*v_drho(1)%array(ia, ir, 1)
2997 : IF (ASSOCIATED(e_drho)) THEN
2998 : vxg(idir, ia, ir, 1) = vxg(idir, ia, ir, 1) + factor2*drho1(idir)%array(ia, ir, 1)*e_drho(ia, ir, 1)
2999 : END IF
3000 : END DO
3001 : END DO
3002 : !$OMP END PARALLEL DO
3003 : END DO
3004 :
3005 : ELSE
3006 : ! partial integration
3007 43624 : DO idir = 1, 3
3008 43624 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(v_drho_r,drho,v_drho,drho1,e_drho,idir)
3009 : v_drho_r(idir, 1)%array(:, :, :) = drho(idir)%array(:, :, :)*v_drho(1)%array(:, :, :) - &
3010 : drho1(idir)%array(:, :, :)*e_drho(:, :, :)
3011 : !$OMP END PARALLEL WORKSHARE
3012 : END DO
3013 :
3014 10906 : CALL xc_pw_divergence(xc_deriv_method_id, v_drho_r(:, 1), tmp_g, vxc_g, v_xc(1))
3015 : END IF
3016 :
3017 : END IF
3018 :
3019 27176 : IF (laplace_f .AND. my_compute_virial) THEN
3020 294530 : virial_pw%array(:, :, :) = -rho(:, :, :)
3021 14 : CALL virial_laplace(virial_pw, pw_pool, virial_xc, v_laplace(1)%array)
3022 : END IF
3023 :
3024 : END IF
3025 :
3026 31770 : IF (laplace_f) THEN
3027 386 : DO ispin = 1, nspins
3028 212 : CALL xc_pw_laplace(v_laplace(ispin), pw_pool, xc_deriv_method_id)
3029 386 : CALL pw_axpy(v_laplace(ispin), v_xc(ispin))
3030 : END DO
3031 : END IF
3032 :
3033 31770 : IF (gradient_f) THEN
3034 :
3035 40826 : DO ispin = 1, nspins
3036 21072 : CALL deallocate_pw(v_drho(ispin), pw_pool)
3037 104042 : DO idir = 1, 3
3038 84288 : CALL deallocate_pw(v_drho_r(idir, ispin), pw_pool)
3039 : END DO
3040 : END DO
3041 19754 : DEALLOCATE (v_drho, v_drho_r)
3042 :
3043 : END IF
3044 :
3045 31770 : IF (laplace_f) THEN
3046 386 : DO ispin = 1, nspins
3047 386 : CALL deallocate_pw(v_laplace(ispin), pw_pool)
3048 : END DO
3049 174 : DEALLOCATE (v_laplace)
3050 : END IF
3051 :
3052 31770 : IF (ASSOCIATED(tmp_g%pw_grid) .AND. ASSOCIATED(pw_pool)) THEN
3053 12060 : CALL pw_pool%give_back_pw(tmp_g)
3054 : END IF
3055 :
3056 31770 : IF (ASSOCIATED(vxc_g%pw_grid) .AND. ASSOCIATED(pw_pool)) THEN
3057 12060 : CALL pw_pool%give_back_pw(vxc_g)
3058 : END IF
3059 :
3060 31770 : IF (my_compute_virial .AND. (gradient_f .OR. laplace_f)) THEN
3061 232 : CALL deallocate_pw(virial_pw, pw_pool)
3062 : END IF
3063 :
3064 31770 : CALL timestop(handle)
3065 :
3066 95310 : END SUBROUTINE xc_calc_2nd_deriv_analytical
3067 :
3068 : ! **************************************************************************************************
3069 : !> \brief allocates grids using pw_pool (if associated) or with bounds
3070 : !> \param pw ...
3071 : !> \param pw_pool ...
3072 : !> \param bo ...
3073 : ! **************************************************************************************************
3074 91842 : SUBROUTINE allocate_pw(pw, pw_pool, bo)
3075 : TYPE(pw_r3d_rs_type), INTENT(OUT) :: pw
3076 : TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool
3077 : INTEGER, DIMENSION(2, 3), INTENT(IN) :: bo
3078 :
3079 91842 : IF (ASSOCIATED(pw_pool)) THEN
3080 60822 : CALL pw_pool%create_pw(pw)
3081 60822 : CALL pw_zero(pw)
3082 : ELSE
3083 155100 : ALLOCATE (pw%array(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
3084 79163040 : pw%array = 0.0_dp
3085 : END IF
3086 :
3087 91842 : END SUBROUTINE allocate_pw
3088 :
3089 : ! **************************************************************************************************
3090 : !> \brief deallocates grid allocated with allocate_pw
3091 : !> \param pw ...
3092 : !> \param pw_pool ...
3093 : ! **************************************************************************************************
3094 91842 : SUBROUTINE deallocate_pw(pw, pw_pool)
3095 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: pw
3096 : TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool
3097 :
3098 91842 : IF (ASSOCIATED(pw_pool)) THEN
3099 60822 : CALL pw_pool%give_back_pw(pw)
3100 : ELSE
3101 31020 : CALL pw%release()
3102 : END IF
3103 :
3104 91842 : END SUBROUTINE deallocate_pw
3105 :
3106 : ! **************************************************************************************************
3107 : !> \brief updates virial from first derivative w.r.t. norm_drho
3108 : !> \param virial_pw ...
3109 : !> \param drho ...
3110 : !> \param drho1 ...
3111 : !> \param deriv_data ...
3112 : !> \param virial_xc ...
3113 : ! **************************************************************************************************
3114 304 : SUBROUTINE virial_drho_drho1(virial_pw, drho, drho1, deriv_data, virial_xc)
3115 : TYPE(pw_r3d_rs_type), INTENT(IN) :: virial_pw
3116 : TYPE(cp_3d_r_cp_type), DIMENSION(3), INTENT(IN) :: drho, drho1
3117 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: deriv_data
3118 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT) :: virial_xc
3119 :
3120 : INTEGER :: idir, jdir
3121 : REAL(KIND=dp) :: tmp
3122 :
3123 1216 : DO idir = 1, 3
3124 912 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(drho,idir,virial_pw,deriv_data)
3125 : virial_pw%array(:, :, :) = drho(idir)%array(:, :, :)*deriv_data(:, :, :)
3126 : !$OMP END PARALLEL WORKSHARE
3127 3952 : DO jdir = 1, 3
3128 : tmp = virial_pw%pw_grid%dvol*accurate_dot_product(virial_pw%array(:, :, :), &
3129 2736 : drho1(jdir)%array(:, :, :))
3130 2736 : virial_xc(jdir, idir) = virial_xc(jdir, idir) + tmp
3131 3648 : virial_xc(idir, jdir) = virial_xc(idir, jdir) + tmp
3132 : END DO
3133 : END DO
3134 :
3135 304 : END SUBROUTINE virial_drho_drho1
3136 :
3137 : ! **************************************************************************************************
3138 : !> \brief Adds virial contribution from second order potential parts
3139 : !> \param virial_pw ...
3140 : !> \param drho ...
3141 : !> \param v_drho ...
3142 : !> \param virial_xc ...
3143 : ! **************************************************************************************************
3144 298 : SUBROUTINE virial_drho_drho(virial_pw, drho, v_drho, virial_xc)
3145 : TYPE(pw_r3d_rs_type), INTENT(IN) :: virial_pw
3146 : TYPE(cp_3d_r_cp_type), DIMENSION(3), INTENT(IN) :: drho
3147 : TYPE(pw_r3d_rs_type), INTENT(IN) :: v_drho
3148 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT) :: virial_xc
3149 :
3150 : INTEGER :: idir, jdir
3151 : REAL(KIND=dp) :: tmp
3152 :
3153 1192 : DO idir = 1, 3
3154 894 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(drho,idir,v_drho,virial_pw)
3155 : virial_pw%array(:, :, :) = drho(idir)%array(:, :, :)*v_drho%array(:, :, :)
3156 : !$OMP END PARALLEL WORKSHARE
3157 2980 : DO jdir = 1, idir
3158 : tmp = -virial_pw%pw_grid%dvol*accurate_dot_product(virial_pw%array(:, :, :), &
3159 1788 : drho(jdir)%array(:, :, :))
3160 1788 : virial_xc(jdir, idir) = virial_xc(jdir, idir) + tmp
3161 2682 : virial_xc(idir, jdir) = virial_xc(jdir, idir)
3162 : END DO
3163 : END DO
3164 :
3165 298 : END SUBROUTINE virial_drho_drho
3166 :
3167 : ! **************************************************************************************************
3168 : !> \brief ...
3169 : !> \param rho_r ...
3170 : !> \param pw_pool ...
3171 : !> \param virial_xc ...
3172 : !> \param deriv_data ...
3173 : ! **************************************************************************************************
3174 150 : SUBROUTINE virial_laplace(rho_r, pw_pool, virial_xc, deriv_data)
3175 : TYPE(pw_r3d_rs_type), TARGET :: rho_r
3176 : TYPE(pw_pool_type), POINTER, INTENT(IN) :: pw_pool
3177 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT) :: virial_xc
3178 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: deriv_data
3179 :
3180 : CHARACTER(len=*), PARAMETER :: routineN = 'virial_laplace'
3181 :
3182 : INTEGER :: handle, idir, jdir
3183 : TYPE(pw_r3d_rs_type), POINTER :: virial_pw
3184 : TYPE(pw_c1d_gs_type), POINTER :: tmp_g, rho_g
3185 : INTEGER, DIMENSION(3) :: my_deriv
3186 :
3187 150 : CALL timeset(routineN, handle)
3188 :
3189 150 : NULLIFY (virial_pw, tmp_g, rho_g)
3190 150 : ALLOCATE (virial_pw, tmp_g, rho_g)
3191 150 : CALL pw_pool%create_pw(virial_pw)
3192 150 : CALL pw_pool%create_pw(tmp_g)
3193 150 : CALL pw_pool%create_pw(rho_g)
3194 150 : CALL pw_zero(virial_pw)
3195 150 : CALL pw_transfer(rho_r, rho_g)
3196 600 : DO idir = 1, 3
3197 1500 : DO jdir = idir, 3
3198 900 : CALL pw_copy(rho_g, tmp_g)
3199 :
3200 900 : my_deriv = 0
3201 900 : my_deriv(idir) = 1
3202 900 : my_deriv(jdir) = my_deriv(jdir) + 1
3203 :
3204 900 : CALL pw_derive(tmp_g, my_deriv)
3205 900 : CALL pw_transfer(tmp_g, virial_pw)
3206 : virial_xc(idir, jdir) = virial_xc(idir, jdir) - 2.0_dp*virial_pw%pw_grid%dvol* &
3207 : accurate_dot_product(virial_pw%array(:, :, :), &
3208 900 : deriv_data(:, :, :))
3209 1350 : virial_xc(jdir, idir) = virial_xc(idir, jdir)
3210 : END DO
3211 : END DO
3212 150 : CALL pw_pool%give_back_pw(virial_pw)
3213 150 : CALL pw_pool%give_back_pw(tmp_g)
3214 150 : CALL pw_pool%give_back_pw(rho_g)
3215 150 : DEALLOCATE (virial_pw, tmp_g, rho_g)
3216 :
3217 150 : CALL timestop(handle)
3218 :
3219 150 : END SUBROUTINE virial_laplace
3220 :
3221 : ! **************************************************************************************************
3222 : !> \brief Prepare objects for the calculation of the 2nd derivatives of the density functional.
3223 : !> The calculation must then be performed with xc_calc_2nd_deriv.
3224 : !> \param deriv_set object containing the XC derivatives (out)
3225 : !> \param rho_set object that will contain the density at which the
3226 : !> derivatives were calculated
3227 : !> \param rho_r the place where you evaluate the derivative
3228 : !> \param pw_pool the pool for the grids
3229 : !> \param xc_section which functional should be used and how to calculate it
3230 : !> \param tau_r kinetic energy density in real space
3231 : ! **************************************************************************************************
3232 12092 : SUBROUTINE xc_prep_2nd_deriv(deriv_set, &
3233 : rho_set, rho_r, pw_pool, xc_section, tau_r)
3234 :
3235 : TYPE(xc_derivative_set_type) :: deriv_set
3236 : TYPE(xc_rho_set_type) :: rho_set
3237 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
3238 : TYPE(pw_pool_type), POINTER :: pw_pool
3239 : TYPE(section_vals_type), POINTER :: xc_section
3240 : TYPE(pw_r3d_rs_type), DIMENSION(:), OPTIONAL, POINTER :: tau_r
3241 :
3242 : CHARACTER(len=*), PARAMETER :: routineN = 'xc_prep_2nd_deriv'
3243 :
3244 : INTEGER :: handle, nspins
3245 : LOGICAL :: lsd
3246 12092 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
3247 12092 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: tau
3248 :
3249 12092 : CALL timeset(routineN, handle)
3250 :
3251 12092 : CPASSERT(ASSOCIATED(xc_section))
3252 12092 : CPASSERT(ASSOCIATED(pw_pool))
3253 :
3254 12092 : nspins = SIZE(rho_r)
3255 12092 : lsd = (nspins /= 1)
3256 :
3257 12092 : NULLIFY (rho_g, tau)
3258 12092 : IF (PRESENT(tau_r)) &
3259 11430 : tau => tau_r
3260 :
3261 12092 : IF (section_get_lval(xc_section, "2ND_DERIV_ANALYTICAL")) THEN
3262 : CALL xc_rho_set_and_dset_create(rho_set, deriv_set, 2, &
3263 : rho_r, rho_g, tau, xc_section, pw_pool, &
3264 12008 : calc_potential=.TRUE.)
3265 : ELSE
3266 : CALL xc_rho_set_and_dset_create(rho_set, deriv_set, 1, &
3267 : rho_r, rho_g, tau, xc_section, pw_pool, &
3268 84 : calc_potential=.TRUE.)
3269 : END IF
3270 :
3271 12092 : CALL timestop(handle)
3272 :
3273 12092 : END SUBROUTINE xc_prep_2nd_deriv
3274 :
3275 : ! **************************************************************************************************
3276 : !> \brief divides derivatives from deriv_set by norm_drho
3277 : !> \param deriv_set ...
3278 : !> \param rho_set ...
3279 : !> \param lsd ...
3280 : ! **************************************************************************************************
3281 140505 : SUBROUTINE divide_by_norm_drho(deriv_set, rho_set, lsd)
3282 :
3283 : TYPE(xc_derivative_set_type), INTENT(INOUT) :: deriv_set
3284 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
3285 : LOGICAL, INTENT(IN) :: lsd
3286 :
3287 140505 : INTEGER, DIMENSION(:), POINTER :: split_desc
3288 : INTEGER :: idesc
3289 : INTEGER, DIMENSION(2, 3) :: bo
3290 : REAL(KIND=dp) :: drho_cutoff
3291 140505 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: norm_drho, norm_drhoa, norm_drhob
3292 1686060 : TYPE(cp_3d_r_cp_type), DIMENSION(3) :: drho, drhoa, drhob
3293 : TYPE(cp_sll_xc_deriv_type), POINTER :: pos
3294 : TYPE(xc_derivative_type), POINTER :: deriv_att
3295 :
3296 : ! check for unknown derivatives and divide by norm_drho where necessary
3297 :
3298 1405050 : bo = rho_set%local_bounds
3299 : CALL xc_rho_set_get(rho_set, drho_cutoff=drho_cutoff, norm_drho=norm_drho, &
3300 : norm_drhoa=norm_drhoa, norm_drhob=norm_drhob, &
3301 140505 : drho=drho, drhoa=drhoa, drhob=drhob, can_return_null=.TRUE.)
3302 :
3303 140505 : pos => deriv_set%derivs
3304 602726 : DO WHILE (cp_sll_xc_deriv_next(pos, el_att=deriv_att))
3305 462221 : CALL xc_derivative_get(deriv_att, split_desc=split_desc)
3306 986138 : DO idesc = 1, SIZE(split_desc)
3307 462221 : SELECT CASE (split_desc(idesc))
3308 : CASE (deriv_norm_drho)
3309 115483 : IF (ASSOCIATED(norm_drho)) THEN
3310 115483 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(deriv_att,norm_drho,drho_cutoff)
3311 : deriv_att%deriv_data(:, :, :) = deriv_att%deriv_data(:, :, :)/ &
3312 : MAX(norm_drho(:, :, :), drho_cutoff)
3313 : !$OMP END PARALLEL WORKSHARE
3314 0 : ELSE IF (ASSOCIATED(drho(1)%array)) THEN
3315 0 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(deriv_att,drho,drho_cutoff)
3316 : deriv_att%deriv_data(:, :, :) = deriv_att%deriv_data(:, :, :)/ &
3317 : MAX(SQRT(drho(1)%array(:, :, :)**2 + &
3318 : drho(2)%array(:, :, :)**2 + &
3319 : drho(3)%array(:, :, :)**2), drho_cutoff)
3320 : !$OMP END PARALLEL WORKSHARE
3321 0 : ELSE IF (ASSOCIATED(drhoa(1)%array) .AND. ASSOCIATED(drhob(1)%array)) THEN
3322 0 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(deriv_att,drhoa,drhob,drho_cutoff)
3323 : deriv_att%deriv_data(:, :, :) = deriv_att%deriv_data(:, :, :)/ &
3324 : MAX(SQRT((drhoa(1)%array(:, :, :) + drhob(1)%array(:, :, :))**2 + &
3325 : (drhoa(2)%array(:, :, :) + drhob(2)%array(:, :, :))**2 + &
3326 : (drhoa(3)%array(:, :, :) + drhob(3)%array(:, :, :))**2), drho_cutoff)
3327 : !$OMP END PARALLEL WORKSHARE
3328 : ELSE
3329 0 : CPABORT("Normalization of derivative requires any of norm_drho, drho or drhoa+drhob!")
3330 : END IF
3331 : CASE (deriv_norm_drhoa)
3332 19440 : IF (ASSOCIATED(norm_drhoa)) THEN
3333 19440 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(deriv_att,norm_drhoa,drho_cutoff)
3334 : deriv_att%deriv_data(:, :, :) = deriv_att%deriv_data(:, :, :)/ &
3335 : MAX(norm_drhoa(:, :, :), drho_cutoff)
3336 : !$OMP END PARALLEL WORKSHARE
3337 0 : ELSE IF (ASSOCIATED(drhoa(1)%array)) THEN
3338 0 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(deriv_att,drhoa,drho_cutoff)
3339 : deriv_att%deriv_data(:, :, :) = deriv_att%deriv_data(:, :, :)/ &
3340 : MAX(SQRT(drhoa(1)%array(:, :, :)**2 + &
3341 : drhoa(2)%array(:, :, :)**2 + &
3342 : drhoa(3)%array(:, :, :)**2), drho_cutoff)
3343 : !$OMP END PARALLEL WORKSHARE
3344 : ELSE
3345 0 : CPABORT("Normalization of derivative requires any of norm_drhoa or drhoa!")
3346 : END IF
3347 : CASE (deriv_norm_drhob)
3348 19436 : IF (ASSOCIATED(norm_drhob)) THEN
3349 19436 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(deriv_att,norm_drhob,drho_cutoff)
3350 : deriv_att%deriv_data(:, :, :) = deriv_att%deriv_data(:, :, :)/ &
3351 : MAX(norm_drhob(:, :, :), drho_cutoff)
3352 : !$OMP END PARALLEL WORKSHARE
3353 0 : ELSE IF (ASSOCIATED(drhob(1)%array)) THEN
3354 0 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(deriv_att,drhob,drho_cutoff)
3355 : deriv_att%deriv_data(:, :, :) = deriv_att%deriv_data(:, :, :)/ &
3356 : MAX(SQRT(drhob(1)%array(:, :, :)**2 + &
3357 : drhob(2)%array(:, :, :)**2 + &
3358 : drhob(3)%array(:, :, :)**2), drho_cutoff)
3359 : !$OMP END PARALLEL WORKSHARE
3360 : ELSE
3361 0 : CPABORT("Normalization of derivative requires any of norm_drhob or drhob!")
3362 : END IF
3363 : CASE (deriv_rho, deriv_tau, deriv_laplace_rho)
3364 160849 : IF (lsd) &
3365 0 : CPABORT(TRIM(id_to_desc(split_desc(idesc)))//" not handled in lsd!'")
3366 : CASE (deriv_rhoa, deriv_rhob, deriv_tau_a, deriv_tau_b, deriv_laplace_rhoa, deriv_laplace_rhob)
3367 : CASE default
3368 383412 : CPABORT("Unknown derivative id")
3369 : END SELECT
3370 : END DO
3371 : END DO
3372 :
3373 140505 : END SUBROUTINE divide_by_norm_drho
3374 :
3375 : ! **************************************************************************************************
3376 : !> \brief allocates and calculates drho from given spin densities drhoa, drhob
3377 : !> \param drho ...
3378 : !> \param drhoa ...
3379 : !> \param drhob ...
3380 : ! **************************************************************************************************
3381 17584 : SUBROUTINE calc_drho_from_ab(drho, drhoa, drhob)
3382 : TYPE(cp_3d_r_cp_type), DIMENSION(3), INTENT(OUT) :: drho
3383 : TYPE(cp_3d_r_cp_type), DIMENSION(3), INTENT(IN) :: drhoa, drhob
3384 :
3385 : CHARACTER(len=*), PARAMETER :: routineN = 'calc_drho_from_ab'
3386 :
3387 : INTEGER :: handle, idir
3388 :
3389 4396 : CALL timeset(routineN, handle)
3390 :
3391 17584 : DO idir = 1, 3
3392 13188 : NULLIFY (drho(idir)%array)
3393 : ALLOCATE (drho(idir)%array(LBOUND(drhoa(1)%array, 1):UBOUND(drhoa(1)%array, 1), &
3394 : LBOUND(drhoa(1)%array, 2):UBOUND(drhoa(1)%array, 2), &
3395 158256 : LBOUND(drhoa(1)%array, 3):UBOUND(drhoa(1)%array, 3)))
3396 17584 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(drho,drhoa,drhob,idir)
3397 : drho(idir)%array(:, :, :) = drhoa(idir)%array(:, :, :) + drhob(idir)%array(:, :, :)
3398 : !$OMP END PARALLEL WORKSHARE
3399 : END DO
3400 :
3401 4396 : CALL timestop(handle)
3402 :
3403 4396 : END SUBROUTINE calc_drho_from_ab
3404 :
3405 : ! **************************************************************************************************
3406 : !> \brief allocates and calculates dot products of two density gradients
3407 : !> \param dr1dr ...
3408 : !> \param drho ...
3409 : !> \param drho1 ...
3410 : ! **************************************************************************************************
3411 23316 : SUBROUTINE prepare_dr1dr(dr1dr, drho, drho1)
3412 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
3413 : INTENT(OUT) :: dr1dr
3414 : TYPE(cp_3d_r_cp_type), DIMENSION(3), INTENT(IN) :: drho, drho1
3415 :
3416 : CHARACTER(len=*), PARAMETER :: routineN = 'prepare_dr1dr'
3417 :
3418 : INTEGER :: handle, idir
3419 :
3420 23316 : CALL timeset(routineN, handle)
3421 :
3422 0 : ALLOCATE (dr1dr(LBOUND(drho(1)%array, 1):UBOUND(drho(1)%array, 1), &
3423 : LBOUND(drho(1)%array, 2):UBOUND(drho(1)%array, 2), &
3424 279792 : LBOUND(drho(1)%array, 3):UBOUND(drho(1)%array, 3)))
3425 :
3426 23316 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(dr1dr,drho,drho1)
3427 : dr1dr(:, :, :) = drho(1)%array(:, :, :)*drho1(1)%array(:, :, :)
3428 : !$OMP END PARALLEL WORKSHARE
3429 69948 : DO idir = 2, 3
3430 69948 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(dr1dr,drho,drho1,idir)
3431 : dr1dr(:, :, :) = dr1dr(:, :, :) + drho(idir)%array(:, :, :)*drho1(idir)%array(:, :, :)
3432 : !$OMP END PARALLEL WORKSHARE
3433 : END DO
3434 :
3435 23316 : CALL timestop(handle)
3436 :
3437 23316 : END SUBROUTINE prepare_dr1dr
3438 :
3439 : ! **************************************************************************************************
3440 : !> \brief allocates and calculates dot product of two densities for triplets
3441 : !> \param dr1dr ...
3442 : !> \param drhoa ...
3443 : !> \param drhob ...
3444 : !> \param drho1a ...
3445 : !> \param drho1b ...
3446 : !> \param fac ...
3447 : ! **************************************************************************************************
3448 870 : SUBROUTINE prepare_dr1dr_ab(dr1dr, drhoa, drhob, drho1a, drho1b, fac)
3449 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
3450 : INTENT(OUT) :: dr1dr
3451 : TYPE(cp_3d_r_cp_type), DIMENSION(3), INTENT(IN) :: drhoa, drhob, drho1a, drho1b
3452 : REAL(KIND=dp), INTENT(IN) :: fac
3453 :
3454 : CHARACTER(len=*), PARAMETER :: routineN = 'prepare_dr1dr_ab'
3455 :
3456 : INTEGER :: handle, idir
3457 :
3458 870 : CALL timeset(routineN, handle)
3459 :
3460 0 : ALLOCATE (dr1dr(LBOUND(drhoa(1)%array, 1):UBOUND(drhoa(1)%array, 1), &
3461 : LBOUND(drhoa(1)%array, 2):UBOUND(drhoa(1)%array, 2), &
3462 10440 : LBOUND(drhoa(1)%array, 3):UBOUND(drhoa(1)%array, 3)))
3463 :
3464 870 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(fac,dr1dr,drho1a,drho1b,drhoa,drhob)
3465 : dr1dr(:, :, :) = drhoa(1)%array(:, :, :)*(drho1a(1)%array(:, :, :) + &
3466 : fac*drho1b(1)%array(:, :, :)) + &
3467 : drhob(1)%array(:, :, :)*(fac*drho1a(1)%array(:, :, :) + &
3468 : drho1b(1)%array(:, :, :))
3469 : !$OMP END PARALLEL WORKSHARE
3470 2610 : DO idir = 2, 3
3471 2610 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(fac,dr1dr,drho1a,drho1b,drhoa,drhob,idir)
3472 : dr1dr(:, :, :) = dr1dr(:, :, :) + &
3473 : drhoa(idir)%array(:, :, :)*(drho1a(idir)%array(:, :, :) + &
3474 : fac*drho1b(idir)%array(:, :, :)) + &
3475 : drhob(idir)%array(:, :, :)*(fac*drho1a(idir)%array(:, :, :) + &
3476 : drho1b(idir)%array(:, :, :))
3477 : !$OMP END PARALLEL WORKSHARE
3478 : END DO
3479 :
3480 870 : CALL timestop(handle)
3481 :
3482 870 : END SUBROUTINE prepare_dr1dr_ab
3483 :
3484 : ! **************************************************************************************************
3485 : !> \brief checks for gradients
3486 : !> \param deriv_set ...
3487 : !> \param lsd ...
3488 : !> \param gradient_f ...
3489 : !> \param tau_f ...
3490 : !> \param laplace_f ...
3491 : ! **************************************************************************************************
3492 142046 : SUBROUTINE check_for_derivatives(deriv_set, lsd, rho_f, gradient_f, tau_f, laplace_f)
3493 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
3494 : LOGICAL, INTENT(IN) :: lsd
3495 : LOGICAL, INTENT(OUT) :: rho_f, gradient_f, tau_f, laplace_f
3496 :
3497 : CHARACTER(len=*), PARAMETER :: routineN = 'check_for_derivatives'
3498 :
3499 : INTEGER :: handle, iorder, order
3500 142046 : INTEGER, DIMENSION(:), POINTER :: split_desc
3501 : TYPE(cp_sll_xc_deriv_type), POINTER :: pos
3502 : TYPE(xc_derivative_type), POINTER :: deriv_att
3503 :
3504 142046 : CALL timeset(routineN, handle)
3505 :
3506 142046 : rho_f = .FALSE.
3507 142046 : gradient_f = .FALSE.
3508 142046 : tau_f = .FALSE.
3509 142046 : laplace_f = .FALSE.
3510 : ! check for unknown derivatives
3511 142046 : pos => deriv_set%derivs
3512 647522 : DO WHILE (cp_sll_xc_deriv_next(pos, el_att=deriv_att))
3513 : CALL xc_derivative_get(deriv_att, order=order, &
3514 505476 : split_desc=split_desc)
3515 647522 : IF (lsd) THEN
3516 309791 : DO iorder = 1, size(split_desc)
3517 153871 : SELECT CASE (split_desc(iorder))
3518 : CASE (deriv_rhoa, deriv_rhob)
3519 80966 : rho_f = .TRUE.
3520 : CASE (deriv_norm_drho, deriv_norm_drhoa, deriv_norm_drhob)
3521 70994 : gradient_f = .TRUE.
3522 : CASE (deriv_tau_a, deriv_tau_b)
3523 2704 : tau_f = .TRUE.
3524 : CASE (deriv_laplace_rhoa, deriv_laplace_rhob)
3525 1256 : laplace_f = .TRUE.
3526 : CASE (deriv_rho, deriv_tau, deriv_laplace_rho)
3527 0 : CPABORT("Derivative not handled in lsd!")
3528 : CASE default
3529 155920 : CPABORT("Unknown derivative id")
3530 : END SELECT
3531 : END DO
3532 : ELSE
3533 652733 : DO iorder = 1, size(split_desc)
3534 351605 : SELECT CASE (split_desc(iorder))
3535 : CASE (deriv_rho)
3536 178505 : rho_f = .TRUE.
3537 : CASE (deriv_tau)
3538 4798 : tau_f = .TRUE.
3539 : CASE (deriv_norm_drho)
3540 116437 : gradient_f = .TRUE.
3541 : CASE (deriv_laplace_rho)
3542 1388 : laplace_f = .TRUE.
3543 : CASE default
3544 301128 : CPABORT("Unknown derivative id")
3545 : END SELECT
3546 : END DO
3547 : END IF
3548 : END DO
3549 :
3550 142046 : CALL timestop(handle)
3551 :
3552 142046 : END SUBROUTINE check_for_derivatives
3553 :
3554 : END MODULE xc
|
