Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief 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 7594 : 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 15188 : calc_potential=.FALSE.)
101 7594 : res = (needs%tau_spin .OR. needs%tau)
102 :
103 7594 : END FUNCTION
104 :
105 : ! **************************************************************************************************
106 : !> \brief ...
107 : !> \param xc_fun_section ...
108 : !> \param lsd ...
109 : !> \return ...
110 : ! **************************************************************************************************
111 7610 : 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 15220 : calc_potential=.FALSE.)
121 7610 : res = (needs%norm_drho .OR. needs%norm_drho_spin)
122 :
123 7610 : 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 217810 : 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 108905 : CPASSERT(ASSOCIATED(rho_r))
164 108905 : CPASSERT(ASSOCIATED(xc_section))
165 108905 : CPASSERT(.NOT. ASSOCIATED(vxc_rho))
166 108905 : CPASSERT(.NOT. ASSOCIATED(vxc_tau))
167 :
168 108905 : virial_xc = 0.0_dp
169 :
170 : CALL section_vals_val_get(xc_section, "FUNCTIONAL_ROUTINE", &
171 108905 : i_val=f_routine)
172 108905 : 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 108905 : 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 108905 : xc_section=xc_section, pw_pool=pw_pool)
187 : CASE default
188 : END SELECT
189 :
190 108905 : 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 266942 : 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 133471 : CALL timeset(routineN, handle)
859 :
860 133471 : CPASSERT(ASSOCIATED(pw_pool))
861 :
862 133471 : nspins = SIZE(rho_r)
863 133471 : lsd = (nspins /= 1)
864 :
865 133471 : xc_fun_sections => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
866 :
867 133471 : 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 133471 : 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 133471 : 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 133471 : deriv_order=deriv_order)
886 :
887 133471 : CALL divide_by_norm_drho(deriv_set, rho_set, lsd)
888 :
889 133471 : CALL timestop(handle)
890 :
891 133471 : 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 253967 : 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 253967 : CPASSERT(ASSOCIATED(pot))
924 1015868 : bo(1, :) = LBOUND(pot)
925 1015868 : bo(2, :) = UBOUND(pot)
926 253967 : my_e_0_scale_factor = 1.0_dp
927 253967 : IF (PRESENT(e_0_scale_factor)) my_e_0_scale_factor = e_0_scale_factor
928 253967 : rho_smooth_cutoff = rho_cutoff*rho_smooth_cutoff_range
929 253967 : rho_smooth_cutoff_2 = (rho_cutoff + rho_smooth_cutoff)/2
930 253967 : rho_smooth_cutoff_range_2 = rho_smooth_cutoff_2 - rho_cutoff
931 :
932 253967 : 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 253967 : 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 4139130 : pot%array(i, j, k) = my_pot/my_rho
1097 229038 : ELSE IF (my_rho < eps2) THEN
1098 229038 : 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 110873 : 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 110873 : REAL(kind=dp), DIMENSION(:, :, :), POINTER :: deriv_data, norm_drho, norm_drho_spin, &
1165 221746 : rho, rhoa, rhob
1166 : TYPE(cp_sll_xc_deriv_type), POINTER :: pos
1167 : TYPE(pw_grid_type), POINTER :: pw_grid
1168 776111 : 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 110873 : CALL timeset(routineN, handle)
1176 110873 : NULLIFY (norm_drho_spin, norm_drho, pos)
1177 :
1178 110873 : pw_grid => rho_r(1)%pw_grid
1179 :
1180 110873 : CPASSERT(ASSOCIATED(xc_section))
1181 110873 : CPASSERT(ASSOCIATED(pw_pool))
1182 110873 : CPASSERT(.NOT. ASSOCIATED(vxc_rho))
1183 110873 : CPASSERT(.NOT. ASSOCIATED(vxc_tau))
1184 110873 : nspins = SIZE(rho_r)
1185 110873 : lsd = (nspins /= 1)
1186 110873 : IF (lsd) THEN
1187 22053 : CPASSERT(nspins == 2)
1188 : END IF
1189 :
1190 110873 : use_virial = compute_virial
1191 110873 : virial_xc = 0.0_dp
1192 :
1193 1108730 : bo = rho_r(1)%pw_grid%bounds_local
1194 110873 : 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 110873 : calc_potential=.TRUE.)
1202 :
1203 : CALL section_vals_val_get(xc_section, "XC_GRID%XC_DERIV", &
1204 110873 : i_val=xc_deriv_method_id)
1205 : CALL section_vals_val_get(xc_section, "XC_GRID%XC_SMOOTH_RHO", &
1206 110873 : i_val=xc_rho_smooth_id)
1207 : CALL section_vals_val_get(xc_section, "DENSITY_SMOOTH_CUTOFF_RANGE", &
1208 110873 : r_val=density_smooth_cut_range)
1209 :
1210 : CALL xc_rho_set_get(rho_set, rho_cutoff=rho_cutoff, &
1211 110873 : drho_cutoff=drho_cutoff)
1212 :
1213 110873 : CALL check_for_derivatives(deriv_set, lsd, has_rho, has_gradient, has_tau, has_laplace)
1214 : ! check for unknown derivatives
1215 110873 : has_derivs = has_rho .OR. has_gradient .OR. has_tau .OR. has_laplace
1216 :
1217 465545 : ALLOCATE (vxc_rho(nspins))
1218 :
1219 : CALL xc_rho_set_get(rho_set, rho=rho, rhoa=rhoa, rhob=rhob, &
1220 110873 : can_return_null=.TRUE.)
1221 :
1222 : ! recover the vxc arrays
1223 110873 : IF (lsd) THEN
1224 22053 : CALL xc_dset_recover_pw(deriv_set, [deriv_rhoa], vxc_rho(1), pw_grid, pw_pool)
1225 22053 : CALL xc_dset_recover_pw(deriv_set, [deriv_rhob], vxc_rho(2), pw_grid, pw_pool)
1226 :
1227 : ELSE
1228 88820 : CALL xc_dset_recover_pw(deriv_set, [deriv_rho], vxc_rho(1), pw_grid, pw_pool)
1229 : END IF
1230 :
1231 110873 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_norm_drho])
1232 110873 : IF (ASSOCIATED(deriv_att)) THEN
1233 59759 : 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 59759 : can_return_null=.TRUE.)
1239 59759 : CALL xc_rho_set_recover_pw(rho_set, pw_grid, pw_pool, dealloc_pw_to_deriv_rho, drho=pw_to_deriv_rho)
1240 :
1241 59759 : CPASSERT(ASSOCIATED(deriv_data))
1242 59759 : IF (use_virial) THEN
1243 1566 : CALL pw_pool%create_pw(virial_pw)
1244 1566 : CALL pw_zero(virial_pw)
1245 6264 : DO idir = 1, 3
1246 4698 : !$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 15660 : DO jdir = 1, idir
1250 : virial_xc(idir, jdir) = -pw_grid%dvol* &
1251 : accurate_dot_product(virial_pw%array(:, :, :), &
1252 9396 : pw_to_deriv_rho(jdir)%array(:, :, :))
1253 14094 : virial_xc(jdir, idir) = virial_xc(idir, jdir)
1254 : END DO
1255 : END DO
1256 1566 : CALL pw_pool%give_back_pw(virial_pw)
1257 : END IF ! use_virial
1258 239036 : DO idir = 1, 3
1259 179277 : CPASSERT(ASSOCIATED(pw_to_deriv_rho(idir)%array))
1260 239036 : !$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 59759 : CALL pw_pool%give_back_cr3d(deriv_att%deriv_data)
1267 :
1268 : END IF
1269 :
1270 110873 : IF ((has_gradient .AND. xc_requires_tmp_g(xc_deriv_method_id)) .OR. pw_grid%spherical) THEN
1271 59529 : CALL pw_pool%create_pw(vxc_g)
1272 59529 : IF (.NOT. pw_grid%spherical) THEN
1273 59529 : CALL pw_pool%create_pw(tmp_g)
1274 : END IF
1275 : END IF
1276 :
1277 243799 : DO ispin = 1, nspins
1278 :
1279 132926 : IF (lsd) THEN
1280 44106 : IF (ispin == 1) THEN
1281 : CALL xc_rho_set_get(rho_set, norm_drhoa=norm_drho_spin, &
1282 22053 : can_return_null=.TRUE.)
1283 22053 : IF (ASSOCIATED(norm_drho_spin)) CALL xc_rho_set_recover_pw( &
1284 14008 : 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 22053 : can_return_null=.TRUE.)
1288 22053 : IF (ASSOCIATED(norm_drho_spin)) CALL xc_rho_set_recover_pw( &
1289 14008 : rho_set, pw_grid, pw_pool, dealloc_pw_to_deriv, drhob=pw_to_deriv)
1290 : END IF
1291 :
1292 88212 : deriv_att => xc_dset_get_derivative(deriv_set, [norm_drho_spin_name(ispin)])
1293 44106 : IF (ASSOCIATED(deriv_att)) THEN
1294 : CPASSERT(lsd)
1295 28016 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
1296 :
1297 28016 : 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 112064 : DO idir = 1, 3
1315 112064 : !$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 132926 : IF (ASSOCIATED(pw_to_deriv_rho(1)%array)) THEN
1324 73009 : IF (.NOT. ASSOCIATED(pw_to_deriv(1)%array)) THEN
1325 46509 : pw_to_deriv = pw_to_deriv_rho
1326 46509 : dealloc_pw_to_deriv = ((.NOT. lsd) .OR. (ispin == 2))
1327 46509 : 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 106000 : DO idir = 1, 3
1332 79500 : CALL pw_axpy(pw_to_deriv_rho(idir), pw_to_deriv(idir))
1333 106000 : IF (ispin == 2) THEN
1334 39750 : IF (dealloc_pw_to_deriv_rho) THEN
1335 39750 : 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 132926 : IF (ASSOCIATED(pw_to_deriv(1)%array)) THEN
1343 298100 : DO idir = 1, 3
1344 298100 : CALL pw_scale(pw_to_deriv(idir), -1.0_dp)
1345 : END DO
1346 :
1347 74525 : CALL xc_pw_divergence(xc_deriv_method_id, pw_to_deriv, tmp_g, vxc_g, vxc_rho(ispin))
1348 :
1349 74525 : IF (dealloc_pw_to_deriv) THEN
1350 298100 : DO idir = 1, 3
1351 298100 : 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 132926 : IF (has_laplace) THEN
1358 1092 : IF (lsd) THEN
1359 472 : IF (ispin == 1) THEN
1360 : deriv_id = deriv_laplace_rhoa
1361 : ELSE
1362 236 : deriv_id = deriv_laplace_rhob
1363 : END IF
1364 : ELSE
1365 : deriv_id = deriv_laplace_rho
1366 : END IF
1367 :
1368 2184 : CALL xc_dset_recover_pw(deriv_set, [deriv_id], pw_to_deriv(1), pw_grid)
1369 :
1370 1092 : IF (use_virial) CALL virial_laplace(rho_r(ispin), pw_pool, virial_xc, pw_to_deriv(1)%array)
1371 :
1372 1092 : CALL xc_pw_laplace(pw_to_deriv(1), pw_pool, xc_deriv_method_id)
1373 :
1374 1092 : CALL pw_axpy(pw_to_deriv(1), vxc_rho(ispin))
1375 :
1376 1092 : CALL pw_pool%give_back_pw(pw_to_deriv(1))
1377 : END IF
1378 :
1379 132926 : 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 132926 : rho_smooth_cutoff_range=density_smooth_cut_range)
1387 :
1388 132926 : v_drho_r = vxc_rho(ispin)
1389 132926 : CALL pw_pool%create_pw(vxc_rho(ispin))
1390 132926 : CALL xc_pw_smooth(v_drho_r, vxc_rho(ispin), xc_rho_smooth_id)
1391 243799 : CALL pw_pool%give_back_pw(v_drho_r)
1392 : END DO
1393 :
1394 110873 : CALL pw_pool%give_back_pw(vxc_g)
1395 110873 : 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 110873 : IF (has_derivs) THEN
1400 110555 : 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 110555 : rho_smooth_cutoff_range=density_smooth_cut_range)
1405 :
1406 110555 : exc = pw_integrate_function(v_drho_r)
1407 : !
1408 : ! return the xc functional value at the grid points
1409 : !
1410 110555 : IF (PRESENT(exc_r)) THEN
1411 96 : exc_r = v_drho_r
1412 : ELSE
1413 110459 : CALL v_drho_r%release()
1414 : END IF
1415 : ELSE
1416 318 : exc = 0.0_dp
1417 : END IF
1418 :
1419 110873 : CALL xc_rho_set_release(rho_set, pw_pool=pw_pool)
1420 :
1421 : ! tau part
1422 110873 : IF (has_tau) THEN
1423 8230 : ALLOCATE (vxc_tau(nspins))
1424 2532 : IF (lsd) THEN
1425 634 : CALL xc_dset_recover_pw(deriv_set, [deriv_tau_a], vxc_tau(1), pw_grid)
1426 634 : CALL xc_dset_recover_pw(deriv_set, [deriv_tau_b], vxc_tau(2), pw_grid)
1427 : ELSE
1428 1898 : CALL xc_dset_recover_pw(deriv_set, [deriv_tau], vxc_tau(1), pw_grid)
1429 : END IF
1430 5698 : DO ispin = 1, nspins
1431 5698 : CPASSERT(ASSOCIATED(vxc_tau(ispin)%array))
1432 : END DO
1433 : END IF
1434 110873 : CALL xc_dset_release(deriv_set)
1435 :
1436 110873 : CALL timestop(handle)
1437 :
1438 2328333 : 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 20968 : 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 10484 : 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 10484 : CALL timeset(routineN, handle)
1473 10484 : NULLIFY (deriv, e_0)
1474 10484 : 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 10484 : calc_potential=.FALSE.)
1482 10484 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::])
1483 :
1484 10484 : IF (ASSOCIATED(deriv)) THEN
1485 10484 : CALL xc_derivative_get(deriv, deriv_data=e_0)
1486 :
1487 : CALL section_vals_val_get(xc_section, "DENSITY_CUTOFF", &
1488 10484 : r_val=rho_cutoff)
1489 : CALL section_vals_val_get(xc_section, "DENSITY_SMOOTH_CUTOFF_RANGE", &
1490 10484 : 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 10484 : rho_smooth_cutoff_range=density_smooth_cut_range)
1495 :
1496 10484 : exc = accurate_sum(e_0)*rho_r(1)%pw_grid%dvol
1497 10484 : IF (rho_r(1)%pw_grid%para%mode == PW_MODE_DISTRIBUTED) THEN
1498 10346 : CALL rho_r(1)%pw_grid%para%group%sum(exc)
1499 : END IF
1500 :
1501 10484 : CALL xc_rho_set_release(rho_set, pw_pool=pw_pool)
1502 10484 : CALL xc_dset_release(deriv_set)
1503 : END IF
1504 10484 : CALL timestop(handle)
1505 199196 : 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 do_excitations ...
1521 : !> \param do_triplet ...
1522 : !> \param compute_virial ...
1523 : !> \param virial_xc virial terms will be collected here
1524 : ! **************************************************************************************************
1525 13850 : SUBROUTINE xc_calc_2nd_deriv(v_xc, v_xc_tau, deriv_set, rho_set, rho1_r, rho1_g, tau1_r, &
1526 : pw_pool, xc_section, gapw, vxg, &
1527 : do_excitations, do_triplet, &
1528 : compute_virial, virial_xc)
1529 :
1530 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_xc, v_xc_tau
1531 : TYPE(xc_derivative_set_type) :: deriv_set
1532 : TYPE(xc_rho_set_type) :: rho_set
1533 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho1_r, tau1_r
1534 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho1_g
1535 : TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool
1536 : TYPE(section_vals_type), INTENT(IN), POINTER :: xc_section
1537 : LOGICAL, INTENT(IN) :: gapw
1538 : REAL(KIND=dp), DIMENSION(:, :, :, :), OPTIONAL, &
1539 : POINTER :: vxg
1540 : LOGICAL, INTENT(IN), OPTIONAL :: do_excitations, &
1541 : do_triplet, compute_virial
1542 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT), &
1543 : OPTIONAL :: virial_xc
1544 :
1545 : CHARACTER(len=*), PARAMETER :: routineN = 'xc_calc_2nd_deriv'
1546 :
1547 : INTEGER :: handle, ispin, nspins
1548 : INTEGER, DIMENSION(2, 3) :: bo
1549 : LOGICAL :: lsd, my_compute_virial, &
1550 : my_do_excitations, &
1551 : my_do_triplet
1552 : REAL(KIND=dp) :: fac
1553 : TYPE(section_vals_type), POINTER :: xc_fun_section
1554 : TYPE(xc_rho_cflags_type) :: needs
1555 : TYPE(xc_rho_set_type) :: rho1_set
1556 :
1557 13850 : CALL timeset(routineN, handle)
1558 :
1559 13850 : my_compute_virial = .FALSE.
1560 13850 : IF (PRESENT(compute_virial)) my_compute_virial = compute_virial
1561 :
1562 13850 : my_do_excitations = .FALSE.
1563 13850 : IF (PRESENT(do_excitations)) my_do_excitations = do_excitations
1564 :
1565 13850 : my_do_triplet = .FALSE.
1566 13850 : IF (PRESENT(do_triplet)) my_do_triplet = do_triplet
1567 :
1568 13850 : nspins = SIZE(rho1_r)
1569 13850 : lsd = (nspins == 2)
1570 13850 : IF (nspins == 1 .AND. my_do_excitations .AND. my_do_triplet) THEN
1571 0 : nspins = 2
1572 0 : lsd = .TRUE.
1573 : END IF
1574 :
1575 13850 : NULLIFY (v_xc, v_xc_tau)
1576 57256 : ALLOCATE (v_xc(nspins))
1577 29556 : DO ispin = 1, nspins
1578 15706 : CALL pw_pool%create_pw(v_xc(ispin))
1579 29556 : CALL pw_zero(v_xc(ispin))
1580 : END DO
1581 :
1582 13850 : xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
1583 13850 : needs = xc_functionals_get_needs(xc_fun_section, lsd, .TRUE.)
1584 :
1585 13850 : IF (needs%tau .OR. needs%tau_spin) THEN
1586 536 : IF (.NOT. ASSOCIATED(tau1_r)) &
1587 0 : CPABORT("Tau-dependent functionals requires allocated kinetic energy density grid")
1588 1736 : ALLOCATE (v_xc_tau(nspins))
1589 1200 : DO ispin = 1, nspins
1590 664 : CALL pw_pool%create_pw(v_xc_tau(ispin))
1591 14514 : CALL pw_zero(v_xc_tau(ispin))
1592 : END DO
1593 : END IF
1594 :
1595 13850 : IF (section_get_lval(xc_section, "2ND_DERIV_ANALYTICAL")) THEN
1596 : !------!
1597 : ! rho1 !
1598 : !------!
1599 135500 : bo = rho1_r(1)%pw_grid%bounds_local
1600 : ! create the place where to store the argument for the functionals
1601 : CALL xc_rho_set_create(rho1_set, bo, &
1602 : rho_cutoff=section_get_rval(xc_section, "DENSITY_CUTOFF"), &
1603 : drho_cutoff=section_get_rval(xc_section, "GRADIENT_CUTOFF"), &
1604 13550 : tau_cutoff=section_get_rval(xc_section, "TAU_CUTOFF"))
1605 :
1606 : ! calculate the arguments needed by the functionals
1607 : CALL xc_rho_set_update(rho1_set, rho1_r, rho1_g, tau1_r, needs, &
1608 : section_get_ival(xc_section, "XC_GRID%XC_DERIV"), &
1609 : section_get_ival(xc_section, "XC_GRID%XC_SMOOTH_RHO"), &
1610 13550 : pw_pool)
1611 :
1612 13550 : fac = 0._dp
1613 13550 : IF (nspins == 1 .AND. my_do_excitations) THEN
1614 1074 : IF (my_do_triplet) fac = -1.0_dp
1615 : END IF
1616 :
1617 : CALL xc_calc_2nd_deriv_analytical(v_xc, v_xc_tau, deriv_set, rho_set, &
1618 : rho1_set, pw_pool, xc_section, gapw, vxg=vxg, &
1619 13550 : tddfpt_fac=fac, compute_virial=compute_virial, virial_xc=virial_xc)
1620 :
1621 13550 : CALL xc_rho_set_release(rho1_set)
1622 :
1623 : ELSE
1624 300 : IF (gapw) CPABORT("Numerical 2nd derivatives not implemented with GAPW")
1625 :
1626 : CALL xc_calc_2nd_deriv_numerical(v_xc, v_xc_tau, rho_set, rho1_r, rho1_g, tau1_r, &
1627 : pw_pool, xc_section, &
1628 : my_do_excitations .AND. my_do_triplet, &
1629 300 : compute_virial, virial_xc, deriv_set)
1630 : END IF
1631 :
1632 13850 : CALL timestop(handle)
1633 :
1634 304700 : END SUBROUTINE xc_calc_2nd_deriv
1635 :
1636 : ! **************************************************************************************************
1637 : !> \brief calculates 2nd derivative numerically
1638 : !> \param v_xc potential to be calculated (has to be allocated already)
1639 : !> \param v_tau tau-part of the potential to be calculated (has to be allocated already)
1640 : !> \param rho_set KS density from xc_prep_2nd_deriv
1641 : !> \param rho1_r first-order density in r-space
1642 : !> \param rho1_g first-order density in g-space
1643 : !> \param tau1_r first-order kinetic-energy density in r-space
1644 : !> \param pw_pool pw pool for new grids
1645 : !> \param xc_section XC section to calculate the derivatives from
1646 : !> \param do_triplet ...
1647 : !> \param calc_virial whether to calculate virial terms
1648 : !> \param virial_xc collects stress tensor components (no metaGGAs!)
1649 : !> \param deriv_set deriv set from xc_prep_2nd_deriv (only for virials)
1650 : ! **************************************************************************************************
1651 340 : SUBROUTINE xc_calc_2nd_deriv_numerical(v_xc, v_tau, rho_set, rho1_r, rho1_g, tau1_r, &
1652 : pw_pool, xc_section, &
1653 : do_triplet, calc_virial, virial_xc, deriv_set)
1654 :
1655 : TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(IN), POINTER :: v_xc, v_tau
1656 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
1657 : TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(IN), POINTER :: rho1_r, tau1_r
1658 : TYPE(pw_c1d_gs_type), DIMENSION(:), INTENT(IN), POINTER :: rho1_g
1659 : TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool
1660 : TYPE(section_vals_type), INTENT(IN), POINTER :: xc_section
1661 : LOGICAL, INTENT(IN) :: do_triplet
1662 : LOGICAL, INTENT(IN), OPTIONAL :: calc_virial
1663 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT), &
1664 : OPTIONAL :: virial_xc
1665 : TYPE(xc_derivative_set_type), OPTIONAL :: deriv_set
1666 :
1667 : CHARACTER(len=*), PARAMETER :: routineN = 'xc_calc_2nd_deriv_numerical'
1668 : REAL(KIND=dp), DIMENSION(-4:4, 4), PARAMETER :: &
1669 : 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, &
1670 : 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, &
1671 : 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, &
1672 : 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])
1673 :
1674 : INTEGER :: handle, idir, ispin, nspins, istep, nsteps
1675 : INTEGER, DIMENSION(2, 3) :: bo
1676 : LOGICAL :: gradient_f, lsd, my_calc_virial, tau_f, laplace_f, rho_f
1677 : REAL(KIND=dp) :: exc, gradient_cut, h, weight, step, rho_cutoff
1678 340 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: dr1dr, dra1dra, drb1drb
1679 : REAL(KIND=dp), DIMENSION(3, 3) :: virial_dummy
1680 340 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: norm_drho, norm_drho2, norm_drho2a, &
1681 340 : norm_drho2b, norm_drhoa, norm_drhob, &
1682 1020 : rho, rho1, rho1a, rho1b, rhoa, rhob, &
1683 680 : tau_a, tau_b, tau, tau1, tau1a, tau1b, laplace, laplace1, &
1684 340 : laplacea, laplaceb, laplace1a, laplace1b, &
1685 680 : laplace2, laplace2a, laplace2b, deriv_data
1686 8160 : TYPE(cp_3d_r_cp_type), DIMENSION(3) :: drho, drho1, drho1a, drho1b, drhoa, drhob
1687 : TYPE(pw_r3d_rs_type) :: v_drho, v_drhoa, v_drhob
1688 340 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: vxc_rho, vxc_tau
1689 340 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
1690 340 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, tau_r
1691 : TYPE(pw_r3d_rs_type) :: virial_pw, v_laplace, v_laplacea, v_laplaceb
1692 : TYPE(section_vals_type), POINTER :: xc_fun_section
1693 : TYPE(xc_derivative_set_type) :: deriv_set1
1694 : TYPE(xc_rho_cflags_type) :: needs
1695 : TYPE(xc_rho_set_type) :: rho1_set, rho2_set
1696 :
1697 340 : CALL timeset(routineN, handle)
1698 :
1699 340 : my_calc_virial = .FALSE.
1700 340 : IF (PRESENT(calc_virial) .AND. PRESENT(virial_xc)) my_calc_virial = calc_virial
1701 :
1702 340 : nspins = SIZE(v_xc)
1703 :
1704 340 : NULLIFY (tau, tau_r, tau_a, tau_b)
1705 :
1706 340 : h = section_get_rval(xc_section, "STEP_SIZE")
1707 340 : nsteps = section_get_ival(xc_section, "NSTEPS")
1708 340 : IF (nsteps < LBOUND(weights, 2) .OR. nspins > UBOUND(weights, 2)) THEN
1709 0 : CPABORT("The number of steps must be a value from 1 to 4.")
1710 : END IF
1711 :
1712 340 : IF (nspins == 2) THEN
1713 148 : NULLIFY (vxc_rho, rho_g, vxc_tau)
1714 444 : ALLOCATE (rho_r(2))
1715 444 : DO ispin = 1, nspins
1716 444 : CALL pw_pool%create_pw(rho_r(ispin))
1717 : END DO
1718 148 : IF (ASSOCIATED(tau1_r) .AND. ASSOCIATED(v_tau)) THEN
1719 162 : ALLOCATE (tau_r(2))
1720 162 : DO ispin = 1, nspins
1721 162 : CALL pw_pool%create_pw(tau_r(ispin))
1722 : END DO
1723 : END IF
1724 148 : CALL xc_rho_set_get(rho_set, can_return_null=.TRUE., rhoa=rhoa, rhob=rhob, tau_a=tau_a, tau_b=tau_b)
1725 1088 : DO istep = -nsteps, nsteps
1726 940 : IF (istep == 0) CYCLE
1727 792 : weight = weights(istep, nsteps)/h
1728 792 : step = REAL(istep, dp)*h
1729 : CALL calc_resp_potential_numer_ab(rho_r, rho_g, rho1_r, rhoa, rhob, vxc_rho, &
1730 792 : tau_r, tau1_r, tau_a, tau_b, vxc_tau, xc_section, pw_pool, step)
1731 2376 : DO ispin = 1, nspins
1732 1584 : CALL pw_axpy(vxc_rho(ispin), v_xc(ispin), weight)
1733 2376 : IF (ASSOCIATED(vxc_tau) .AND. ASSOCIATED(v_tau)) THEN
1734 456 : CALL pw_axpy(vxc_tau(ispin), v_tau(ispin), weight)
1735 : END IF
1736 : END DO
1737 2376 : DO ispin = 1, nspins
1738 2376 : CALL vxc_rho(ispin)%release()
1739 : END DO
1740 792 : DEALLOCATE (vxc_rho)
1741 940 : IF (ASSOCIATED(vxc_tau)) THEN
1742 684 : DO ispin = 1, nspins
1743 684 : CALL vxc_tau(ispin)%release()
1744 : END DO
1745 228 : DEALLOCATE (vxc_tau)
1746 : END IF
1747 : END DO
1748 192 : ELSE IF (nspins == 1 .AND. do_triplet) THEN
1749 20 : NULLIFY (vxc_rho, vxc_tau, rho_g)
1750 60 : ALLOCATE (rho_r(2))
1751 60 : DO ispin = 1, 2
1752 60 : CALL pw_pool%create_pw(rho_r(ispin))
1753 : END DO
1754 20 : IF (ASSOCIATED(tau1_r) .AND. ASSOCIATED(v_tau)) THEN
1755 0 : ALLOCATE (tau_r(2))
1756 0 : DO ispin = 1, nspins
1757 0 : CALL pw_pool%create_pw(tau_r(ispin))
1758 : END DO
1759 : END IF
1760 20 : CALL xc_rho_set_get(rho_set, can_return_null=.TRUE., rhoa=rhoa, rhob=rhob, tau_a=tau_a, tau_b=tau_b)
1761 160 : DO istep = -nsteps, nsteps
1762 140 : IF (istep == 0) CYCLE
1763 120 : weight = weights(istep, nsteps)/h
1764 120 : step = REAL(istep, dp)*h
1765 : ! K(alpha,alpha)
1766 120 : !$OMP PARALLEL DEFAULT(NONE) SHARED(rho_r,rhoa,rhob,step,rho1_r,tau_r,tau_a,tau_b,tau1_r)
1767 : !$OMP WORKSHARE
1768 : rho_r(1)%array(:, :, :) = rhoa(:, :, :) + step*rho1_r(1)%array(:, :, :)
1769 : !$OMP END WORKSHARE NOWAIT
1770 : !$OMP WORKSHARE
1771 : rho_r(2)%array(:, :, :) = rhob(:, :, :)
1772 : !$OMP END WORKSHARE NOWAIT
1773 : IF (ASSOCIATED(tau1_r)) THEN
1774 : !$OMP WORKSHARE
1775 : tau_r(1)%array(:, :, :) = tau_a(:, :, :) + step*tau1_r(1)%array(:, :, :)
1776 : !$OMP END WORKSHARE NOWAIT
1777 : !$OMP WORKSHARE
1778 : tau_r(2)%array(:, :, :) = tau_b(:, :, :)
1779 : !$OMP END WORKSHARE NOWAIT
1780 : END IF
1781 : !$OMP END PARALLEL
1782 : CALL xc_vxc_pw_create(vxc_rho, vxc_tau, exc, rho_r, rho_g, tau_r, xc_section, &
1783 120 : pw_pool, .FALSE., virial_dummy)
1784 120 : CALL pw_axpy(vxc_rho(1), v_xc(1), weight)
1785 120 : IF (ASSOCIATED(vxc_tau) .AND. ASSOCIATED(v_tau)) THEN
1786 0 : CALL pw_axpy(vxc_tau(1), v_tau(1), weight)
1787 : END IF
1788 360 : DO ispin = 1, 2
1789 360 : CALL vxc_rho(ispin)%release()
1790 : END DO
1791 120 : DEALLOCATE (vxc_rho)
1792 120 : IF (ASSOCIATED(vxc_tau)) THEN
1793 0 : DO ispin = 1, 2
1794 0 : CALL vxc_tau(ispin)%release()
1795 : END DO
1796 0 : DEALLOCATE (vxc_tau)
1797 : END IF
1798 120 : !$OMP PARALLEL DEFAULT(NONE) SHARED(rho_r,rhoa,rhob,step,rho1_r,tau_r,tau_a,tau_b,tau1_r)
1799 : !$OMP WORKSHARE
1800 : ! K(alpha,beta)
1801 : rho_r(1)%array(:, :, :) = rhoa(:, :, :)
1802 : !$OMP END WORKSHARE NOWAIT
1803 : !$OMP WORKSHARE
1804 : rho_r(2)%array(:, :, :) = rhob(:, :, :) + step*rho1_r(1)%array(:, :, :)
1805 : !$OMP END WORKSHARE NOWAIT
1806 : IF (ASSOCIATED(tau1_r)) THEN
1807 : !$OMP WORKSHARE
1808 : tau_r(1)%array(:, :, :) = tau_a(:, :, :)
1809 : !$OMP END WORKSHARE NOWAIT
1810 : !$OMP WORKSHARE
1811 : tau_r(2)%array(:, :, :) = tau_b(:, :, :) + step*tau1_r(1)%array(:, :, :)
1812 : !$OMP END WORKSHARE NOWAIT
1813 : END IF
1814 : !$OMP END PARALLEL
1815 : CALL xc_vxc_pw_create(vxc_rho, vxc_tau, exc, rho_r, rho_g, tau_r, xc_section, &
1816 120 : pw_pool, .FALSE., virial_dummy)
1817 120 : CALL pw_axpy(vxc_rho(1), v_xc(1), weight)
1818 120 : IF (ASSOCIATED(vxc_tau) .AND. ASSOCIATED(v_tau)) THEN
1819 0 : CALL pw_axpy(vxc_tau(1), v_tau(1), weight)
1820 : END IF
1821 360 : DO ispin = 1, 2
1822 360 : CALL vxc_rho(ispin)%release()
1823 : END DO
1824 120 : DEALLOCATE (vxc_rho)
1825 140 : IF (ASSOCIATED(vxc_tau)) THEN
1826 0 : DO ispin = 1, 2
1827 0 : CALL vxc_tau(ispin)%release()
1828 : END DO
1829 0 : DEALLOCATE (vxc_tau)
1830 : END IF
1831 : END DO
1832 : ELSE
1833 172 : NULLIFY (vxc_rho, rho_r, rho_g, vxc_tau, tau_r, tau)
1834 344 : ALLOCATE (rho_r(1))
1835 172 : CALL pw_pool%create_pw(rho_r(1))
1836 172 : IF (ASSOCIATED(tau1_r) .AND. ASSOCIATED(v_tau)) THEN
1837 92 : ALLOCATE (tau_r(1))
1838 46 : CALL pw_pool%create_pw(tau_r(1))
1839 : END IF
1840 172 : CALL xc_rho_set_get(rho_set, can_return_null=.TRUE., rho=rho, tau=tau)
1841 1280 : DO istep = -nsteps, nsteps
1842 1108 : IF (istep == 0) CYCLE
1843 936 : weight = weights(istep, nsteps)/h
1844 936 : step = REAL(istep, dp)*h
1845 936 : !$OMP PARALLEL DEFAULT(NONE) SHARED(rho_r,rho,step,rho1_r,tau1_r,tau,tau_r)
1846 : !$OMP WORKSHARE
1847 : rho_r(1)%array(:, :, :) = rho(:, :, :) + step*rho1_r(1)%array(:, :, :)
1848 : !$OMP END WORKSHARE NOWAIT
1849 : IF (ASSOCIATED(tau1_r) .AND. ASSOCIATED(tau) .AND. ASSOCIATED(tau1_r)) THEN
1850 : !$OMP WORKSHARE
1851 : tau_r(1)%array(:, :, :) = tau(:, :, :) + step*tau1_r(1)%array(:, :, :)
1852 : !$OMP END WORKSHARE NOWAIT
1853 : END IF
1854 : !$OMP END PARALLEL
1855 : CALL xc_vxc_pw_create(vxc_rho, vxc_tau, exc, rho_r, rho_g, tau_r, xc_section, &
1856 936 : pw_pool, .FALSE., virial_dummy)
1857 936 : CALL pw_axpy(vxc_rho(1), v_xc(1), weight)
1858 936 : IF (ASSOCIATED(vxc_tau) .AND. ASSOCIATED(v_tau)) THEN
1859 276 : CALL pw_axpy(vxc_tau(1), v_tau(1), weight)
1860 : END IF
1861 936 : CALL vxc_rho(1)%release()
1862 936 : DEALLOCATE (vxc_rho)
1863 1108 : IF (ASSOCIATED(vxc_tau)) THEN
1864 276 : CALL vxc_tau(1)%release()
1865 276 : DEALLOCATE (vxc_tau)
1866 : END IF
1867 : END DO
1868 : END IF
1869 :
1870 340 : IF (my_calc_virial) THEN
1871 36 : lsd = (nspins == 2)
1872 36 : IF (nspins == 1 .AND. do_triplet) THEN
1873 0 : lsd = .TRUE.
1874 : END IF
1875 :
1876 36 : CALL check_for_derivatives(deriv_set, (nspins == 2), rho_f, gradient_f, tau_f, laplace_f)
1877 :
1878 : ! Calculate the virial terms
1879 : ! Those arising from the first derivatives are treated like in xc_calc_2nd_deriv_analytical
1880 : ! Those arising from the second derivatives are calculated numerically
1881 : ! We assume that all metaGGA functionals require the gradient
1882 36 : IF (gradient_f) THEN
1883 360 : bo = rho_set%local_bounds
1884 :
1885 : ! Create the work grid for the virial terms
1886 36 : CALL allocate_pw(virial_pw, pw_pool, bo)
1887 :
1888 36 : gradient_cut = section_get_rval(xc_section, "GRADIENT_CUTOFF")
1889 :
1890 : ! create the container to store the argument of the functionals
1891 : CALL xc_rho_set_create(rho1_set, bo, &
1892 : rho_cutoff=section_get_rval(xc_section, "DENSITY_CUTOFF"), &
1893 : drho_cutoff=gradient_cut, &
1894 36 : tau_cutoff=section_get_rval(xc_section, "TAU_CUTOFF"))
1895 :
1896 36 : xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
1897 36 : needs = xc_functionals_get_needs(xc_fun_section, lsd, .TRUE.)
1898 :
1899 : ! calculate the arguments needed by the functionals
1900 : CALL xc_rho_set_update(rho1_set, rho1_r, rho1_g, tau1_r, needs, &
1901 : section_get_ival(xc_section, "XC_GRID%XC_DERIV"), &
1902 : section_get_ival(xc_section, "XC_GRID%XC_SMOOTH_RHO"), &
1903 36 : pw_pool)
1904 :
1905 36 : IF (lsd) THEN
1906 : CALL xc_rho_set_get(rho_set, drhoa=drhoa, drhob=drhob, norm_drho=norm_drho, &
1907 : norm_drhoa=norm_drhoa, norm_drhob=norm_drhob, tau_a=tau_a, tau_b=tau_b, &
1908 10 : laplace_rhoa=laplacea, laplace_rhob=laplaceb, can_return_null=.TRUE.)
1909 : CALL xc_rho_set_get(rho1_set, rhoa=rho1a, rhob=rho1b, drhoa=drho1a, drhob=drho1b, laplace_rhoa=laplace1a, &
1910 10 : laplace_rhob=laplace1b, can_return_null=.TRUE.)
1911 :
1912 10 : CALL calc_drho_from_ab(drho, drhoa, drhob)
1913 10 : CALL calc_drho_from_ab(drho1, drho1a, drho1b)
1914 : ELSE
1915 26 : CALL xc_rho_set_get(rho_set, drho=drho, norm_drho=norm_drho, tau=tau, laplace_rho=laplace, can_return_null=.TRUE.)
1916 26 : CALL xc_rho_set_get(rho1_set, rho=rho1, drho=drho1, laplace_rho=laplace1, can_return_null=.TRUE.)
1917 : END IF
1918 :
1919 36 : CALL prepare_dr1dr(dr1dr, drho, drho1)
1920 :
1921 36 : IF (lsd) THEN
1922 10 : CALL prepare_dr1dr(dra1dra, drhoa, drho1a)
1923 10 : CALL prepare_dr1dr(drb1drb, drhob, drho1b)
1924 :
1925 10 : CALL allocate_pw(v_drho, pw_pool, bo)
1926 10 : CALL allocate_pw(v_drhoa, pw_pool, bo)
1927 10 : CALL allocate_pw(v_drhob, pw_pool, bo)
1928 :
1929 10 : IF (ASSOCIATED(norm_drhoa)) CALL apply_drho(deriv_set, [deriv_norm_drhoa], virial_pw, drhoa, drho1a, virial_xc, &
1930 10 : norm_drhoa, gradient_cut, dra1dra, v_drhoa%array)
1931 10 : IF (ASSOCIATED(norm_drhob)) CALL apply_drho(deriv_set, [deriv_norm_drhob], virial_pw, drhob, drho1b, virial_xc, &
1932 10 : norm_drhob, gradient_cut, drb1drb, v_drhob%array)
1933 10 : IF (ASSOCIATED(norm_drho)) CALL apply_drho(deriv_set, [deriv_norm_drho], virial_pw, drho, drho1, virial_xc, &
1934 6 : norm_drho, gradient_cut, dr1dr, v_drho%array)
1935 10 : IF (laplace_f) THEN
1936 2 : CALL xc_derivative_get(xc_dset_get_derivative(deriv_set, [deriv_laplace_rhoa]), deriv_data=deriv_data)
1937 2 : CPASSERT(ASSOCIATED(deriv_data))
1938 15026 : virial_pw%array(:, :, :) = -rho1a(:, :, :)
1939 2 : CALL virial_laplace(virial_pw, pw_pool, virial_xc, deriv_data)
1940 :
1941 2 : CALL allocate_pw(v_laplacea, pw_pool, bo)
1942 :
1943 2 : CALL xc_derivative_get(xc_dset_get_derivative(deriv_set, [deriv_laplace_rhob]), deriv_data=deriv_data)
1944 2 : CPASSERT(ASSOCIATED(deriv_data))
1945 15026 : virial_pw%array(:, :, :) = -rho1b(:, :, :)
1946 2 : CALL virial_laplace(virial_pw, pw_pool, virial_xc, deriv_data)
1947 :
1948 2 : CALL allocate_pw(v_laplaceb, pw_pool, bo)
1949 : END IF
1950 :
1951 : ELSE
1952 :
1953 : ! Create the work grid for the potential of the gradient part
1954 26 : CALL allocate_pw(v_drho, pw_pool, bo)
1955 :
1956 : CALL apply_drho(deriv_set, [deriv_norm_drho], virial_pw, drho, drho1, virial_xc, &
1957 26 : norm_drho, gradient_cut, dr1dr, v_drho%array)
1958 26 : IF (laplace_f) THEN
1959 2 : CALL xc_derivative_get(xc_dset_get_derivative(deriv_set, [deriv_laplace_rho]), deriv_data=deriv_data)
1960 2 : CPASSERT(ASSOCIATED(deriv_data))
1961 28862 : virial_pw%array(:, :, :) = -rho1(:, :, :)
1962 2 : CALL virial_laplace(virial_pw, pw_pool, virial_xc, deriv_data)
1963 :
1964 2 : CALL allocate_pw(v_laplace, pw_pool, bo)
1965 : END IF
1966 :
1967 : END IF
1968 :
1969 36 : IF (lsd) THEN
1970 150260 : rho_r(1)%array = rhoa
1971 150260 : rho_r(2)%array = rhob
1972 : ELSE
1973 701552 : rho_r(1)%array = rho
1974 : END IF
1975 36 : IF (ASSOCIATED(tau1_r)) THEN
1976 8 : IF (lsd) THEN
1977 60104 : tau_r(1)%array = tau_a
1978 60104 : tau_r(2)%array = tau_b
1979 : ELSE
1980 115448 : tau_r(1)%array = tau
1981 : END IF
1982 : END IF
1983 :
1984 : ! Create deriv sets with same densities but different gradients
1985 36 : CALL xc_dset_create(deriv_set1, pw_pool)
1986 :
1987 36 : rho_cutoff = section_get_rval(xc_section, "DENSITY_CUTOFF")
1988 :
1989 : ! create the place where to store the argument for the functionals
1990 : CALL xc_rho_set_create(rho2_set, bo, &
1991 : rho_cutoff=rho_cutoff, &
1992 : drho_cutoff=section_get_rval(xc_section, "GRADIENT_CUTOFF"), &
1993 36 : tau_cutoff=section_get_rval(xc_section, "TAU_CUTOFF"))
1994 :
1995 : ! calculate the arguments needed by the functionals
1996 : CALL xc_rho_set_update(rho2_set, rho_r, rho_g, tau_r, needs, &
1997 : section_get_ival(xc_section, "XC_GRID%XC_DERIV"), &
1998 : section_get_ival(xc_section, "XC_GRID%XC_SMOOTH_RHO"), &
1999 36 : pw_pool)
2000 :
2001 36 : IF (lsd) THEN
2002 : CALL xc_rho_set_get(rho1_set, rhoa=rho1a, rhob=rho1b, tau_a=tau1a, tau_b=tau1b, &
2003 10 : laplace_rhoa=laplace1a, laplace_rhob=laplace1b, can_return_null=.TRUE.)
2004 : CALL xc_rho_set_get(rho2_set, norm_drhoa=norm_drho2a, norm_drhob=norm_drho2b, &
2005 10 : norm_drho=norm_drho2, laplace_rhoa=laplace2a, laplace_rhob=laplace2b, can_return_null=.TRUE.)
2006 :
2007 64 : DO istep = -nsteps, nsteps
2008 54 : IF (istep == 0) CYCLE
2009 44 : weight = weights(istep, nsteps)/h
2010 44 : step = REAL(istep, dp)*h
2011 44 : IF (ASSOCIATED(norm_drhoa)) THEN
2012 44 : CALL get_derivs_rho(norm_drho2a, norm_drhoa, step, xc_fun_section, lsd, rho2_set, deriv_set1)
2013 : CALL update_deriv_rho(deriv_set1, [deriv_rhoa], bo, &
2014 44 : norm_drhoa, gradient_cut, weight, rho1a, v_drhoa%array)
2015 : CALL update_deriv_rho(deriv_set1, [deriv_rhob], bo, &
2016 44 : norm_drhoa, gradient_cut, weight, rho1b, v_drhoa%array)
2017 : CALL update_deriv_rho(deriv_set1, [deriv_norm_drhoa], bo, &
2018 44 : norm_drhoa, gradient_cut, weight, dra1dra, v_drhoa%array)
2019 : CALL update_deriv_drho_ab(deriv_set1, [deriv_norm_drhob], bo, &
2020 44 : norm_drhoa, gradient_cut, weight, dra1dra, drb1drb, v_drhoa%array, v_drhob%array)
2021 : CALL update_deriv_drho_ab(deriv_set1, [deriv_norm_drho], bo, &
2022 44 : norm_drhoa, gradient_cut, weight, dra1dra, dr1dr, v_drhoa%array, v_drho%array)
2023 44 : IF (tau_f) THEN
2024 : CALL update_deriv_rho(deriv_set1, [deriv_tau_a], bo, &
2025 8 : norm_drhoa, gradient_cut, weight, tau1a, v_drhoa%array)
2026 : CALL update_deriv_rho(deriv_set1, [deriv_tau_b], bo, &
2027 8 : norm_drhoa, gradient_cut, weight, tau1b, v_drhoa%array)
2028 : END IF
2029 44 : IF (laplace_f) THEN
2030 : CALL update_deriv_rho(deriv_set1, [deriv_laplace_rhoa], bo, &
2031 4 : norm_drhoa, gradient_cut, weight, laplace1a, v_drhoa%array)
2032 : CALL update_deriv_rho(deriv_set1, [deriv_laplace_rhob], bo, &
2033 4 : norm_drhoa, gradient_cut, weight, laplace1b, v_drhoa%array)
2034 : END IF
2035 : END IF
2036 :
2037 44 : IF (ASSOCIATED(norm_drhob)) THEN
2038 44 : CALL get_derivs_rho(norm_drho2b, norm_drhob, step, xc_fun_section, lsd, rho2_set, deriv_set1)
2039 : CALL update_deriv_rho(deriv_set1, [deriv_rhoa], bo, &
2040 44 : norm_drhob, gradient_cut, weight, rho1a, v_drhob%array)
2041 : CALL update_deriv_rho(deriv_set1, [deriv_rhob], bo, &
2042 44 : norm_drhob, gradient_cut, weight, rho1b, v_drhob%array)
2043 : CALL update_deriv_rho(deriv_set1, [deriv_norm_drhob], bo, &
2044 44 : norm_drhob, gradient_cut, weight, drb1drb, v_drhob%array)
2045 : CALL update_deriv_drho_ab(deriv_set1, [deriv_norm_drhoa], bo, &
2046 44 : norm_drhob, gradient_cut, weight, drb1drb, dra1dra, v_drhob%array, v_drhoa%array)
2047 : CALL update_deriv_drho_ab(deriv_set1, [deriv_norm_drho], bo, &
2048 44 : norm_drhob, gradient_cut, weight, drb1drb, dr1dr, v_drhob%array, v_drho%array)
2049 44 : IF (tau_f) THEN
2050 : CALL update_deriv_rho(deriv_set1, [deriv_tau_a], bo, &
2051 8 : norm_drhob, gradient_cut, weight, tau1a, v_drhob%array)
2052 : CALL update_deriv_rho(deriv_set1, [deriv_tau_b], bo, &
2053 8 : norm_drhob, gradient_cut, weight, tau1b, v_drhob%array)
2054 : END IF
2055 44 : IF (laplace_f) THEN
2056 : CALL update_deriv_rho(deriv_set1, [deriv_laplace_rhoa], bo, &
2057 4 : norm_drhob, gradient_cut, weight, laplace1a, v_drhob%array)
2058 : CALL update_deriv_rho(deriv_set1, [deriv_laplace_rhob], bo, &
2059 4 : norm_drhob, gradient_cut, weight, laplace1b, v_drhob%array)
2060 : END IF
2061 : END IF
2062 :
2063 44 : IF (ASSOCIATED(norm_drho)) THEN
2064 20 : CALL get_derivs_rho(norm_drho2, norm_drho, step, xc_fun_section, lsd, rho2_set, deriv_set1)
2065 : CALL update_deriv_rho(deriv_set1, [deriv_rhoa], bo, &
2066 20 : norm_drho, gradient_cut, weight, rho1a, v_drho%array)
2067 : CALL update_deriv_rho(deriv_set1, [deriv_rhob], bo, &
2068 20 : norm_drho, gradient_cut, weight, rho1b, v_drho%array)
2069 : CALL update_deriv_rho(deriv_set1, [deriv_norm_drho], bo, &
2070 20 : norm_drho, gradient_cut, weight, dr1dr, v_drho%array)
2071 : CALL update_deriv_drho_ab(deriv_set1, [deriv_norm_drhoa], bo, &
2072 20 : norm_drho, gradient_cut, weight, dr1dr, dra1dra, v_drho%array, v_drhoa%array)
2073 : CALL update_deriv_drho_ab(deriv_set1, [deriv_norm_drhob], bo, &
2074 20 : norm_drho, gradient_cut, weight, dr1dr, drb1drb, v_drho%array, v_drhob%array)
2075 20 : IF (tau_f) THEN
2076 : CALL update_deriv_rho(deriv_set1, [deriv_tau_a], bo, &
2077 8 : norm_drho, gradient_cut, weight, tau1a, v_drho%array)
2078 : CALL update_deriv_rho(deriv_set1, [deriv_tau_b], bo, &
2079 8 : norm_drho, gradient_cut, weight, tau1b, v_drho%array)
2080 : END IF
2081 20 : IF (laplace_f) THEN
2082 : CALL update_deriv_rho(deriv_set1, [deriv_laplace_rhoa], bo, &
2083 4 : norm_drho, gradient_cut, weight, laplace1a, v_drho%array)
2084 : CALL update_deriv_rho(deriv_set1, [deriv_laplace_rhob], bo, &
2085 4 : norm_drho, gradient_cut, weight, laplace1b, v_drho%array)
2086 : END IF
2087 : END IF
2088 :
2089 54 : IF (laplace_f) THEN
2090 :
2091 4 : CALL get_derivs_rho(laplace2a, laplacea, step, xc_fun_section, lsd, rho2_set, deriv_set1)
2092 :
2093 : ! Obtain the numerical 2nd derivatives w.r.t. to drho and collect the potential
2094 : CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [deriv_rhoa], bo, &
2095 4 : weight, rho1a, v_laplacea%array)
2096 : CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [deriv_rhob], bo, &
2097 4 : weight, rho1b, v_laplacea%array)
2098 4 : IF (ASSOCIATED(norm_drho)) THEN
2099 : CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [deriv_norm_drho], bo, &
2100 4 : weight, dr1dr, v_laplacea%array)
2101 : END IF
2102 4 : IF (ASSOCIATED(norm_drhoa)) THEN
2103 : CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [deriv_norm_drhoa], bo, &
2104 4 : weight, dra1dra, v_laplacea%array)
2105 : END IF
2106 4 : IF (ASSOCIATED(norm_drhob)) THEN
2107 : CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [deriv_norm_drhob], bo, &
2108 4 : weight, drb1drb, v_laplacea%array)
2109 : END IF
2110 :
2111 4 : IF (ASSOCIATED(tau1a)) THEN
2112 : CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [deriv_tau_a], bo, &
2113 4 : weight, tau1a, v_laplacea%array)
2114 : END IF
2115 4 : IF (ASSOCIATED(tau1b)) THEN
2116 : CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [deriv_tau_b], bo, &
2117 4 : weight, tau1b, v_laplacea%array)
2118 : END IF
2119 :
2120 : CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [deriv_laplace_rhoa], bo, &
2121 4 : weight, laplace1a, v_laplacea%array)
2122 :
2123 : CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [deriv_laplace_rhob], bo, &
2124 4 : weight, laplace1b, v_laplacea%array)
2125 :
2126 : ! The same for the beta spin
2127 4 : CALL get_derivs_rho(laplace2b, laplaceb, step, xc_fun_section, lsd, rho2_set, deriv_set1)
2128 :
2129 : ! Obtain the numerical 2nd derivatives w.r.t. to drho and collect the potential
2130 : CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [deriv_rhoa], bo, &
2131 4 : weight, rho1a, v_laplaceb%array)
2132 : CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [deriv_rhob], bo, &
2133 4 : weight, rho1b, v_laplaceb%array)
2134 4 : IF (ASSOCIATED(norm_drho)) THEN
2135 : CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [deriv_norm_drho], bo, &
2136 4 : weight, dr1dr, v_laplaceb%array)
2137 : END IF
2138 4 : IF (ASSOCIATED(norm_drhoa)) THEN
2139 : CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [deriv_norm_drhoa], bo, &
2140 4 : weight, dra1dra, v_laplaceb%array)
2141 : END IF
2142 4 : IF (ASSOCIATED(norm_drhob)) THEN
2143 : CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [deriv_norm_drhob], bo, &
2144 4 : weight, drb1drb, v_laplaceb%array)
2145 : END IF
2146 :
2147 4 : IF (tau_f) THEN
2148 : CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [deriv_tau_a], bo, &
2149 4 : weight, tau1a, v_laplaceb%array)
2150 : CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [deriv_tau_b], bo, &
2151 4 : weight, tau1b, v_laplaceb%array)
2152 : END IF
2153 :
2154 : CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [deriv_laplace_rhoa], bo, &
2155 4 : weight, laplace1a, v_laplaceb%array)
2156 :
2157 : CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [deriv_laplace_rhob], bo, &
2158 4 : weight, laplace1b, v_laplaceb%array)
2159 : END IF
2160 : END DO
2161 :
2162 10 : CALL virial_drho_drho(virial_pw, drhoa, v_drhoa, virial_xc)
2163 10 : CALL virial_drho_drho(virial_pw, drhob, v_drhob, virial_xc)
2164 10 : CALL virial_drho_drho(virial_pw, drho, v_drho, virial_xc)
2165 :
2166 10 : CALL deallocate_pw(v_drho, pw_pool)
2167 10 : CALL deallocate_pw(v_drhoa, pw_pool)
2168 10 : CALL deallocate_pw(v_drhob, pw_pool)
2169 :
2170 10 : IF (laplace_f) THEN
2171 15026 : virial_pw%array(:, :, :) = -rhoa(:, :, :)
2172 2 : CALL virial_laplace(virial_pw, pw_pool, virial_xc, v_laplacea%array)
2173 2 : CALL deallocate_pw(v_laplacea, pw_pool)
2174 :
2175 15026 : virial_pw%array(:, :, :) = -rhob(:, :, :)
2176 2 : CALL virial_laplace(virial_pw, pw_pool, virial_xc, v_laplaceb%array)
2177 2 : CALL deallocate_pw(v_laplaceb, pw_pool)
2178 : END IF
2179 :
2180 10 : CALL deallocate_pw(virial_pw, pw_pool)
2181 :
2182 40 : DO idir = 1, 3
2183 30 : DEALLOCATE (drho(idir)%array)
2184 40 : DEALLOCATE (drho1(idir)%array)
2185 : END DO
2186 10 : DEALLOCATE (dra1dra, drb1drb)
2187 :
2188 : ELSE
2189 26 : CALL xc_rho_set_get(rho1_set, rho=rho1, tau=tau1, laplace_rho=laplace1, can_return_null=.TRUE.)
2190 26 : CALL xc_rho_set_get(rho2_set, norm_drho=norm_drho2, laplace_rho=laplace2, can_return_null=.TRUE.)
2191 :
2192 200 : DO istep = -nsteps, nsteps
2193 174 : IF (istep == 0) CYCLE
2194 148 : weight = weights(istep, nsteps)/h
2195 148 : step = REAL(istep, dp)*h
2196 148 : CALL get_derivs_rho(norm_drho2, norm_drho, step, xc_fun_section, lsd, rho2_set, deriv_set1)
2197 :
2198 : ! Obtain the numerical 2nd derivatives w.r.t. to drho and collect the potential
2199 : CALL update_deriv_rho(deriv_set1, [deriv_rho], bo, &
2200 148 : norm_drho, gradient_cut, weight, rho1, v_drho%array)
2201 : CALL update_deriv_rho(deriv_set1, [deriv_norm_drho], bo, &
2202 148 : norm_drho, gradient_cut, weight, dr1dr, v_drho%array)
2203 :
2204 148 : IF (tau_f) THEN
2205 : CALL update_deriv_rho(deriv_set1, [deriv_tau], bo, &
2206 24 : norm_drho, gradient_cut, weight, tau1, v_drho%array)
2207 : END IF
2208 174 : IF (laplace_f) THEN
2209 : CALL update_deriv_rho(deriv_set1, [deriv_laplace_rho], bo, &
2210 12 : norm_drho, gradient_cut, weight, laplace1, v_drho%array)
2211 :
2212 12 : CALL get_derivs_rho(laplace2, laplace, step, xc_fun_section, lsd, rho2_set, deriv_set1)
2213 :
2214 : ! Obtain the numerical 2nd derivatives w.r.t. to drho and collect the potential
2215 : CALL update_deriv(deriv_set1, laplace, rho_cutoff, [deriv_rho], bo, &
2216 12 : weight, rho1, v_laplace%array)
2217 : CALL update_deriv(deriv_set1, laplace, rho_cutoff, [deriv_norm_drho], bo, &
2218 12 : weight, dr1dr, v_laplace%array)
2219 :
2220 12 : IF (tau_f) THEN
2221 : CALL update_deriv(deriv_set1, laplace, rho_cutoff, [deriv_tau], bo, &
2222 12 : weight, tau1, v_laplace%array)
2223 : END IF
2224 :
2225 : CALL update_deriv(deriv_set1, laplace, rho_cutoff, [deriv_laplace_rho], bo, &
2226 12 : weight, laplace1, v_laplace%array)
2227 : END IF
2228 : END DO
2229 :
2230 : ! Calculate the virial contribution from the potential
2231 26 : CALL virial_drho_drho(virial_pw, drho, v_drho, virial_xc)
2232 :
2233 26 : CALL deallocate_pw(v_drho, pw_pool)
2234 :
2235 26 : IF (laplace_f) THEN
2236 28862 : virial_pw%array(:, :, :) = -rho(:, :, :)
2237 2 : CALL virial_laplace(virial_pw, pw_pool, virial_xc, v_laplace%array)
2238 2 : CALL deallocate_pw(v_laplace, pw_pool)
2239 : END IF
2240 :
2241 26 : CALL deallocate_pw(virial_pw, pw_pool)
2242 : END IF
2243 :
2244 : END IF
2245 :
2246 36 : CALL xc_dset_release(deriv_set1)
2247 :
2248 36 : DEALLOCATE (dr1dr)
2249 :
2250 36 : CALL xc_rho_set_release(rho1_set)
2251 36 : CALL xc_rho_set_release(rho2_set)
2252 : END IF
2253 :
2254 848 : DO ispin = 1, SIZE(rho_r)
2255 848 : CALL pw_pool%give_back_pw(rho_r(ispin))
2256 : END DO
2257 340 : DEALLOCATE (rho_r)
2258 :
2259 340 : IF (ASSOCIATED(tau_r)) THEN
2260 254 : DO ispin = 1, SIZE(tau_r)
2261 254 : CALL pw_pool%give_back_pw(tau_r(ispin))
2262 : END DO
2263 100 : DEALLOCATE (tau_r)
2264 : END IF
2265 :
2266 340 : CALL timestop(handle)
2267 :
2268 15300 : END SUBROUTINE xc_calc_2nd_deriv_numerical
2269 :
2270 : ! **************************************************************************************************
2271 : !> \brief ...
2272 : !> \param rho_r ...
2273 : !> \param rho_g ...
2274 : !> \param rho1_r ...
2275 : !> \param rhoa ...
2276 : !> \param rhob ...
2277 : !> \param vxc_rho ...
2278 : !> \param tau_r ...
2279 : !> \param tau1_r ...
2280 : !> \param tau_a ...
2281 : !> \param tau_b ...
2282 : !> \param vxc_tau ...
2283 : !> \param xc_section ...
2284 : !> \param pw_pool ...
2285 : !> \param step ...
2286 : ! **************************************************************************************************
2287 792 : SUBROUTINE calc_resp_potential_numer_ab(rho_r, rho_g, rho1_r, rhoa, rhob, vxc_rho, &
2288 : tau_r, tau1_r, tau_a, tau_b, vxc_tau, &
2289 : xc_section, pw_pool, step)
2290 :
2291 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER, INTENT(IN) :: vxc_rho, vxc_tau
2292 : TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(IN) :: rho1_r
2293 : TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(IN), POINTER :: tau1_r
2294 : TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool
2295 : TYPE(section_vals_type), INTENT(IN), POINTER :: xc_section
2296 : REAL(KIND=dp), INTENT(IN) :: step
2297 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER, INTENT(IN) :: rhoa, rhob, tau_a, tau_b
2298 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER, INTENT(IN) :: rho_r
2299 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
2300 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: tau_r
2301 :
2302 : CHARACTER(len=*), PARAMETER :: routineN = 'calc_resp_potential_numer_ab'
2303 :
2304 : INTEGER :: handle
2305 : REAL(KIND=dp) :: exc
2306 : REAL(KIND=dp), DIMENSION(3, 3) :: virial_dummy
2307 :
2308 792 : CALL timeset(routineN, handle)
2309 :
2310 792 : !$OMP PARALLEL DEFAULT(NONE) SHARED(rho_r,rhoa,rhob,step,rho1_r,tau_r,tau_a,tau_b,tau1_r)
2311 : !$OMP WORKSHARE
2312 : rho_r(1)%array(:, :, :) = rhoa(:, :, :) + step*rho1_r(1)%array(:, :, :)
2313 : !$OMP END WORKSHARE NOWAIT
2314 : !$OMP WORKSHARE
2315 : rho_r(2)%array(:, :, :) = rhob(:, :, :) + step*rho1_r(2)%array(:, :, :)
2316 : !$OMP END WORKSHARE NOWAIT
2317 : IF (ASSOCIATED(tau1_r) .AND. ASSOCIATED(tau_r) .AND. ASSOCIATED(tau_a) .AND. ASSOCIATED(tau_b)) THEN
2318 : !$OMP WORKSHARE
2319 : tau_r(1)%array(:, :, :) = tau_a(:, :, :) + step*tau1_r(1)%array(:, :, :)
2320 : !$OMP END WORKSHARE NOWAIT
2321 : !$OMP WORKSHARE
2322 : tau_r(2)%array(:, :, :) = tau_b(:, :, :) + step*tau1_r(2)%array(:, :, :)
2323 : !$OMP END WORKSHARE NOWAIT
2324 : END IF
2325 : !$OMP END PARALLEL
2326 : CALL xc_vxc_pw_create(vxc_rho, vxc_tau, exc, rho_r, rho_g, tau_r, xc_section, &
2327 792 : pw_pool, .FALSE., virial_dummy)
2328 :
2329 792 : CALL timestop(handle)
2330 :
2331 792 : END SUBROUTINE calc_resp_potential_numer_ab
2332 :
2333 : ! **************************************************************************************************
2334 : !> \brief calculates stress tensor and potential contributions from the first derivative
2335 : !> \param deriv_set ...
2336 : !> \param description ...
2337 : !> \param virial_pw ...
2338 : !> \param drho ...
2339 : !> \param drho1 ...
2340 : !> \param virial_xc ...
2341 : !> \param norm_drho ...
2342 : !> \param gradient_cut ...
2343 : !> \param dr1dr ...
2344 : !> \param v_drho ...
2345 : ! **************************************************************************************************
2346 52 : SUBROUTINE apply_drho(deriv_set, description, virial_pw, drho, drho1, virial_xc, norm_drho, gradient_cut, dr1dr, v_drho)
2347 :
2348 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
2349 : INTEGER, DIMENSION(:), INTENT(in) :: description
2350 : TYPE(pw_r3d_rs_type), INTENT(IN) :: virial_pw
2351 : TYPE(cp_3d_r_cp_type), DIMENSION(3), INTENT(IN) :: drho, drho1
2352 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT) :: virial_xc
2353 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: norm_drho
2354 : REAL(KIND=dp), INTENT(IN) :: gradient_cut
2355 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: dr1dr
2356 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: v_drho
2357 :
2358 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_drho'
2359 :
2360 : INTEGER :: handle
2361 52 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: deriv_data
2362 : TYPE(xc_derivative_type), POINTER :: deriv_att
2363 :
2364 52 : CALL timeset(routineN, handle)
2365 :
2366 52 : deriv_att => xc_dset_get_derivative(deriv_set, description)
2367 52 : IF (ASSOCIATED(deriv_att)) THEN
2368 52 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
2369 52 : CALL virial_drho_drho1(virial_pw, drho, drho1, deriv_data, virial_xc)
2370 :
2371 52 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(dr1dr,gradient_cut,norm_drho,v_drho,deriv_data)
2372 : v_drho(:, :, :) = v_drho(:, :, :) + &
2373 : deriv_data(:, :, :)*dr1dr(:, :, :)/MAX(gradient_cut, norm_drho(:, :, :))**2
2374 : !$OMP END PARALLEL WORKSHARE
2375 : END IF
2376 :
2377 52 : CALL timestop(handle)
2378 :
2379 52 : END SUBROUTINE apply_drho
2380 :
2381 : ! **************************************************************************************************
2382 : !> \brief adds potential contributions from derivatives of rho or diagonal terms of norm_drho
2383 : !> \param deriv_set1 ...
2384 : !> \param description ...
2385 : !> \param bo ...
2386 : !> \param norm_drho norm_drho of which derivative is calculated
2387 : !> \param gradient_cut ...
2388 : !> \param h ...
2389 : !> \param rho1 function to contract the derivative with (rho1 for rho, dr1dr for norm_drho)
2390 : !> \param v_drho ...
2391 : ! **************************************************************************************************
2392 728 : SUBROUTINE update_deriv_rho(deriv_set1, description, bo, norm_drho, gradient_cut, weight, rho1, v_drho)
2393 :
2394 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set1
2395 : INTEGER, DIMENSION(:), INTENT(in) :: description
2396 : INTEGER, DIMENSION(2, 3), INTENT(IN) :: bo
2397 : REAL(KIND=dp), DIMENSION(bo(1, 1):bo(2, 1), bo(1, &
2398 : 2):bo(2, 2), bo(1, 3):bo(2, 3)), INTENT(IN) :: norm_drho
2399 : REAL(KIND=dp), INTENT(IN) :: gradient_cut, weight
2400 : REAL(KIND=dp), DIMENSION(bo(1, 1):bo(2, 1), bo(1, &
2401 : 2):bo(2, 2), bo(1, 3):bo(2, 3)), INTENT(IN) :: rho1
2402 : REAL(KIND=dp), DIMENSION(bo(1, 1):bo(2, 1), bo(1, &
2403 : 2):bo(2, 2), bo(1, 3):bo(2, 3)), INTENT(INOUT) :: v_drho
2404 :
2405 : CHARACTER(len=*), PARAMETER :: routineN = 'update_deriv_rho'
2406 :
2407 : INTEGER :: handle, i, j, k
2408 : REAL(KIND=dp) :: de
2409 728 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: deriv_data1
2410 : TYPE(xc_derivative_type), POINTER :: deriv_att1
2411 :
2412 728 : CALL timeset(routineN, handle)
2413 :
2414 : ! Obtain the numerical 2nd derivatives w.r.t. to drho and collect the potential
2415 728 : deriv_att1 => xc_dset_get_derivative(deriv_set1, description)
2416 728 : IF (ASSOCIATED(deriv_att1)) THEN
2417 728 : CALL xc_derivative_get(deriv_att1, deriv_data=deriv_data1)
2418 : !$OMP PARALLEL DO DEFAULT(NONE) &
2419 : !$OMP SHARED(bo,deriv_data1,weight,norm_drho,v_drho,rho1,gradient_cut) &
2420 : !$OMP PRIVATE(i,j,k,de) &
2421 728 : !$OMP COLLAPSE(3)
2422 : DO k = bo(1, 3), bo(2, 3)
2423 : DO j = bo(1, 2), bo(2, 2)
2424 : DO i = bo(1, 1), bo(2, 1)
2425 : de = weight*deriv_data1(i, j, k)/MAX(gradient_cut, norm_drho(i, j, k))**2
2426 : v_drho(i, j, k) = v_drho(i, j, k) - de*rho1(i, j, k)
2427 : END DO
2428 : END DO
2429 : END DO
2430 : !$OMP END PARALLEL DO
2431 : END IF
2432 :
2433 728 : CALL timestop(handle)
2434 :
2435 728 : END SUBROUTINE update_deriv_rho
2436 :
2437 : ! **************************************************************************************************
2438 : !> \brief adds potential contributions from derivatives of a component with positive and negative values
2439 : !> \param deriv_set1 ...
2440 : !> \param description ...
2441 : !> \param bo ...
2442 : !> \param h ...
2443 : !> \param rho1 function to contract the derivative with (rho1 for rho, dr1dr for norm_drho)
2444 : !> \param v ...
2445 : ! **************************************************************************************************
2446 120 : SUBROUTINE update_deriv(deriv_set1, rho, rho_cutoff, description, bo, weight, rho1, v)
2447 :
2448 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set1
2449 : INTEGER, DIMENSION(:), INTENT(in) :: description
2450 : INTEGER, DIMENSION(2, 3), INTENT(IN) :: bo
2451 : REAL(KIND=dp), INTENT(IN) :: weight, rho_cutoff
2452 : REAL(KIND=dp), DIMENSION(bo(1, 1):bo(2, 1), bo(1, &
2453 : 2):bo(2, 2), bo(1, 3):bo(2, 3)), INTENT(IN) :: rho, rho1
2454 : REAL(KIND=dp), DIMENSION(bo(1, 1):bo(2, 1), bo(1, &
2455 : 2):bo(2, 2), bo(1, 3):bo(2, 3)), INTENT(INOUT) :: v
2456 :
2457 : CHARACTER(len=*), PARAMETER :: routineN = 'update_deriv'
2458 :
2459 : INTEGER :: handle, i, j, k
2460 : REAL(KIND=dp) :: de
2461 120 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: deriv_data1
2462 : TYPE(xc_derivative_type), POINTER :: deriv_att1
2463 :
2464 120 : CALL timeset(routineN, handle)
2465 :
2466 : ! Obtain the numerical 2nd derivatives w.r.t. to drho and collect the potential
2467 120 : deriv_att1 => xc_dset_get_derivative(deriv_set1, description)
2468 120 : IF (ASSOCIATED(deriv_att1)) THEN
2469 120 : CALL xc_derivative_get(deriv_att1, deriv_data=deriv_data1)
2470 : !$OMP PARALLEL DO DEFAULT(NONE) &
2471 : !$OMP SHARED(bo,deriv_data1,weight,v,rho1,rho, rho_cutoff) &
2472 : !$OMP PRIVATE(i,j,k,de) &
2473 120 : !$OMP COLLAPSE(3)
2474 : DO k = bo(1, 3), bo(2, 3)
2475 : DO j = bo(1, 2), bo(2, 2)
2476 : DO i = bo(1, 1), bo(2, 1)
2477 : ! We have to consider that the given density (mostly the Laplacian) may have positive and negative values
2478 : de = weight*deriv_data1(i, j, k)/SIGN(MAX(ABS(rho(i, j, k)), rho_cutoff), rho(i, j, k))
2479 : v(i, j, k) = v(i, j, k) + de*rho1(i, j, k)
2480 : END DO
2481 : END DO
2482 : END DO
2483 : !$OMP END PARALLEL DO
2484 : END IF
2485 :
2486 120 : CALL timestop(handle)
2487 :
2488 120 : END SUBROUTINE update_deriv
2489 :
2490 : ! **************************************************************************************************
2491 : !> \brief adds mixed derivatives of norm_drho
2492 : !> \param deriv_set1 ...
2493 : !> \param description ...
2494 : !> \param bo ...
2495 : !> \param norm_drhoa norm_drho of which derivatives is calculated
2496 : !> \param gradient_cut ...
2497 : !> \param h ...
2498 : !> \param dra1dra dr1dr corresponding to norm_drho
2499 : !> \param drb1drb ...
2500 : !> \param v_drhoa potential corresponding to norm_drho
2501 : !> \param v_drhob ...
2502 : ! **************************************************************************************************
2503 216 : SUBROUTINE update_deriv_drho_ab(deriv_set1, description, bo, &
2504 216 : norm_drhoa, gradient_cut, weight, dra1dra, drb1drb, v_drhoa, v_drhob)
2505 :
2506 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set1
2507 : INTEGER, DIMENSION(:), INTENT(in) :: description
2508 : INTEGER, DIMENSION(2, 3), INTENT(IN) :: bo
2509 : REAL(KIND=dp), DIMENSION(bo(1, 1):bo(2, 1), bo(1, &
2510 : 2):bo(2, 2), bo(1, 3):bo(2, 3)), INTENT(IN) :: norm_drhoa
2511 : REAL(KIND=dp), INTENT(IN) :: gradient_cut, weight
2512 : REAL(KIND=dp), DIMENSION(bo(1, 1):bo(2, 1), bo(1, &
2513 : 2):bo(2, 2), bo(1, 3):bo(2, 3)), INTENT(IN) :: dra1dra, drb1drb
2514 : REAL(KIND=dp), DIMENSION(bo(1, 1):bo(2, 1), bo(1, &
2515 : 2):bo(2, 2), bo(1, 3):bo(2, 3)), INTENT(INOUT) :: v_drhoa, v_drhob
2516 :
2517 : CHARACTER(len=*), PARAMETER :: routineN = 'update_deriv_drho_ab'
2518 :
2519 : INTEGER :: handle, i, j, k
2520 : REAL(KIND=dp) :: de
2521 216 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: deriv_data1
2522 : TYPE(xc_derivative_type), POINTER :: deriv_att1
2523 :
2524 216 : CALL timeset(routineN, handle)
2525 :
2526 216 : deriv_att1 => xc_dset_get_derivative(deriv_set1, description)
2527 216 : IF (ASSOCIATED(deriv_att1)) THEN
2528 168 : CALL xc_derivative_get(deriv_att1, deriv_data=deriv_data1)
2529 : !$OMP PARALLEL DO DEFAULT(NONE) &
2530 : !$OMP PRIVATE(k,j,i,de) &
2531 : !$OMP SHARED(bo,drb1drb,dra1dra,deriv_data1,weight,gradient_cut,norm_drhoa,v_drhoa,v_drhob) &
2532 168 : !$OMP COLLAPSE(3)
2533 : DO k = bo(1, 3), bo(2, 3)
2534 : DO j = bo(1, 2), bo(2, 2)
2535 : DO i = bo(1, 1), bo(2, 1)
2536 : ! We introduce a factor of two because we will average between both numerical derivatives
2537 : de = 0.5_dp*weight*deriv_data1(i, j, k)/MAX(gradient_cut, norm_drhoa(i, j, k))**2
2538 : v_drhoa(i, j, k) = v_drhoa(i, j, k) - de*drb1drb(i, j, k)
2539 : v_drhob(i, j, k) = v_drhob(i, j, k) - de*dra1dra(i, j, k)
2540 : END DO
2541 : END DO
2542 : END DO
2543 : !$OMP END PARALLEL DO
2544 : END IF
2545 :
2546 216 : CALL timestop(handle)
2547 :
2548 216 : END SUBROUTINE update_deriv_drho_ab
2549 :
2550 : ! **************************************************************************************************
2551 : !> \brief calculate derivative sets for helper points
2552 : !> \param norm_drho2 norm_drho of new points
2553 : !> \param norm_drho norm_drho of KS density
2554 : !> \param h ...
2555 : !> \param xc_fun_section ...
2556 : !> \param lsd ...
2557 : !> \param rho2_set rho_set for new points
2558 : !> \param deriv_set1 will contain derivatives of the perturbed density
2559 : ! **************************************************************************************************
2560 276 : SUBROUTINE get_derivs_rho(norm_drho2, norm_drho, step, xc_fun_section, lsd, rho2_set, deriv_set1)
2561 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: norm_drho2
2562 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: norm_drho
2563 : REAL(KIND=dp), INTENT(IN) :: step
2564 : TYPE(section_vals_type), INTENT(IN), POINTER :: xc_fun_section
2565 : LOGICAL, INTENT(IN) :: lsd
2566 : TYPE(xc_rho_set_type), INTENT(INOUT) :: rho2_set
2567 : TYPE(xc_derivative_set_type) :: deriv_set1
2568 :
2569 : CHARACTER(len=*), PARAMETER :: routineN = 'get_derivs_rho'
2570 :
2571 : INTEGER :: handle
2572 :
2573 276 : CALL timeset(routineN, handle)
2574 :
2575 : ! Copy the densities, do one step into the direction of drho
2576 276 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(norm_drho,norm_drho2,step)
2577 : norm_drho2 = norm_drho*(1.0_dp + step)
2578 : !$OMP END PARALLEL WORKSHARE
2579 :
2580 276 : CALL xc_dset_zero_all(deriv_set1)
2581 :
2582 : ! Calculate the derivatives of the functional
2583 : CALL xc_functionals_eval(xc_fun_section, &
2584 : lsd=lsd, &
2585 : rho_set=rho2_set, &
2586 : deriv_set=deriv_set1, &
2587 276 : deriv_order=1)
2588 :
2589 : ! Return to the original values
2590 276 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(norm_drho,norm_drho2)
2591 : norm_drho2 = norm_drho
2592 : !$OMP END PARALLEL WORKSHARE
2593 :
2594 276 : CALL divide_by_norm_drho(deriv_set1, rho2_set, lsd)
2595 :
2596 276 : CALL timestop(handle)
2597 :
2598 276 : END SUBROUTINE get_derivs_rho
2599 :
2600 : ! **************************************************************************************************
2601 : !> \brief Calculates the second derivative of E_xc at rho in the direction
2602 : !> rho1 (if you see the second derivative as bilinear form)
2603 : !> partial_rho|_(rho=rho) partial_rho|_(rho=rho) E_xc drho(rho1)drho
2604 : !> The other direction is still undetermined, thus it returns
2605 : !> a potential (partial integration is performed to reduce it to
2606 : !> function of rho, removing the dependence from its partial derivs)
2607 : !> Has to be called after the setup by xc_prep_2nd_deriv.
2608 : !> \param v_xc exchange-correlation potential
2609 : !> \param v_xc_tau ...
2610 : !> \param deriv_set derivatives of the exchange-correlation potential
2611 : !> \param rho_set object containing the density at which the derivatives were calculated
2612 : !> \param rho1_set object containing the density with which to fold
2613 : !> \param pw_pool the pool for the grids
2614 : !> \param xc_section XC parameters
2615 : !> \param gapw Gaussian and augmented plane waves calculation
2616 : !> \param vxg ...
2617 : !> \param tddfpt_fac factor that multiplies the crossterms (tddfpt triplets
2618 : !> on a closed shell system it should be -1, defaults to 1)
2619 : !> \param compute_virial ...
2620 : !> \param virial_xc ...
2621 : !> \note
2622 : !> The old version of this routine was smarter: it handled split_desc(1)
2623 : !> and split_desc(2) separately, thus the code automatically handled all
2624 : !> possible cross terms (you only had to check if it was diagonal to avoid
2625 : !> double counting). I think that is the way to go if you want to add more
2626 : !> terms (tau,rho in LSD,...). The problem with the old code was that it
2627 : !> because of the old functional structure it sometime guessed wrongly
2628 : !> which derivative was where. There were probably still bugs with gradient
2629 : !> corrected functionals (never tested), and it didn't contain first
2630 : !> derivatives with respect to drho (that contribute also to the second
2631 : !> derivative wrt. rho).
2632 : !> The code was a little complex because it really tried to handle any
2633 : !> functional derivative in the most efficient way with the given contents of
2634 : !> rho_set.
2635 : !> Anyway I strongly encourage whoever wants to modify this code to give a
2636 : !> look to the old version. [fawzi]
2637 : ! **************************************************************************************************
2638 30342 : SUBROUTINE xc_calc_2nd_deriv_analytical(v_xc, v_xc_tau, deriv_set, rho_set, rho1_set, &
2639 : pw_pool, xc_section, gapw, vxg, tddfpt_fac, &
2640 : compute_virial, virial_xc)
2641 :
2642 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_xc, v_xc_tau
2643 : TYPE(xc_derivative_set_type) :: deriv_set
2644 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set, rho1_set
2645 : TYPE(pw_pool_type), POINTER :: pw_pool
2646 : TYPE(section_vals_type), POINTER :: xc_section
2647 : LOGICAL, INTENT(IN), OPTIONAL :: gapw
2648 : REAL(kind=dp), DIMENSION(:, :, :, :), OPTIONAL, &
2649 : POINTER :: vxg
2650 : REAL(kind=dp), INTENT(in), OPTIONAL :: tddfpt_fac
2651 : LOGICAL, INTENT(IN), OPTIONAL :: compute_virial
2652 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT), &
2653 : OPTIONAL :: virial_xc
2654 :
2655 : CHARACTER(len=*), PARAMETER :: routineN = 'xc_calc_2nd_deriv_analytical'
2656 :
2657 : INTEGER :: handle, i, ia, idir, ir, ispin, j, jdir, &
2658 : k, nspins, xc_deriv_method_id
2659 : INTEGER, DIMENSION(2, 3) :: bo
2660 : LOGICAL :: gradient_f, lsd, my_compute_virial, &
2661 : my_gapw, tau_f, laplace_f, rho_f
2662 : REAL(KIND=dp) :: fac, gradient_cut, tmp, factor2
2663 30342 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: dr1dr, dra1dra, drb1drb
2664 60684 : REAL(kind=dp), DIMENSION(:, :, :), POINTER :: deriv_data, e_drhoa, e_drhob, &
2665 30342 : e_drho, norm_drho, norm_drhoa, &
2666 30342 : norm_drhob, rho1, rho1a, rho1b, &
2667 30342 : tau1, tau1a, tau1b, laplace1, laplace1a, laplace1b, &
2668 30342 : rho, rhoa, rhob
2669 576498 : TYPE(cp_3d_r_cp_type), DIMENSION(3) :: drho, drho1, drho1a, drho1b, drhoa, drhob
2670 30342 : TYPE(pw_r3d_rs_type), DIMENSION(:), ALLOCATABLE :: v_drhoa, v_drhob, v_drho, v_laplace
2671 30342 : TYPE(pw_r3d_rs_type), DIMENSION(:, :), ALLOCATABLE :: v_drho_r
2672 : TYPE(pw_r3d_rs_type) :: virial_pw
2673 : TYPE(pw_c1d_gs_type) :: tmp_g, vxc_g
2674 : TYPE(xc_derivative_type), POINTER :: deriv_att
2675 :
2676 30342 : CALL timeset(routineN, handle)
2677 :
2678 30342 : NULLIFY (e_drhoa, e_drhob, e_drho)
2679 :
2680 30342 : my_gapw = .FALSE.
2681 30342 : IF (PRESENT(gapw)) my_gapw = gapw
2682 :
2683 30342 : my_compute_virial = .FALSE.
2684 30342 : IF (PRESENT(compute_virial)) my_compute_virial = compute_virial
2685 :
2686 30342 : CPASSERT(ASSOCIATED(v_xc))
2687 30342 : CPASSERT(ASSOCIATED(xc_section))
2688 30342 : IF (my_gapw) THEN
2689 9602 : CPASSERT(PRESENT(vxg))
2690 : END IF
2691 30342 : IF (my_compute_virial) THEN
2692 358 : CPASSERT(PRESENT(virial_xc))
2693 : END IF
2694 :
2695 : CALL section_vals_val_get(xc_section, "XC_GRID%XC_DERIV", &
2696 30342 : i_val=xc_deriv_method_id)
2697 30342 : CALL xc_rho_set_get(rho_set, drho_cutoff=gradient_cut)
2698 30342 : nspins = SIZE(v_xc)
2699 30342 : lsd = ASSOCIATED(rho_set%rhoa)
2700 30342 : fac = 0.0_dp
2701 30342 : factor2 = 1.0_dp
2702 30342 : IF (PRESENT(tddfpt_fac)) fac = tddfpt_fac
2703 30342 : IF (PRESENT(tddfpt_fac)) factor2 = tddfpt_fac
2704 :
2705 303420 : bo = rho_set%local_bounds
2706 :
2707 30342 : CALL check_for_derivatives(deriv_set, lsd, rho_f, gradient_f, tau_f, laplace_f)
2708 :
2709 30342 : IF (tau_f) THEN
2710 656 : CPASSERT(ASSOCIATED(v_xc_tau))
2711 : END IF
2712 :
2713 30342 : IF (gradient_f) THEN
2714 195870 : ALLOCATE (v_drho_r(3, nspins), v_drho(nspins))
2715 39174 : DO ispin = 1, nspins
2716 80928 : DO idir = 1, 3
2717 80928 : CALL allocate_pw(v_drho_r(idir, ispin), pw_pool, bo)
2718 : END DO
2719 39174 : CALL allocate_pw(v_drho(ispin), pw_pool, bo)
2720 : END DO
2721 :
2722 18942 : IF (xc_requires_tmp_g(xc_deriv_method_id) .AND. .NOT. my_gapw) THEN
2723 11482 : IF (ASSOCIATED(pw_pool)) THEN
2724 11482 : CALL pw_pool%create_pw(tmp_g)
2725 11482 : CALL pw_pool%create_pw(vxc_g)
2726 : ELSE
2727 : ! remember to refix for gapw
2728 0 : CPABORT("XC_DERIV method is not implemented in GAPW")
2729 : END IF
2730 : END IF
2731 : END IF
2732 :
2733 63878 : DO ispin = 1, nspins
2734 945940500 : v_xc(ispin)%array = 0.0_dp
2735 : END DO
2736 :
2737 30342 : IF (tau_f) THEN
2738 1386 : DO ispin = 1, nspins
2739 14897758 : v_xc_tau(ispin)%array = 0.0_dp
2740 : END DO
2741 : END IF
2742 :
2743 30342 : IF (laplace_f .AND. my_gapw) &
2744 0 : CPABORT("Laplace-dependent functional not implemented with GAPW!")
2745 :
2746 30342 : IF (my_compute_virial .AND. (gradient_f .OR. laplace_f)) CALL allocate_pw(virial_pw, pw_pool, bo)
2747 :
2748 30342 : IF (lsd) THEN
2749 :
2750 : !-------------------!
2751 : ! UNrestricted case !
2752 : !-------------------!
2753 :
2754 4232 : CALL xc_rho_set_get(rho1_set, rhoa=rho1a, rhob=rho1b)
2755 :
2756 4232 : IF (gradient_f) THEN
2757 : CALL xc_rho_set_get(rho_set, drhoa=drhoa, drhob=drhob, &
2758 2000 : norm_drho=norm_drho, norm_drhoa=norm_drhoa, norm_drhob=norm_drhob)
2759 2000 : CALL xc_rho_set_get(rho1_set, drhoa=drho1a, drhob=drho1b)
2760 :
2761 2000 : CALL calc_drho_from_ab(drho, drhoa, drhob)
2762 2000 : CALL calc_drho_from_ab(drho1, drho1a, drho1b)
2763 :
2764 2000 : CALL prepare_dr1dr(dra1dra, drhoa, drho1a)
2765 2000 : IF (nspins /= 1) THEN
2766 1290 : CALL prepare_dr1dr(drb1drb, drhob, drho1b)
2767 1290 : CALL prepare_dr1dr(dr1dr, drho, drho1)
2768 : ELSE
2769 710 : CALL prepare_dr1dr(drb1drb, drhob, drho1b)
2770 710 : CALL prepare_dr1dr_ab(dr1dr, drhoa, drhob, drho1a, drho1b, fac)
2771 : END IF
2772 :
2773 14580 : ALLOCATE (v_drhoa(nspins), v_drhob(nspins))
2774 5290 : DO ispin = 1, nspins
2775 3290 : CALL allocate_pw(v_drhoa(ispin), pw_pool, bo)
2776 5290 : CALL allocate_pw(v_drhob(ispin), pw_pool, bo)
2777 : END DO
2778 :
2779 : END IF
2780 :
2781 4232 : IF (laplace_f) THEN
2782 38 : CALL xc_rho_set_get(rho1_set, laplace_rhoa=laplace1a, laplace_rhob=laplace1b)
2783 :
2784 190 : ALLOCATE (v_laplace(nspins))
2785 114 : DO ispin = 1, nspins
2786 114 : CALL allocate_pw(v_laplace(ispin), pw_pool, bo)
2787 : END DO
2788 :
2789 38 : IF (my_compute_virial) CALL xc_rho_set_get(rho_set, rhoa=rhoa, rhob=rhob)
2790 : END IF
2791 :
2792 4232 : IF (tau_f) THEN
2793 74 : CALL xc_rho_set_get(rho1_set, tau_a=tau1a, tau_b=tau1b)
2794 : END IF
2795 :
2796 4232 : IF (nspins /= 1) THEN
2797 :
2798 33242 : $:add_2nd_derivative_terms(arguments_openshell)
2799 :
2800 : ELSE
2801 :
2802 1038 : $:add_2nd_derivative_terms(arguments_triplet_outer, arguments_triplet_inner)
2803 :
2804 : END IF
2805 :
2806 4232 : IF (gradient_f) THEN
2807 :
2808 2000 : IF (my_compute_virial) THEN
2809 10 : CALL virial_drho_drho(virial_pw, drhoa, v_drhoa(1), virial_xc)
2810 10 : CALL virial_drho_drho(virial_pw, drhob, v_drhob(2), virial_xc)
2811 40 : DO idir = 1, 3
2812 30 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(drho,idir,v_drho,virial_pw)
2813 : virial_pw%array(:, :, :) = drho(idir)%array(:, :, :)*(v_drho(1)%array(:, :, :) + v_drho(2)%array(:, :, :))
2814 : !$OMP END PARALLEL WORKSHARE
2815 100 : DO jdir = 1, idir
2816 : tmp = -0.5_dp*virial_pw%pw_grid%dvol*accurate_dot_product(virial_pw%array(:, :, :), &
2817 60 : drho(jdir)%array(:, :, :))
2818 60 : virial_xc(jdir, idir) = virial_xc(jdir, idir) + tmp
2819 90 : virial_xc(idir, jdir) = virial_xc(jdir, idir)
2820 : END DO
2821 : END DO
2822 : END IF ! my_compute_virial
2823 :
2824 2000 : IF (my_gapw) THEN
2825 : !$OMP PARALLEL DO DEFAULT(NONE) &
2826 : !$OMP PRIVATE(ia,idir,ispin,ir) &
2827 : !$OMP SHARED(bo,nspins,vxg,drhoa,drhob,v_drhoa,v_drhob,v_drho, &
2828 : !$OMP e_drhoa,e_drhob,e_drho,drho1a,drho1b,fac,drho,drho1) &
2829 524 : !$OMP COLLAPSE(3)
2830 : DO ir = bo(1, 2), bo(2, 2)
2831 : DO ia = bo(1, 1), bo(2, 1)
2832 : DO idir = 1, 3
2833 : DO ispin = 1, nspins
2834 : vxg(idir, ia, ir, ispin) = &
2835 : -(v_drhoa(ispin)%array(ia, ir, 1)*drhoa(idir)%array(ia, ir, 1) + &
2836 : v_drhob(ispin)%array(ia, ir, 1)*drhob(idir)%array(ia, ir, 1) + &
2837 : v_drho(ispin)%array(ia, ir, 1)*drho(idir)%array(ia, ir, 1))
2838 : END DO
2839 : IF (ASSOCIATED(e_drhoa)) THEN
2840 : vxg(idir, ia, ir, 1) = vxg(idir, ia, ir, 1) + &
2841 : e_drhoa(ia, ir, 1)*drho1a(idir)%array(ia, ir, 1)
2842 : END IF
2843 : IF (nspins /= 1 .AND. ASSOCIATED(e_drhob)) THEN
2844 : vxg(idir, ia, ir, 2) = vxg(idir, ia, ir, 2) + &
2845 : e_drhob(ia, ir, 1)*drho1b(idir)%array(ia, ir, 1)
2846 : END IF
2847 : IF (ASSOCIATED(e_drho)) THEN
2848 : IF (nspins /= 1) THEN
2849 : vxg(idir, ia, ir, 1) = vxg(idir, ia, ir, 1) + &
2850 : e_drho(ia, ir, 1)*drho1(idir)%array(ia, ir, 1)
2851 : vxg(idir, ia, ir, 2) = vxg(idir, ia, ir, 2) + &
2852 : e_drho(ia, ir, 1)*drho1(idir)%array(ia, ir, 1)
2853 : ELSE
2854 : vxg(idir, ia, ir, 1) = vxg(idir, ia, ir, 1) + &
2855 : e_drho(ia, ir, 1)*(drho1a(idir)%array(ia, ir, 1) + &
2856 : fac*drho1b(idir)%array(ia, ir, 1))
2857 : END IF
2858 : END IF
2859 : END DO
2860 : END DO
2861 : END DO
2862 : !$OMP END PARALLEL DO
2863 : ELSE
2864 :
2865 : ! partial integration
2866 5904 : DO idir = 1, 3
2867 :
2868 12342 : DO ispin = 1, nspins
2869 12342 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(v_drho_r,v_drhoa,v_drhob,v_drho,drhoa,drhob,drho,ispin,idir)
2870 : v_drho_r(idir, ispin)%array(:, :, :) = &
2871 : v_drhoa(ispin)%array(:, :, :)*drhoa(idir)%array(:, :, :) + &
2872 : v_drhob(ispin)%array(:, :, :)*drhob(idir)%array(:, :, :) + &
2873 : v_drho(ispin)%array(:, :, :)*drho(idir)%array(:, :, :)
2874 : !$OMP END PARALLEL WORKSHARE
2875 : END DO
2876 4428 : IF (ASSOCIATED(e_drhoa)) THEN
2877 4428 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(v_drho_r,e_drhoa,drho1a,idir)
2878 : v_drho_r(idir, 1)%array(:, :, :) = v_drho_r(idir, 1)%array(:, :, :) - &
2879 : e_drhoa(:, :, :)*drho1a(idir)%array(:, :, :)
2880 : !$OMP END PARALLEL WORKSHARE
2881 : END IF
2882 4428 : IF (nspins /= 1 .AND. ASSOCIATED(e_drhob)) THEN
2883 3486 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(v_drho_r,e_drhob,drho1b,idir)
2884 : v_drho_r(idir, 2)%array(:, :, :) = v_drho_r(idir, 2)%array(:, :, :) - &
2885 : e_drhob(:, :, :)*drho1b(idir)%array(:, :, :)
2886 : !$OMP END PARALLEL WORKSHARE
2887 : END IF
2888 5904 : IF (ASSOCIATED(e_drho)) THEN
2889 : !$OMP PARALLEL DO DEFAULT(NONE) &
2890 : !$OMP PRIVATE(k,j,i) &
2891 : !$OMP SHARED(bo,v_drho_r,e_drho,drho1a,drho1b,drho1,fac,idir,nspins) &
2892 4428 : !$OMP COLLAPSE(3)
2893 : DO k = bo(1, 3), bo(2, 3)
2894 : DO j = bo(1, 2), bo(2, 2)
2895 : DO i = bo(1, 1), bo(2, 1)
2896 : IF (nspins /= 1) THEN
2897 : v_drho_r(idir, 1)%array(i, j, k) = v_drho_r(idir, 1)%array(i, j, k) - &
2898 : e_drho(i, j, k)*drho1(idir)%array(i, j, k)
2899 : v_drho_r(idir, 2)%array(i, j, k) = v_drho_r(idir, 2)%array(i, j, k) - &
2900 : e_drho(i, j, k)*drho1(idir)%array(i, j, k)
2901 : ELSE
2902 : v_drho_r(idir, 1)%array(i, j, k) = v_drho_r(idir, 1)%array(i, j, k) - &
2903 : e_drho(i, j, k)*(drho1a(idir)%array(i, j, k) + &
2904 : fac*drho1b(idir)%array(i, j, k))
2905 : END IF
2906 : END DO
2907 : END DO
2908 : END DO
2909 : !$OMP END PARALLEL DO
2910 : END IF
2911 : END DO
2912 :
2913 4114 : DO ispin = 1, nspins
2914 : ! partial integration
2915 4114 : CALL xc_pw_divergence(xc_deriv_method_id, v_drho_r(:, ispin), tmp_g, vxc_g, v_xc(ispin))
2916 : END DO ! ispin
2917 :
2918 : END IF
2919 :
2920 8000 : DO idir = 1, 3
2921 6000 : DEALLOCATE (drho(idir)%array)
2922 8000 : DEALLOCATE (drho1(idir)%array)
2923 : END DO
2924 :
2925 5290 : DO ispin = 1, nspins
2926 3290 : CALL deallocate_pw(v_drhoa(ispin), pw_pool)
2927 5290 : CALL deallocate_pw(v_drhob(ispin), pw_pool)
2928 : END DO
2929 :
2930 2000 : DEALLOCATE (v_drhoa, v_drhob)
2931 :
2932 : END IF ! gradient_f
2933 :
2934 4232 : IF (laplace_f .AND. my_compute_virial) THEN
2935 15026 : virial_pw%array(:, :, :) = -rhoa(:, :, :)
2936 2 : CALL virial_laplace(virial_pw, pw_pool, virial_xc, v_laplace(1)%array)
2937 15026 : virial_pw%array(:, :, :) = -rhob(:, :, :)
2938 2 : CALL virial_laplace(virial_pw, pw_pool, virial_xc, v_laplace(2)%array)
2939 : END IF
2940 :
2941 : ELSE
2942 :
2943 : !-----------------!
2944 : ! restricted case !
2945 : !-----------------!
2946 :
2947 26110 : CALL xc_rho_set_get(rho1_set, rho=rho1)
2948 :
2949 26110 : IF (gradient_f) THEN
2950 16942 : CALL xc_rho_set_get(rho_set, drho=drho, norm_drho=norm_drho)
2951 16942 : CALL xc_rho_set_get(rho1_set, drho=drho1)
2952 16942 : CALL prepare_dr1dr(dr1dr, drho, drho1)
2953 : END IF
2954 :
2955 26110 : IF (laplace_f) THEN
2956 136 : CALL xc_rho_set_get(rho1_set, laplace_rho=laplace1)
2957 :
2958 544 : ALLOCATE (v_laplace(nspins))
2959 272 : DO ispin = 1, nspins
2960 272 : CALL allocate_pw(v_laplace(ispin), pw_pool, bo)
2961 : END DO
2962 :
2963 136 : IF (my_compute_virial) CALL xc_rho_set_get(rho_set, rho=rho)
2964 : END IF
2965 :
2966 26110 : IF (tau_f) THEN
2967 582 : CALL xc_rho_set_get(rho1_set, tau=tau1)
2968 : END IF
2969 :
2970 320626 : $:add_2nd_derivative_terms(arguments_closedshell)
2971 :
2972 26110 : IF (gradient_f) THEN
2973 :
2974 16942 : IF (my_compute_virial) THEN
2975 222 : CALL virial_drho_drho(virial_pw, drho, v_drho(1), virial_xc)
2976 : END IF ! my_compute_virial
2977 :
2978 16942 : IF (my_gapw) THEN
2979 :
2980 26080 : DO idir = 1, 3
2981 : !$OMP PARALLEL DO DEFAULT(NONE) &
2982 : !$OMP PRIVATE(ia,ir) &
2983 : !$OMP SHARED(bo,vxg,drho,v_drho,e_drho,drho1,idir,factor2) &
2984 26080 : !$OMP COLLAPSE(2)
2985 : DO ia = bo(1, 1), bo(2, 1)
2986 : DO ir = bo(1, 2), bo(2, 2)
2987 : vxg(idir, ia, ir, 1) = -drho(idir)%array(ia, ir, 1)*v_drho(1)%array(ia, ir, 1)
2988 : IF (ASSOCIATED(e_drho)) THEN
2989 : vxg(idir, ia, ir, 1) = vxg(idir, ia, ir, 1) + factor2*drho1(idir)%array(ia, ir, 1)*e_drho(ia, ir, 1)
2990 : END IF
2991 : END DO
2992 : END DO
2993 : !$OMP END PARALLEL DO
2994 : END DO
2995 :
2996 : ELSE
2997 : ! partial integration
2998 41688 : DO idir = 1, 3
2999 41688 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(v_drho_r,drho,v_drho,drho1,e_drho,idir)
3000 : v_drho_r(idir, 1)%array(:, :, :) = drho(idir)%array(:, :, :)*v_drho(1)%array(:, :, :) - &
3001 : drho1(idir)%array(:, :, :)*e_drho(:, :, :)
3002 : !$OMP END PARALLEL WORKSHARE
3003 : END DO
3004 :
3005 10422 : CALL xc_pw_divergence(xc_deriv_method_id, v_drho_r(:, 1), tmp_g, vxc_g, v_xc(1))
3006 : END IF
3007 :
3008 : END IF
3009 :
3010 26110 : IF (laplace_f .AND. my_compute_virial) THEN
3011 294530 : virial_pw%array(:, :, :) = -rho(:, :, :)
3012 14 : CALL virial_laplace(virial_pw, pw_pool, virial_xc, v_laplace(1)%array)
3013 : END IF
3014 :
3015 : END IF
3016 :
3017 30342 : IF (laplace_f) THEN
3018 386 : DO ispin = 1, nspins
3019 212 : CALL xc_pw_laplace(v_laplace(ispin), pw_pool, xc_deriv_method_id)
3020 386 : CALL pw_axpy(v_laplace(ispin), v_xc(ispin))
3021 : END DO
3022 : END IF
3023 :
3024 30342 : IF (gradient_f) THEN
3025 :
3026 39174 : DO ispin = 1, nspins
3027 20232 : CALL deallocate_pw(v_drho(ispin), pw_pool)
3028 99870 : DO idir = 1, 3
3029 80928 : CALL deallocate_pw(v_drho_r(idir, ispin), pw_pool)
3030 : END DO
3031 : END DO
3032 18942 : DEALLOCATE (v_drho, v_drho_r)
3033 :
3034 : END IF
3035 :
3036 30342 : IF (laplace_f) THEN
3037 386 : DO ispin = 1, nspins
3038 386 : CALL deallocate_pw(v_laplace(ispin), pw_pool)
3039 : END DO
3040 174 : DEALLOCATE (v_laplace)
3041 : END IF
3042 :
3043 30342 : IF (ASSOCIATED(tmp_g%pw_grid) .AND. ASSOCIATED(pw_pool)) THEN
3044 11482 : CALL pw_pool%give_back_pw(tmp_g)
3045 : END IF
3046 :
3047 30342 : IF (ASSOCIATED(vxc_g%pw_grid) .AND. ASSOCIATED(pw_pool)) THEN
3048 11482 : CALL pw_pool%give_back_pw(vxc_g)
3049 : END IF
3050 :
3051 30342 : IF (my_compute_virial .AND. (gradient_f .OR. laplace_f)) THEN
3052 232 : CALL deallocate_pw(virial_pw, pw_pool)
3053 : END IF
3054 :
3055 30342 : CALL timestop(handle)
3056 :
3057 91026 : END SUBROUTINE xc_calc_2nd_deriv_analytical
3058 :
3059 : ! **************************************************************************************************
3060 : !> \brief allocates grids using pw_pool (if associated) or with bounds
3061 : !> \param pw ...
3062 : !> \param pw_pool ...
3063 : !> \param bo ...
3064 : ! **************************************************************************************************
3065 88050 : SUBROUTINE allocate_pw(pw, pw_pool, bo)
3066 : TYPE(pw_r3d_rs_type), INTENT(OUT) :: pw
3067 : TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool
3068 : INTEGER, DIMENSION(2, 3), INTENT(IN) :: bo
3069 :
3070 88050 : IF (ASSOCIATED(pw_pool)) THEN
3071 58058 : CALL pw_pool%create_pw(pw)
3072 58058 : CALL pw_zero(pw)
3073 : ELSE
3074 149960 : ALLOCATE (pw%array(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
3075 76539584 : pw%array = 0.0_dp
3076 : END IF
3077 :
3078 88050 : END SUBROUTINE allocate_pw
3079 :
3080 : ! **************************************************************************************************
3081 : !> \brief deallocates grid allocated with allocate_pw
3082 : !> \param pw ...
3083 : !> \param pw_pool ...
3084 : ! **************************************************************************************************
3085 88050 : SUBROUTINE deallocate_pw(pw, pw_pool)
3086 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: pw
3087 : TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool
3088 :
3089 88050 : IF (ASSOCIATED(pw_pool)) THEN
3090 58058 : CALL pw_pool%give_back_pw(pw)
3091 : ELSE
3092 29992 : CALL pw%release()
3093 : END IF
3094 :
3095 88050 : END SUBROUTINE deallocate_pw
3096 :
3097 : ! **************************************************************************************************
3098 : !> \brief updates virial from first derivative w.r.t. norm_drho
3099 : !> \param virial_pw ...
3100 : !> \param drho ...
3101 : !> \param drho1 ...
3102 : !> \param deriv_data ...
3103 : !> \param virial_xc ...
3104 : ! **************************************************************************************************
3105 304 : SUBROUTINE virial_drho_drho1(virial_pw, drho, drho1, deriv_data, virial_xc)
3106 : TYPE(pw_r3d_rs_type), INTENT(IN) :: virial_pw
3107 : TYPE(cp_3d_r_cp_type), DIMENSION(3), INTENT(IN) :: drho, drho1
3108 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: deriv_data
3109 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT) :: virial_xc
3110 :
3111 : INTEGER :: idir, jdir
3112 : REAL(KIND=dp) :: tmp
3113 :
3114 1216 : DO idir = 1, 3
3115 912 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(drho,idir,virial_pw,deriv_data)
3116 : virial_pw%array(:, :, :) = drho(idir)%array(:, :, :)*deriv_data(:, :, :)
3117 : !$OMP END PARALLEL WORKSHARE
3118 3952 : DO jdir = 1, 3
3119 : tmp = virial_pw%pw_grid%dvol*accurate_dot_product(virial_pw%array(:, :, :), &
3120 2736 : drho1(jdir)%array(:, :, :))
3121 2736 : virial_xc(jdir, idir) = virial_xc(jdir, idir) + tmp
3122 3648 : virial_xc(idir, jdir) = virial_xc(idir, jdir) + tmp
3123 : END DO
3124 : END DO
3125 :
3126 304 : END SUBROUTINE virial_drho_drho1
3127 :
3128 : ! **************************************************************************************************
3129 : !> \brief Adds virial contribution from second order potential parts
3130 : !> \param virial_pw ...
3131 : !> \param drho ...
3132 : !> \param v_drho ...
3133 : !> \param virial_xc ...
3134 : ! **************************************************************************************************
3135 298 : SUBROUTINE virial_drho_drho(virial_pw, drho, v_drho, virial_xc)
3136 : TYPE(pw_r3d_rs_type), INTENT(IN) :: virial_pw
3137 : TYPE(cp_3d_r_cp_type), DIMENSION(3), INTENT(IN) :: drho
3138 : TYPE(pw_r3d_rs_type), INTENT(IN) :: v_drho
3139 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT) :: virial_xc
3140 :
3141 : INTEGER :: idir, jdir
3142 : REAL(KIND=dp) :: tmp
3143 :
3144 1192 : DO idir = 1, 3
3145 894 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(drho,idir,v_drho,virial_pw)
3146 : virial_pw%array(:, :, :) = drho(idir)%array(:, :, :)*v_drho%array(:, :, :)
3147 : !$OMP END PARALLEL WORKSHARE
3148 2980 : DO jdir = 1, idir
3149 : tmp = -virial_pw%pw_grid%dvol*accurate_dot_product(virial_pw%array(:, :, :), &
3150 1788 : drho(jdir)%array(:, :, :))
3151 1788 : virial_xc(jdir, idir) = virial_xc(jdir, idir) + tmp
3152 2682 : virial_xc(idir, jdir) = virial_xc(jdir, idir)
3153 : END DO
3154 : END DO
3155 :
3156 298 : END SUBROUTINE virial_drho_drho
3157 :
3158 : ! **************************************************************************************************
3159 : !> \brief ...
3160 : !> \param rho_r ...
3161 : !> \param pw_pool ...
3162 : !> \param virial_xc ...
3163 : !> \param deriv_data ...
3164 : ! **************************************************************************************************
3165 150 : SUBROUTINE virial_laplace(rho_r, pw_pool, virial_xc, deriv_data)
3166 : TYPE(pw_r3d_rs_type), TARGET :: rho_r
3167 : TYPE(pw_pool_type), POINTER, INTENT(IN) :: pw_pool
3168 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT) :: virial_xc
3169 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: deriv_data
3170 :
3171 : CHARACTER(len=*), PARAMETER :: routineN = 'virial_laplace'
3172 :
3173 : INTEGER :: handle, idir, jdir
3174 : TYPE(pw_r3d_rs_type), POINTER :: virial_pw
3175 : TYPE(pw_c1d_gs_type), POINTER :: tmp_g, rho_g
3176 : INTEGER, DIMENSION(3) :: my_deriv
3177 :
3178 150 : CALL timeset(routineN, handle)
3179 :
3180 150 : NULLIFY (virial_pw, tmp_g, rho_g)
3181 150 : ALLOCATE (virial_pw, tmp_g, rho_g)
3182 150 : CALL pw_pool%create_pw(virial_pw)
3183 150 : CALL pw_pool%create_pw(tmp_g)
3184 150 : CALL pw_pool%create_pw(rho_g)
3185 150 : CALL pw_zero(virial_pw)
3186 150 : CALL pw_transfer(rho_r, rho_g)
3187 600 : DO idir = 1, 3
3188 1500 : DO jdir = idir, 3
3189 900 : CALL pw_copy(rho_g, tmp_g)
3190 :
3191 900 : my_deriv = 0
3192 900 : my_deriv(idir) = 1
3193 900 : my_deriv(jdir) = my_deriv(jdir) + 1
3194 :
3195 900 : CALL pw_derive(tmp_g, my_deriv)
3196 900 : CALL pw_transfer(tmp_g, virial_pw)
3197 : virial_xc(idir, jdir) = virial_xc(idir, jdir) - 2.0_dp*virial_pw%pw_grid%dvol* &
3198 : accurate_dot_product(virial_pw%array(:, :, :), &
3199 900 : deriv_data(:, :, :))
3200 1350 : virial_xc(jdir, idir) = virial_xc(idir, jdir)
3201 : END DO
3202 : END DO
3203 150 : CALL pw_pool%give_back_pw(virial_pw)
3204 150 : CALL pw_pool%give_back_pw(tmp_g)
3205 150 : CALL pw_pool%give_back_pw(rho_g)
3206 150 : DEALLOCATE (virial_pw, tmp_g, rho_g)
3207 :
3208 150 : CALL timestop(handle)
3209 :
3210 150 : END SUBROUTINE virial_laplace
3211 :
3212 : ! **************************************************************************************************
3213 : !> \brief Prepare objects for the calculation of the 2nd derivatives of the density functional.
3214 : !> The calculation must then be performed with xc_calc_2nd_deriv.
3215 : !> \param deriv_set object containing the XC derivatives (out)
3216 : !> \param rho_set object that will contain the density at which the
3217 : !> derivatives were calculated
3218 : !> \param rho_r the place where you evaluate the derivative
3219 : !> \param pw_pool the pool for the grids
3220 : !> \param xc_section which functional should be used and how to calculate it
3221 : !> \param tau_r kinetic energy density in real space
3222 : ! **************************************************************************************************
3223 12114 : SUBROUTINE xc_prep_2nd_deriv(deriv_set, &
3224 : rho_set, rho_r, pw_pool, xc_section, tau_r)
3225 :
3226 : TYPE(xc_derivative_set_type) :: deriv_set
3227 : TYPE(xc_rho_set_type) :: rho_set
3228 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
3229 : TYPE(pw_pool_type), POINTER :: pw_pool
3230 : TYPE(section_vals_type), POINTER :: xc_section
3231 : TYPE(pw_r3d_rs_type), DIMENSION(:), OPTIONAL, POINTER :: tau_r
3232 :
3233 : CHARACTER(len=*), PARAMETER :: routineN = 'xc_prep_2nd_deriv'
3234 :
3235 : INTEGER :: handle, nspins
3236 : LOGICAL :: lsd
3237 12114 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
3238 12114 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: tau
3239 :
3240 12114 : CALL timeset(routineN, handle)
3241 :
3242 12114 : CPASSERT(ASSOCIATED(xc_section))
3243 12114 : CPASSERT(ASSOCIATED(pw_pool))
3244 :
3245 12114 : nspins = SIZE(rho_r)
3246 12114 : lsd = (nspins /= 1)
3247 :
3248 12114 : NULLIFY (rho_g, tau)
3249 12114 : IF (PRESENT(tau_r)) &
3250 11452 : tau => tau_r
3251 :
3252 12114 : IF (section_get_lval(xc_section, "2ND_DERIV_ANALYTICAL")) THEN
3253 : CALL xc_rho_set_and_dset_create(rho_set, deriv_set, 2, &
3254 : rho_r, rho_g, tau, xc_section, pw_pool, &
3255 12030 : calc_potential=.TRUE.)
3256 : ELSE
3257 : CALL xc_rho_set_and_dset_create(rho_set, deriv_set, 1, &
3258 : rho_r, rho_g, tau, xc_section, pw_pool, &
3259 84 : calc_potential=.TRUE.)
3260 : END IF
3261 :
3262 12114 : CALL timestop(handle)
3263 :
3264 12114 : END SUBROUTINE xc_prep_2nd_deriv
3265 :
3266 : ! **************************************************************************************************
3267 : !> \brief divides derivatives from deriv_set by norm_drho
3268 : !> \param deriv_set ...
3269 : !> \param rho_set ...
3270 : !> \param lsd ...
3271 : ! **************************************************************************************************
3272 143367 : SUBROUTINE divide_by_norm_drho(deriv_set, rho_set, lsd)
3273 :
3274 : TYPE(xc_derivative_set_type), INTENT(INOUT) :: deriv_set
3275 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
3276 : LOGICAL, INTENT(IN) :: lsd
3277 :
3278 143367 : INTEGER, DIMENSION(:), POINTER :: split_desc
3279 : INTEGER :: idesc
3280 : INTEGER, DIMENSION(2, 3) :: bo
3281 : REAL(KIND=dp) :: drho_cutoff
3282 143367 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: norm_drho, norm_drhoa, norm_drhob
3283 1720404 : TYPE(cp_3d_r_cp_type), DIMENSION(3) :: drho, drhoa, drhob
3284 : TYPE(cp_sll_xc_deriv_type), POINTER :: pos
3285 : TYPE(xc_derivative_type), POINTER :: deriv_att
3286 :
3287 : ! check for unknown derivatives and divide by norm_drho where necessary
3288 :
3289 1433670 : bo = rho_set%local_bounds
3290 : CALL xc_rho_set_get(rho_set, drho_cutoff=drho_cutoff, norm_drho=norm_drho, &
3291 : norm_drhoa=norm_drhoa, norm_drhob=norm_drhob, &
3292 143367 : drho=drho, drhoa=drhoa, drhob=drhob, can_return_null=.TRUE.)
3293 :
3294 143367 : pos => deriv_set%derivs
3295 607947 : DO WHILE (cp_sll_xc_deriv_next(pos, el_att=deriv_att))
3296 464580 : CALL xc_derivative_get(deriv_att, split_desc=split_desc)
3297 989406 : DO idesc = 1, SIZE(split_desc)
3298 464580 : SELECT CASE (split_desc(idesc))
3299 : CASE (deriv_norm_drho)
3300 115657 : IF (ASSOCIATED(norm_drho)) THEN
3301 115657 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(deriv_att,norm_drho,drho_cutoff)
3302 : deriv_att%deriv_data(:, :, :) = deriv_att%deriv_data(:, :, :)/ &
3303 : MAX(norm_drho(:, :, :), drho_cutoff)
3304 : !$OMP END PARALLEL WORKSHARE
3305 0 : ELSE IF (ASSOCIATED(drho(1)%array)) THEN
3306 0 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(deriv_att,drho,drho_cutoff)
3307 : deriv_att%deriv_data(:, :, :) = deriv_att%deriv_data(:, :, :)/ &
3308 : MAX(SQRT(drho(1)%array(:, :, :)**2 + &
3309 : drho(2)%array(:, :, :)**2 + &
3310 : drho(3)%array(:, :, :)**2), drho_cutoff)
3311 : !$OMP END PARALLEL WORKSHARE
3312 0 : ELSE IF (ASSOCIATED(drhoa(1)%array) .AND. ASSOCIATED(drhob(1)%array)) THEN
3313 0 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(deriv_att,drhoa,drhob,drho_cutoff)
3314 : deriv_att%deriv_data(:, :, :) = deriv_att%deriv_data(:, :, :)/ &
3315 : MAX(SQRT((drhoa(1)%array(:, :, :) + drhob(1)%array(:, :, :))**2 + &
3316 : (drhoa(2)%array(:, :, :) + drhob(2)%array(:, :, :))**2 + &
3317 : (drhoa(3)%array(:, :, :) + drhob(3)%array(:, :, :))**2), drho_cutoff)
3318 : !$OMP END PARALLEL WORKSHARE
3319 : ELSE
3320 0 : CPABORT("Normalization of derivative requires any of norm_drho, drho or drhoa+drhob!")
3321 : END IF
3322 : CASE (deriv_norm_drhoa)
3323 19150 : IF (ASSOCIATED(norm_drhoa)) THEN
3324 19150 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(deriv_att,norm_drhoa,drho_cutoff)
3325 : deriv_att%deriv_data(:, :, :) = deriv_att%deriv_data(:, :, :)/ &
3326 : MAX(norm_drhoa(:, :, :), drho_cutoff)
3327 : !$OMP END PARALLEL WORKSHARE
3328 0 : ELSE IF (ASSOCIATED(drhoa(1)%array)) THEN
3329 0 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(deriv_att,drhoa,drho_cutoff)
3330 : deriv_att%deriv_data(:, :, :) = deriv_att%deriv_data(:, :, :)/ &
3331 : MAX(SQRT(drhoa(1)%array(:, :, :)**2 + &
3332 : drhoa(2)%array(:, :, :)**2 + &
3333 : drhoa(3)%array(:, :, :)**2), drho_cutoff)
3334 : !$OMP END PARALLEL WORKSHARE
3335 : ELSE
3336 0 : CPABORT("Normalization of derivative requires any of norm_drhoa or drhoa!")
3337 : END IF
3338 : CASE (deriv_norm_drhob)
3339 19146 : IF (ASSOCIATED(norm_drhob)) THEN
3340 19146 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(deriv_att,norm_drhob,drho_cutoff)
3341 : deriv_att%deriv_data(:, :, :) = deriv_att%deriv_data(:, :, :)/ &
3342 : MAX(norm_drhob(:, :, :), drho_cutoff)
3343 : !$OMP END PARALLEL WORKSHARE
3344 0 : ELSE IF (ASSOCIATED(drhob(1)%array)) THEN
3345 0 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(deriv_att,drhob,drho_cutoff)
3346 : deriv_att%deriv_data(:, :, :) = deriv_att%deriv_data(:, :, :)/ &
3347 : MAX(SQRT(drhob(1)%array(:, :, :)**2 + &
3348 : drhob(2)%array(:, :, :)**2 + &
3349 : drhob(3)%array(:, :, :)**2), drho_cutoff)
3350 : !$OMP END PARALLEL WORKSHARE
3351 : ELSE
3352 0 : CPABORT("Normalization of derivative requires any of norm_drhob or drhob!")
3353 : END IF
3354 : CASE (deriv_rho, deriv_tau, deriv_laplace_rho)
3355 161030 : IF (lsd) &
3356 0 : CPABORT(TRIM(id_to_desc(split_desc(idesc)))//" not handled in lsd!'")
3357 : CASE (deriv_rhoa, deriv_rhob, deriv_tau_a, deriv_tau_b, deriv_laplace_rhoa, deriv_laplace_rhob)
3358 : CASE default
3359 381459 : CPABORT("Unknown derivative id")
3360 : END SELECT
3361 : END DO
3362 : END DO
3363 :
3364 143367 : END SUBROUTINE divide_by_norm_drho
3365 :
3366 : ! **************************************************************************************************
3367 : !> \brief allocates and calculates drho from given spin densities drhoa, drhob
3368 : !> \param drho ...
3369 : !> \param drhoa ...
3370 : !> \param drhob ...
3371 : ! **************************************************************************************************
3372 16080 : SUBROUTINE calc_drho_from_ab(drho, drhoa, drhob)
3373 : TYPE(cp_3d_r_cp_type), DIMENSION(3), INTENT(OUT) :: drho
3374 : TYPE(cp_3d_r_cp_type), DIMENSION(3), INTENT(IN) :: drhoa, drhob
3375 :
3376 : CHARACTER(len=*), PARAMETER :: routineN = 'calc_drho_from_ab'
3377 :
3378 : INTEGER :: handle, idir
3379 :
3380 4020 : CALL timeset(routineN, handle)
3381 :
3382 16080 : DO idir = 1, 3
3383 12060 : NULLIFY (drho(idir)%array)
3384 : ALLOCATE (drho(idir)%array(LBOUND(drhoa(1)%array, 1):UBOUND(drhoa(1)%array, 1), &
3385 : LBOUND(drhoa(1)%array, 2):UBOUND(drhoa(1)%array, 2), &
3386 144720 : LBOUND(drhoa(1)%array, 3):UBOUND(drhoa(1)%array, 3)))
3387 16080 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(drho,drhoa,drhob,idir)
3388 : drho(idir)%array(:, :, :) = drhoa(idir)%array(:, :, :) + drhob(idir)%array(:, :, :)
3389 : !$OMP END PARALLEL WORKSHARE
3390 : END DO
3391 :
3392 4020 : CALL timestop(handle)
3393 :
3394 4020 : END SUBROUTINE calc_drho_from_ab
3395 :
3396 : ! **************************************************************************************************
3397 : !> \brief allocates and calculates dot products of two density gradients
3398 : !> \param dr1dr ...
3399 : !> \param drho ...
3400 : !> \param drho1 ...
3401 : ! **************************************************************************************************
3402 22288 : SUBROUTINE prepare_dr1dr(dr1dr, drho, drho1)
3403 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
3404 : INTENT(OUT) :: dr1dr
3405 : TYPE(cp_3d_r_cp_type), DIMENSION(3), INTENT(IN) :: drho, drho1
3406 :
3407 : CHARACTER(len=*), PARAMETER :: routineN = 'prepare_dr1dr'
3408 :
3409 : INTEGER :: handle, idir
3410 :
3411 22288 : CALL timeset(routineN, handle)
3412 :
3413 0 : ALLOCATE (dr1dr(LBOUND(drho(1)%array, 1):UBOUND(drho(1)%array, 1), &
3414 : LBOUND(drho(1)%array, 2):UBOUND(drho(1)%array, 2), &
3415 267456 : LBOUND(drho(1)%array, 3):UBOUND(drho(1)%array, 3)))
3416 :
3417 22288 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(dr1dr,drho,drho1)
3418 : dr1dr(:, :, :) = drho(1)%array(:, :, :)*drho1(1)%array(:, :, :)
3419 : !$OMP END PARALLEL WORKSHARE
3420 66864 : DO idir = 2, 3
3421 66864 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(dr1dr,drho,drho1,idir)
3422 : dr1dr(:, :, :) = dr1dr(:, :, :) + drho(idir)%array(:, :, :)*drho1(idir)%array(:, :, :)
3423 : !$OMP END PARALLEL WORKSHARE
3424 : END DO
3425 :
3426 22288 : CALL timestop(handle)
3427 :
3428 22288 : END SUBROUTINE prepare_dr1dr
3429 :
3430 : ! **************************************************************************************************
3431 : !> \brief allocates and calculates dot product of two densities for triplets
3432 : !> \param dr1dr ...
3433 : !> \param drhoa ...
3434 : !> \param drhob ...
3435 : !> \param drho1a ...
3436 : !> \param drho1b ...
3437 : !> \param fac ...
3438 : ! **************************************************************************************************
3439 710 : SUBROUTINE prepare_dr1dr_ab(dr1dr, drhoa, drhob, drho1a, drho1b, fac)
3440 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
3441 : INTENT(OUT) :: dr1dr
3442 : TYPE(cp_3d_r_cp_type), DIMENSION(3), INTENT(IN) :: drhoa, drhob, drho1a, drho1b
3443 : REAL(KIND=dp), INTENT(IN) :: fac
3444 :
3445 : CHARACTER(len=*), PARAMETER :: routineN = 'prepare_dr1dr_ab'
3446 :
3447 : INTEGER :: handle, idir
3448 :
3449 710 : CALL timeset(routineN, handle)
3450 :
3451 0 : ALLOCATE (dr1dr(LBOUND(drhoa(1)%array, 1):UBOUND(drhoa(1)%array, 1), &
3452 : LBOUND(drhoa(1)%array, 2):UBOUND(drhoa(1)%array, 2), &
3453 8520 : LBOUND(drhoa(1)%array, 3):UBOUND(drhoa(1)%array, 3)))
3454 :
3455 710 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(fac,dr1dr,drho1a,drho1b,drhoa,drhob)
3456 : dr1dr(:, :, :) = drhoa(1)%array(:, :, :)*(drho1a(1)%array(:, :, :) + &
3457 : fac*drho1b(1)%array(:, :, :)) + &
3458 : drhob(1)%array(:, :, :)*(fac*drho1a(1)%array(:, :, :) + &
3459 : drho1b(1)%array(:, :, :))
3460 : !$OMP END PARALLEL WORKSHARE
3461 2130 : DO idir = 2, 3
3462 2130 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(fac,dr1dr,drho1a,drho1b,drhoa,drhob,idir)
3463 : dr1dr(:, :, :) = dr1dr(:, :, :) + &
3464 : drhoa(idir)%array(:, :, :)*(drho1a(idir)%array(:, :, :) + &
3465 : fac*drho1b(idir)%array(:, :, :)) + &
3466 : drhob(idir)%array(:, :, :)*(fac*drho1a(idir)%array(:, :, :) + &
3467 : drho1b(idir)%array(:, :, :))
3468 : !$OMP END PARALLEL WORKSHARE
3469 : END DO
3470 :
3471 710 : CALL timestop(handle)
3472 :
3473 710 : END SUBROUTINE prepare_dr1dr_ab
3474 :
3475 : ! **************************************************************************************************
3476 : !> \brief checks for gradients
3477 : !> \param deriv_set ...
3478 : !> \param lsd ...
3479 : !> \param gradient_f ...
3480 : !> \param tau_f ...
3481 : !> \param laplace_f ...
3482 : ! **************************************************************************************************
3483 141251 : SUBROUTINE check_for_derivatives(deriv_set, lsd, rho_f, gradient_f, tau_f, laplace_f)
3484 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
3485 : LOGICAL, INTENT(IN) :: lsd
3486 : LOGICAL, INTENT(OUT) :: rho_f, gradient_f, tau_f, laplace_f
3487 :
3488 : CHARACTER(len=*), PARAMETER :: routineN = 'check_for_derivatives'
3489 :
3490 : INTEGER :: handle, iorder, order
3491 141251 : INTEGER, DIMENSION(:), POINTER :: split_desc
3492 : TYPE(cp_sll_xc_deriv_type), POINTER :: pos
3493 : TYPE(xc_derivative_type), POINTER :: deriv_att
3494 :
3495 141251 : CALL timeset(routineN, handle)
3496 :
3497 141251 : rho_f = .FALSE.
3498 141251 : gradient_f = .FALSE.
3499 141251 : tau_f = .FALSE.
3500 141251 : laplace_f = .FALSE.
3501 : ! check for unknown derivatives
3502 141251 : pos => deriv_set%derivs
3503 641379 : DO WHILE (cp_sll_xc_deriv_next(pos, el_att=deriv_att))
3504 : CALL xc_derivative_get(deriv_att, order=order, &
3505 500128 : split_desc=split_desc)
3506 641379 : IF (lsd) THEN
3507 300269 : DO iorder = 1, size(split_desc)
3508 149813 : SELECT CASE (split_desc(iorder))
3509 : CASE (deriv_rhoa, deriv_rhob)
3510 77768 : rho_f = .TRUE.
3511 : CASE (deriv_norm_drho, deriv_norm_drhoa, deriv_norm_drhob)
3512 68616 : gradient_f = .TRUE.
3513 : CASE (deriv_tau_a, deriv_tau_b)
3514 2760 : tau_f = .TRUE.
3515 : CASE (deriv_laplace_rhoa, deriv_laplace_rhob)
3516 1312 : laplace_f = .TRUE.
3517 : CASE (deriv_rho, deriv_tau, deriv_laplace_rho)
3518 0 : CPABORT("Derivative not handled in lsd!")
3519 : CASE default
3520 150456 : CPABORT("Unknown derivative id")
3521 : END SELECT
3522 : END DO
3523 : ELSE
3524 647980 : DO iorder = 1, size(split_desc)
3525 350315 : SELECT CASE (split_desc(iorder))
3526 : CASE (deriv_rho)
3527 176258 : rho_f = .TRUE.
3528 : CASE (deriv_tau)
3529 4948 : tau_f = .TRUE.
3530 : CASE (deriv_norm_drho)
3531 115021 : gradient_f = .TRUE.
3532 : CASE (deriv_laplace_rho)
3533 1438 : laplace_f = .TRUE.
3534 : CASE default
3535 297665 : CPABORT("Unknown derivative id")
3536 : END SELECT
3537 : END DO
3538 : END IF
3539 : END DO
3540 :
3541 141251 : CALL timestop(handle)
3542 :
3543 141251 : END SUBROUTINE check_for_derivatives
3544 :
3545 : END MODULE xc
|
