Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Exchange and Correlation functional calculations
10 : !> \par History
11 : !> (13-Feb-2001) JGH, based on earlier version of apsi
12 : !> 02.2003 Many many changes [fawzi]
13 : !> 03.2004 new xc interface [fawzi]
14 : !> 04.2004 kinetic functionals [fawzi]
15 : !> \author fawzi
16 : ! **************************************************************************************************
17 : MODULE xc
18 : #:include 'xc.fypp'
19 : USE cp_array_utils, ONLY: cp_3d_r_cp_type
20 : USE cp_linked_list_xc_deriv, ONLY: cp_sll_xc_deriv_next, &
21 : cp_sll_xc_deriv_type
22 : USE cp_log_handling, ONLY: cp_get_default_logger, &
23 : cp_logger_get_default_unit_nr, &
24 : cp_logger_type, &
25 : cp_to_string
26 : USE input_section_types, ONLY: section_get_ival, &
27 : section_get_lval, &
28 : section_get_rval, &
29 : section_vals_get_subs_vals, &
30 : section_vals_type, &
31 : section_vals_val_get
32 : USE kahan_sum, ONLY: accurate_dot_product, &
33 : accurate_sum
34 : USE kinds, ONLY: default_path_length, &
35 : dp
36 : USE pw_grid_types, ONLY: PW_MODE_DISTRIBUTED, &
37 : pw_grid_type
38 : USE pw_methods, ONLY: pw_axpy, &
39 : pw_copy, &
40 : pw_derive, &
41 : pw_scale, &
42 : pw_transfer, &
43 : pw_zero, pw_integrate_function
44 : USE pw_pool_types, ONLY: &
45 : pw_pool_type
46 : USE pw_types, ONLY: &
47 : pw_c1d_gs_type, pw_r3d_rs_type
48 : USE xc_derivative_desc, ONLY: &
49 : deriv_rho, deriv_rhoa, deriv_rhob, &
50 : deriv_norm_drhoa, deriv_norm_drhob, deriv_norm_drho, deriv_tau_a, deriv_tau_b, deriv_tau, &
51 : deriv_laplace_rho, deriv_laplace_rhoa, deriv_laplace_rhob, id_to_desc
52 : USE xc_derivative_set_types, ONLY: xc_derivative_set_type, &
53 : xc_dset_create, &
54 : xc_dset_get_derivative, &
55 : xc_dset_release, &
56 : xc_dset_zero_all, xc_dset_recover_pw
57 : USE xc_derivative_types, ONLY: xc_derivative_get, &
58 : xc_derivative_type
59 : USE xc_derivatives, ONLY: xc_functionals_eval, &
60 : xc_functionals_get_needs
61 : USE xc_input_constants, ONLY: &
62 : xc_debug_new_routine, xc_default_f_routine, xc_test_lsd_f_routine
63 : USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
64 : USE xc_rho_set_types, ONLY: xc_rho_set_create, &
65 : xc_rho_set_get, &
66 : xc_rho_set_release, &
67 : xc_rho_set_type, &
68 : xc_rho_set_update, xc_rho_set_recover_pw
69 : USE xc_util, ONLY: xc_pw_smooth, xc_pw_laplace, xc_pw_divergence, xc_requires_tmp_g
70 : #include "../base/base_uses.f90"
71 :
72 : IMPLICIT NONE
73 : PRIVATE
74 : PUBLIC :: xc_vxc_pw_create1, xc_vxc_pw_create, &
75 : xc_exc_calc, xc_calc_2nd_deriv_analytical, xc_calc_2nd_deriv_numerical, &
76 : xc_calc_2nd_deriv, xc_prep_2nd_deriv, divide_by_norm_drho, smooth_cutoff, &
77 : xc_uses_kinetic_energy_density, xc_uses_norm_drho
78 : PUBLIC :: calc_xc_density
79 :
80 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
81 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc'
82 :
83 : CONTAINS
84 :
85 : ! **************************************************************************************************
86 : !> \brief ...
87 : !> \param xc_fun_section ...
88 : !> \param lsd ...
89 : !> \return ...
90 : ! **************************************************************************************************
91 6766 : 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 13532 : calc_potential=.FALSE.)
101 6766 : res = (needs%tau_spin .OR. needs%tau)
102 :
103 6766 : END FUNCTION
104 :
105 : ! **************************************************************************************************
106 : !> \brief ...
107 : !> \param xc_fun_section ...
108 : !> \param lsd ...
109 : !> \return ...
110 : ! **************************************************************************************************
111 6782 : 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 13564 : calc_potential=.FALSE.)
121 6782 : res = (needs%norm_drho .OR. needs%norm_drho_spin)
122 :
123 6782 : 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 207832 : 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 103916 : CPASSERT(ASSOCIATED(rho_r))
164 103916 : CPASSERT(ASSOCIATED(xc_section))
165 103916 : CPASSERT(.NOT. ASSOCIATED(vxc_rho))
166 103916 : CPASSERT(.NOT. ASSOCIATED(vxc_tau))
167 :
168 103916 : virial_xc = 0.0_dp
169 :
170 : CALL section_vals_val_get(xc_section, "FUNCTIONAL_ROUTINE", &
171 103916 : i_val=f_routine)
172 103916 : 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 103916 : 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 103916 : xc_section=xc_section, pw_pool=pw_pool)
187 : CASE default
188 : END SELECT
189 :
190 103916 : 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 128281 : 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 128281 : CALL timeset(routineN, handle)
859 :
860 128281 : CPASSERT(ASSOCIATED(pw_pool))
861 :
862 128281 : nspins = SIZE(rho_r)
863 128281 : lsd = (nspins /= 1)
864 :
865 128281 : xc_fun_sections => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
866 :
867 128281 : 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 128281 : 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 128281 : 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 128281 : deriv_order=deriv_order)
886 :
887 128281 : CALL divide_by_norm_drho(deriv_set, rho_set, lsd)
888 :
889 128281 : CALL timestop(handle)
890 :
891 128281 : 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 247514 : 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 247514 : CPASSERT(ASSOCIATED(pot))
924 990056 : bo(1, :) = LBOUND(pot)
925 990056 : bo(2, :) = UBOUND(pot)
926 247514 : my_e_0_scale_factor = 1.0_dp
927 247514 : IF (PRESENT(e_0_scale_factor)) my_e_0_scale_factor = e_0_scale_factor
928 247514 : rho_smooth_cutoff = rho_cutoff*rho_smooth_cutoff_range
929 247514 : rho_smooth_cutoff_2 = (rho_cutoff + rho_smooth_cutoff)/2
930 247514 : rho_smooth_cutoff_range_2 = rho_smooth_cutoff_2 - rho_cutoff
931 :
932 247514 : 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) SHARED(bo,e_0,pot,rho,&
937 : !$OMP rho_cutoff,rho_smooth_cutoff,rho_smooth_cutoff_2,&
938 : !$OMP rho_smooth_cutoff_range_2,my_e_0_scale_factor)&
939 0 : !$OMP PRIVATE(k,j,i,my_rho,my_rho_n,my_rho_n2) COLLAPSE(3)
940 : DO k = bo(1, 3), bo(2, 3)
941 : DO j = bo(1, 2), bo(2, 2)
942 : DO i = bo(1, 1), bo(2, 1)
943 : my_rho = rho(i, j, k)
944 : IF (my_rho < rho_smooth_cutoff) THEN
945 : IF (my_rho < rho_cutoff) THEN
946 : pot(i, j, k) = 0.0_dp
947 : ELSEIF (my_rho < rho_smooth_cutoff_2) THEN
948 : my_rho_n = (my_rho - rho_cutoff)/rho_smooth_cutoff_range_2
949 : my_rho_n2 = my_rho_n*my_rho_n
950 : pot(i, j, k) = pot(i, j, k)* &
951 : my_rho_n2*(my_rho_n - 0.5_dp*my_rho_n2) + &
952 : my_e_0_scale_factor*e_0(i, j, k)* &
953 : my_rho_n2*(3.0_dp - 2.0_dp*my_rho_n) &
954 : /rho_smooth_cutoff_range_2
955 : ELSE
956 : my_rho_n = 2.0_dp - (my_rho - rho_cutoff)/rho_smooth_cutoff_range_2
957 : my_rho_n2 = my_rho_n*my_rho_n
958 : pot(i, j, k) = pot(i, j, k)* &
959 : (1.0_dp - my_rho_n2*(my_rho_n - 0.5_dp*my_rho_n2)) &
960 : + my_e_0_scale_factor*e_0(i, j, k)* &
961 : my_rho_n2*(3.0_dp - 2.0_dp*my_rho_n) &
962 : /rho_smooth_cutoff_range_2
963 : END IF
964 : END IF
965 : END DO
966 : END DO
967 : END DO
968 : ELSE
969 : !$OMP PARALLEL DO DEFAULT(NONE) SHARED(bo,pot,e_0,rhoa,rhob,&
970 : !$OMP rho_cutoff,rho_smooth_cutoff,rho_smooth_cutoff_2,&
971 : !$OMP rho_smooth_cutoff_range_2,my_e_0_scale_factor)&
972 0 : !$OMP PRIVATE(k,j,i,my_rho,my_rho_n,my_rho_n2) COLLAPSE(3)
973 : DO k = bo(1, 3), bo(2, 3)
974 : DO j = bo(1, 2), bo(2, 2)
975 : DO i = bo(1, 1), bo(2, 1)
976 : my_rho = rhoa(i, j, k) + rhob(i, j, k)
977 : IF (my_rho < rho_smooth_cutoff) THEN
978 : IF (my_rho < rho_cutoff) THEN
979 : pot(i, j, k) = 0.0_dp
980 : ELSEIF (my_rho < rho_smooth_cutoff_2) THEN
981 : my_rho_n = (my_rho - rho_cutoff)/rho_smooth_cutoff_range_2
982 : my_rho_n2 = my_rho_n*my_rho_n
983 : pot(i, j, k) = pot(i, j, k)* &
984 : my_rho_n2*(my_rho_n - 0.5_dp*my_rho_n2) + &
985 : my_e_0_scale_factor*e_0(i, j, k)* &
986 : my_rho_n2*(3.0_dp - 2.0_dp*my_rho_n) &
987 : /rho_smooth_cutoff_range_2
988 : ELSE
989 : my_rho_n = 2.0_dp - (my_rho - rho_cutoff)/rho_smooth_cutoff_range_2
990 : my_rho_n2 = my_rho_n*my_rho_n
991 : pot(i, j, k) = pot(i, j, k)* &
992 : (1.0_dp - my_rho_n2*(my_rho_n - 0.5_dp*my_rho_n2)) &
993 : + my_e_0_scale_factor*e_0(i, j, k)* &
994 : my_rho_n2*(3.0_dp - 2.0_dp*my_rho_n) &
995 : /rho_smooth_cutoff_range_2
996 : END IF
997 : END IF
998 : END DO
999 : END DO
1000 : END DO
1001 : END IF
1002 : ELSE
1003 2 : IF (ASSOCIATED(rho)) THEN
1004 : !$OMP PARALLEL DO DEFAULT(NONE) SHARED(bo,pot,&
1005 : !$OMP rho_cutoff,rho_smooth_cutoff,rho_smooth_cutoff_2,&
1006 : !$OMP rho_smooth_cutoff_range_2,rho)&
1007 2 : !$OMP PRIVATE(k,j,i,my_rho,my_rho_n,my_rho_n2) COLLAPSE(3)
1008 : DO k = bo(1, 3), bo(2, 3)
1009 : DO j = bo(1, 2), bo(2, 2)
1010 : DO i = bo(1, 1), bo(2, 1)
1011 : my_rho = rho(i, j, k)
1012 : IF (my_rho < rho_smooth_cutoff) THEN
1013 : IF (my_rho < rho_cutoff) THEN
1014 : pot(i, j, k) = 0.0_dp
1015 : ELSEIF (my_rho < rho_smooth_cutoff_2) THEN
1016 : my_rho_n = (my_rho - rho_cutoff)/rho_smooth_cutoff_range_2
1017 : my_rho_n2 = my_rho_n*my_rho_n
1018 : pot(i, j, k) = pot(i, j, k)* &
1019 : my_rho_n2*(my_rho_n - 0.5_dp*my_rho_n2)
1020 : ELSE
1021 : my_rho_n = 2.0_dp - (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 : (1.0_dp - my_rho_n2*(my_rho_n - 0.5_dp*my_rho_n2))
1025 : END IF
1026 : END IF
1027 : END DO
1028 : END DO
1029 : END DO
1030 : ELSE
1031 : !$OMP PARALLEL DO DEFAULT(NONE) SHARED(bo,pot,&
1032 : !$OMP rho_cutoff,rho_smooth_cutoff,rho_smooth_cutoff_2,&
1033 : !$OMP rho_smooth_cutoff_range_2,rhoa,rhob)&
1034 0 : !$OMP PRIVATE(k,j,i,my_rho,my_rho_n,my_rho_n2) COLLAPSE(3)
1035 : DO k = bo(1, 3), bo(2, 3)
1036 : DO j = bo(1, 2), bo(2, 2)
1037 : DO i = bo(1, 1), bo(2, 1)
1038 : my_rho = rhoa(i, j, k) + rhob(i, j, k)
1039 : IF (my_rho < rho_smooth_cutoff) THEN
1040 : IF (my_rho < rho_cutoff) THEN
1041 : pot(i, j, k) = 0.0_dp
1042 : ELSEIF (my_rho < rho_smooth_cutoff_2) THEN
1043 : my_rho_n = (my_rho - rho_cutoff)/rho_smooth_cutoff_range_2
1044 : my_rho_n2 = my_rho_n*my_rho_n
1045 : pot(i, j, k) = pot(i, j, k)* &
1046 : my_rho_n2*(my_rho_n - 0.5_dp*my_rho_n2)
1047 : ELSE
1048 : my_rho_n = 2.0_dp - (my_rho - rho_cutoff)/rho_smooth_cutoff_range_2
1049 : my_rho_n2 = my_rho_n*my_rho_n
1050 : pot(i, j, k) = pot(i, j, k)* &
1051 : (1.0_dp - my_rho_n2*(my_rho_n - 0.5_dp*my_rho_n2))
1052 : END IF
1053 : END IF
1054 : END DO
1055 : END DO
1056 : END DO
1057 : END IF
1058 : END IF
1059 : END IF
1060 247514 : END SUBROUTINE smooth_cutoff
1061 :
1062 64 : SUBROUTINE calc_xc_density(pot, rho, rho_cutoff)
1063 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: pot
1064 : TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(INOUT) :: rho
1065 : REAL(kind=dp), INTENT(in) :: rho_cutoff
1066 :
1067 : INTEGER :: i, j, k, nspins
1068 : INTEGER, DIMENSION(2, 3) :: bo
1069 : REAL(kind=dp) :: eps1, eps2, my_rho, my_pot
1070 :
1071 256 : bo(1, :) = LBOUND(pot%array)
1072 256 : bo(2, :) = UBOUND(pot%array)
1073 64 : nspins = SIZE(rho)
1074 :
1075 64 : eps1 = rho_cutoff*1.E-4_dp
1076 64 : eps2 = rho_cutoff
1077 :
1078 3160 : DO k = bo(1, 3), bo(2, 3)
1079 161272 : DO j = bo(1, 2), bo(2, 2)
1080 4529376 : DO i = bo(1, 1), bo(2, 1)
1081 4368168 : my_pot = pot%array(i, j, k)
1082 4368168 : IF (nspins == 2) THEN
1083 339714 : my_rho = rho(1)%array(i, j, k) + rho(2)%array(i, j, k)
1084 : ELSE
1085 4028454 : my_rho = rho(1)%array(i, j, k)
1086 : END IF
1087 4526280 : IF (my_rho > eps1) THEN
1088 4139834 : pot%array(i, j, k) = my_pot/my_rho
1089 228334 : ELSE IF (my_rho < eps2) THEN
1090 228334 : pot%array(i, j, k) = 0.0_dp
1091 : ELSE
1092 0 : pot%array(i, j, k) = MIN(my_pot/my_rho, my_rho**(1._dp/3._dp))
1093 : END IF
1094 : END DO
1095 : END DO
1096 : END DO
1097 :
1098 64 : END SUBROUTINE calc_xc_density
1099 :
1100 : ! **************************************************************************************************
1101 : !> \brief Exchange and Correlation functional calculations
1102 : !> \param vxc_rho will contain the v_xc part that depend on rho
1103 : !> (if one of the chosen xc functionals has it it is allocated and you
1104 : !> are responsible for it)
1105 : !> \param vxc_tau will contain the kinetic tau part of v_xc
1106 : !> (if one of the chosen xc functionals has it it is allocated and you
1107 : !> are responsible for it)
1108 : !> \param exc the xc energy
1109 : !> \param rho_r the value of the density in the real space
1110 : !> \param rho_g value of the density in the g space (needs to be associated
1111 : !> only for gradient corrections)
1112 : !> \param tau value of the kinetic density tau on the grid (can be null,
1113 : !> used only with meta functionals)
1114 : !> \param xc_section which functional to calculate, and how to do it
1115 : !> \param pw_pool the pool for the grids
1116 : !> \param compute_virial ...
1117 : !> \param virial_xc ...
1118 : !> \param exc_r the value of the xc functional in the real space
1119 : !> \par History
1120 : !> JGH (13-Jun-2002): adaptation to new functionals
1121 : !> Fawzi (11.2002): drho_g(1:3)->drho_g
1122 : !> Fawzi (1.2003). lsd version
1123 : !> Fawzi (11.2003): version using the new xc interface
1124 : !> Fawzi (03.2004): fft free for smoothed density and derivs, gga lsd
1125 : !> Fawzi (04.2004): metafunctionals
1126 : !> mguidon (12.2008) : laplace functionals
1127 : !> \author fawzi; based LDA version of JGH, based on earlier version of apsi
1128 : !> \note
1129 : !> Beware: some really dirty pointer handling!
1130 : !> energy should be kept consistent with xc_exc_calc
1131 : ! **************************************************************************************************
1132 108860 : SUBROUTINE xc_vxc_pw_create(vxc_rho, vxc_tau, exc, rho_r, rho_g, tau, xc_section, &
1133 : pw_pool, compute_virial, virial_xc, exc_r)
1134 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: vxc_rho, vxc_tau
1135 : REAL(KIND=dp), INTENT(out) :: exc
1136 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, tau
1137 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
1138 : TYPE(section_vals_type), POINTER :: xc_section
1139 : TYPE(pw_pool_type), POINTER :: pw_pool
1140 : LOGICAL :: compute_virial
1141 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT) :: virial_xc
1142 : TYPE(pw_r3d_rs_type), INTENT(INOUT), OPTIONAL :: exc_r
1143 :
1144 : CHARACTER(len=*), PARAMETER :: routineN = 'xc_vxc_pw_create'
1145 : INTEGER, DIMENSION(2), PARAMETER :: norm_drho_spin_name = [deriv_norm_drhoa, deriv_norm_drhob]
1146 :
1147 : INTEGER :: handle, idir, ispin, jdir, &
1148 : npoints, nspins, &
1149 : xc_deriv_method_id, xc_rho_smooth_id, deriv_id
1150 : INTEGER, DIMENSION(2, 3) :: bo
1151 : LOGICAL :: dealloc_pw_to_deriv, has_laplace, &
1152 : has_tau, lsd, use_virial, has_gradient, &
1153 : has_derivs, has_rho, dealloc_pw_to_deriv_rho
1154 : REAL(KIND=dp) :: density_smooth_cut_range, drho_cutoff, &
1155 : rho_cutoff
1156 108860 : REAL(kind=dp), DIMENSION(:, :, :), POINTER :: deriv_data, norm_drho, norm_drho_spin, &
1157 217720 : rho, rhoa, rhob
1158 : TYPE(cp_sll_xc_deriv_type), POINTER :: pos
1159 : TYPE(pw_grid_type), POINTER :: pw_grid
1160 762020 : TYPE(pw_r3d_rs_type), DIMENSION(3) :: pw_to_deriv, pw_to_deriv_rho
1161 : TYPE(pw_c1d_gs_type) :: tmp_g, vxc_g
1162 : TYPE(pw_r3d_rs_type) :: v_drho_r, virial_pw
1163 : TYPE(xc_derivative_set_type) :: deriv_set
1164 : TYPE(xc_derivative_type), POINTER :: deriv_att
1165 : TYPE(xc_rho_set_type) :: rho_set
1166 :
1167 108860 : CALL timeset(routineN, handle)
1168 108860 : NULLIFY (norm_drho_spin, norm_drho, pos)
1169 :
1170 108860 : pw_grid => rho_r(1)%pw_grid
1171 :
1172 108860 : CPASSERT(ASSOCIATED(xc_section))
1173 108860 : CPASSERT(ASSOCIATED(pw_pool))
1174 108860 : CPASSERT(.NOT. ASSOCIATED(vxc_rho))
1175 108860 : CPASSERT(.NOT. ASSOCIATED(vxc_tau))
1176 108860 : nspins = SIZE(rho_r)
1177 108860 : lsd = (nspins /= 1)
1178 108860 : IF (lsd) THEN
1179 22337 : CPASSERT(nspins == 2)
1180 : END IF
1181 :
1182 108860 : use_virial = compute_virial
1183 108860 : virial_xc = 0.0_dp
1184 :
1185 1088600 : bo = rho_r(1)%pw_grid%bounds_local
1186 108860 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
1187 :
1188 : ! calculate the potential derivatives
1189 : CALL xc_rho_set_and_dset_create(rho_set=rho_set, deriv_set=deriv_set, &
1190 : deriv_order=1, rho_r=rho_r, rho_g=rho_g, tau=tau, &
1191 : xc_section=xc_section, &
1192 : pw_pool=pw_pool, &
1193 108860 : calc_potential=.TRUE.)
1194 :
1195 : CALL section_vals_val_get(xc_section, "XC_GRID%XC_DERIV", &
1196 108860 : i_val=xc_deriv_method_id)
1197 : CALL section_vals_val_get(xc_section, "XC_GRID%XC_SMOOTH_RHO", &
1198 108860 : i_val=xc_rho_smooth_id)
1199 : CALL section_vals_val_get(xc_section, "DENSITY_SMOOTH_CUTOFF_RANGE", &
1200 108860 : r_val=density_smooth_cut_range)
1201 :
1202 : CALL xc_rho_set_get(rho_set, rho_cutoff=rho_cutoff, &
1203 108860 : drho_cutoff=drho_cutoff)
1204 :
1205 108860 : CALL check_for_derivatives(deriv_set, lsd, has_rho, has_gradient, has_tau, has_laplace)
1206 : ! check for unknown derivatives
1207 108860 : has_derivs = has_rho .OR. has_gradient .OR. has_tau .OR. has_laplace
1208 :
1209 457777 : ALLOCATE (vxc_rho(nspins))
1210 :
1211 : CALL xc_rho_set_get(rho_set, rho=rho, rhoa=rhoa, rhob=rhob, &
1212 108860 : can_return_null=.TRUE.)
1213 :
1214 : ! recover the vxc arrays
1215 108860 : IF (lsd) THEN
1216 22337 : CALL xc_dset_recover_pw(deriv_set, [deriv_rhoa], vxc_rho(1), pw_grid, pw_pool)
1217 22337 : CALL xc_dset_recover_pw(deriv_set, [deriv_rhob], vxc_rho(2), pw_grid, pw_pool)
1218 :
1219 : ELSE
1220 86523 : CALL xc_dset_recover_pw(deriv_set, [deriv_rho], vxc_rho(1), pw_grid, pw_pool)
1221 : END IF
1222 :
1223 108860 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_norm_drho])
1224 108860 : IF (ASSOCIATED(deriv_att)) THEN
1225 57785 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
1226 :
1227 : CALL xc_rho_set_get(rho_set, norm_drho=norm_drho, &
1228 : rho_cutoff=rho_cutoff, &
1229 : drho_cutoff=drho_cutoff, &
1230 57785 : can_return_null=.TRUE.)
1231 57785 : CALL xc_rho_set_recover_pw(rho_set, pw_grid, pw_pool, dealloc_pw_to_deriv_rho, drho=pw_to_deriv_rho)
1232 :
1233 57785 : CPASSERT(ASSOCIATED(deriv_data))
1234 57785 : IF (use_virial) THEN
1235 1522 : CALL pw_pool%create_pw(virial_pw)
1236 1522 : CALL pw_zero(virial_pw)
1237 6088 : DO idir = 1, 3
1238 4566 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(virial_pw,pw_to_deriv_rho,deriv_data,idir)
1239 : virial_pw%array(:, :, :) = pw_to_deriv_rho(idir)%array(:, :, :)*deriv_data(:, :, :)
1240 : !$OMP END PARALLEL WORKSHARE
1241 15220 : DO jdir = 1, idir
1242 : virial_xc(idir, jdir) = -pw_grid%dvol* &
1243 : accurate_dot_product(virial_pw%array(:, :, :), &
1244 9132 : pw_to_deriv_rho(jdir)%array(:, :, :))
1245 13698 : virial_xc(jdir, idir) = virial_xc(idir, jdir)
1246 : END DO
1247 : END DO
1248 1522 : CALL pw_pool%give_back_pw(virial_pw)
1249 : END IF ! use_virial
1250 231140 : DO idir = 1, 3
1251 173355 : CPASSERT(ASSOCIATED(pw_to_deriv_rho(idir)%array))
1252 231140 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(deriv_data,pw_to_deriv_rho,idir)
1253 : pw_to_deriv_rho(idir)%array(:, :, :) = pw_to_deriv_rho(idir)%array(:, :, :)*deriv_data(:, :, :)
1254 : !$OMP END PARALLEL WORKSHARE
1255 : END DO
1256 :
1257 : ! Deallocate pw to save memory
1258 57785 : CALL pw_pool%give_back_cr3d(deriv_att%deriv_data)
1259 :
1260 : END IF
1261 :
1262 108860 : IF ((has_gradient .AND. xc_requires_tmp_g(xc_deriv_method_id)) .OR. pw_grid%spherical) THEN
1263 57059 : CALL pw_pool%create_pw(vxc_g)
1264 57059 : IF (.NOT. pw_grid%spherical) THEN
1265 57059 : CALL pw_pool%create_pw(tmp_g)
1266 : END IF
1267 : END IF
1268 :
1269 240057 : DO ispin = 1, nspins
1270 :
1271 131197 : IF (lsd) THEN
1272 44674 : IF (ispin == 1) THEN
1273 : CALL xc_rho_set_get(rho_set, norm_drhoa=norm_drho_spin, &
1274 22337 : can_return_null=.TRUE.)
1275 22337 : IF (ASSOCIATED(norm_drho_spin)) CALL xc_rho_set_recover_pw( &
1276 13848 : rho_set, pw_grid, pw_pool, dealloc_pw_to_deriv, drhoa=pw_to_deriv)
1277 : ELSE
1278 : CALL xc_rho_set_get(rho_set, norm_drhob=norm_drho_spin, &
1279 22337 : can_return_null=.TRUE.)
1280 22337 : IF (ASSOCIATED(norm_drho_spin)) CALL xc_rho_set_recover_pw( &
1281 13848 : rho_set, pw_grid, pw_pool, dealloc_pw_to_deriv, drhob=pw_to_deriv)
1282 : END IF
1283 :
1284 89348 : deriv_att => xc_dset_get_derivative(deriv_set, [norm_drho_spin_name(ispin)])
1285 44674 : IF (ASSOCIATED(deriv_att)) THEN
1286 : CPASSERT(lsd)
1287 27696 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
1288 :
1289 27696 : IF (use_virial) THEN
1290 96 : CALL pw_pool%create_pw(virial_pw)
1291 96 : CALL pw_zero(virial_pw)
1292 384 : DO idir = 1, 3
1293 288 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(deriv_data,pw_to_deriv,virial_pw,idir)
1294 : virial_pw%array(:, :, :) = pw_to_deriv(idir)%array(:, :, :)*deriv_data(:, :, :)
1295 : !$OMP END PARALLEL WORKSHARE
1296 960 : DO jdir = 1, idir
1297 : virial_xc(idir, jdir) = virial_xc(idir, jdir) - pw_grid%dvol* &
1298 : accurate_dot_product(virial_pw%array(:, :, :), &
1299 576 : pw_to_deriv(jdir)%array(:, :, :))
1300 864 : virial_xc(jdir, idir) = virial_xc(idir, jdir)
1301 : END DO
1302 : END DO
1303 96 : CALL pw_pool%give_back_pw(virial_pw)
1304 : END IF ! use_virial
1305 :
1306 110784 : DO idir = 1, 3
1307 110784 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(deriv_data,idir,pw_to_deriv)
1308 : pw_to_deriv(idir)%array(:, :, :) = deriv_data(:, :, :)* &
1309 : pw_to_deriv(idir)%array(:, :, :)
1310 : !$OMP END PARALLEL WORKSHARE
1311 : END DO
1312 : END IF ! deriv_att
1313 :
1314 : END IF ! LSD
1315 :
1316 131197 : IF (ASSOCIATED(pw_to_deriv_rho(1)%array)) THEN
1317 70939 : IF (.NOT. ASSOCIATED(pw_to_deriv(1)%array)) THEN
1318 44631 : pw_to_deriv = pw_to_deriv_rho
1319 44631 : dealloc_pw_to_deriv = ((.NOT. lsd) .OR. (ispin == 2))
1320 44631 : dealloc_pw_to_deriv = dealloc_pw_to_deriv .AND. dealloc_pw_to_deriv_rho
1321 : ELSE
1322 : ! This branch is called in case of open-shell systems
1323 : ! Add the contributions from norm_drho and norm_drho_spin
1324 105232 : DO idir = 1, 3
1325 78924 : CALL pw_axpy(pw_to_deriv_rho(idir), pw_to_deriv(idir))
1326 105232 : IF (ispin == 2) THEN
1327 39462 : IF (dealloc_pw_to_deriv_rho) THEN
1328 39462 : CALL pw_pool%give_back_pw(pw_to_deriv_rho(idir))
1329 : END IF
1330 : END IF
1331 : END DO
1332 : END IF
1333 : END IF
1334 :
1335 131197 : IF (ASSOCIATED(pw_to_deriv(1)%array)) THEN
1336 289308 : DO idir = 1, 3
1337 289308 : CALL pw_scale(pw_to_deriv(idir), -1.0_dp)
1338 : END DO
1339 :
1340 72327 : CALL xc_pw_divergence(xc_deriv_method_id, pw_to_deriv, tmp_g, vxc_g, vxc_rho(ispin))
1341 :
1342 72327 : IF (dealloc_pw_to_deriv) THEN
1343 289308 : DO idir = 1, 3
1344 289308 : CALL pw_pool%give_back_pw(pw_to_deriv(idir))
1345 : END DO
1346 : END IF
1347 : END IF
1348 :
1349 : ! ** Add laplace part to vxc_rho
1350 131197 : IF (has_laplace) THEN
1351 896 : IF (lsd) THEN
1352 352 : IF (ispin == 1) THEN
1353 : deriv_id = deriv_laplace_rhoa
1354 : ELSE
1355 176 : deriv_id = deriv_laplace_rhob
1356 : END IF
1357 : ELSE
1358 : deriv_id = deriv_laplace_rho
1359 : END IF
1360 :
1361 1792 : CALL xc_dset_recover_pw(deriv_set, [deriv_id], pw_to_deriv(1), pw_grid)
1362 :
1363 896 : IF (use_virial) CALL virial_laplace(rho_r(ispin), pw_pool, virial_xc, pw_to_deriv(1)%array)
1364 :
1365 896 : CALL xc_pw_laplace(pw_to_deriv(1), pw_pool, xc_deriv_method_id)
1366 :
1367 896 : CALL pw_axpy(pw_to_deriv(1), vxc_rho(ispin))
1368 :
1369 896 : CALL pw_pool%give_back_pw(pw_to_deriv(1))
1370 : END IF
1371 :
1372 131197 : IF (pw_grid%spherical) THEN
1373 : ! filter vxc
1374 0 : CALL pw_transfer(vxc_rho(ispin), vxc_g)
1375 0 : CALL pw_transfer(vxc_g, vxc_rho(ispin))
1376 : END IF
1377 : CALL smooth_cutoff(pot=vxc_rho(ispin)%array, rho=rho, rhoa=rhoa, rhob=rhob, &
1378 : rho_cutoff=rho_cutoff*density_smooth_cut_range, &
1379 131197 : rho_smooth_cutoff_range=density_smooth_cut_range)
1380 :
1381 131197 : v_drho_r = vxc_rho(ispin)
1382 131197 : CALL pw_pool%create_pw(vxc_rho(ispin))
1383 131197 : CALL xc_pw_smooth(v_drho_r, vxc_rho(ispin), xc_rho_smooth_id)
1384 240057 : CALL pw_pool%give_back_pw(v_drho_r)
1385 : END DO
1386 :
1387 108860 : CALL pw_pool%give_back_pw(vxc_g)
1388 108860 : CALL pw_pool%give_back_pw(tmp_g)
1389 :
1390 : ! 0-deriv -> value of exc
1391 : ! this has to be kept consistent with xc_exc_calc
1392 108860 : IF (has_derivs) THEN
1393 108540 : CALL xc_dset_recover_pw(deriv_set, [INTEGER::], v_drho_r, pw_grid)
1394 :
1395 : CALL smooth_cutoff(pot=v_drho_r%array, rho=rho, rhoa=rhoa, rhob=rhob, &
1396 : rho_cutoff=rho_cutoff, &
1397 108540 : rho_smooth_cutoff_range=density_smooth_cut_range)
1398 :
1399 108540 : exc = pw_integrate_function(v_drho_r)
1400 : !
1401 : ! return the xc functional value at the grid points
1402 : !
1403 108540 : IF (PRESENT(exc_r)) THEN
1404 96 : exc_r = v_drho_r
1405 : ELSE
1406 108444 : CALL v_drho_r%release()
1407 : END IF
1408 : ELSE
1409 320 : exc = 0.0_dp
1410 : END IF
1411 :
1412 108860 : CALL xc_rho_set_release(rho_set, pw_pool=pw_pool)
1413 :
1414 : ! tau part
1415 108860 : IF (has_tau) THEN
1416 10158 : ALLOCATE (vxc_tau(nspins))
1417 2396 : IF (lsd) THEN
1418 574 : CALL xc_dset_recover_pw(deriv_set, [deriv_tau_a], vxc_tau(1), pw_grid)
1419 574 : CALL xc_dset_recover_pw(deriv_set, [deriv_tau_b], vxc_tau(2), pw_grid)
1420 : ELSE
1421 1822 : CALL xc_dset_recover_pw(deriv_set, [deriv_tau], vxc_tau(1), pw_grid)
1422 : END IF
1423 5366 : DO ispin = 1, nspins
1424 5366 : CPASSERT(ASSOCIATED(vxc_tau(ispin)%array))
1425 : END DO
1426 : END IF
1427 108860 : CALL xc_dset_release(deriv_set)
1428 :
1429 108860 : CALL timestop(handle)
1430 :
1431 2286060 : END SUBROUTINE xc_vxc_pw_create
1432 :
1433 : ! **************************************************************************************************
1434 : !> \brief calculates just the exchange and correlation energy
1435 : !> (no vxc)
1436 : !> \param rho_r realspace density on the grid
1437 : !> \param rho_g g-space density on the grid
1438 : !> \param tau kinetic energy density on the grid
1439 : !> \param xc_section XC parameters
1440 : !> \param pw_pool pool of plain-wave grids
1441 : !> \return the XC energy
1442 : !> \par History
1443 : !> 11.2003 created [fawzi]
1444 : !> \author fawzi
1445 : !> \note
1446 : !> has to be kept consistent with xc_vxc_pw_create
1447 : ! **************************************************************************************************
1448 15550 : FUNCTION xc_exc_calc(rho_r, rho_g, tau, xc_section, pw_pool) &
1449 : RESULT(exc)
1450 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, tau
1451 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
1452 : TYPE(section_vals_type), POINTER :: xc_section
1453 : TYPE(pw_pool_type), POINTER :: pw_pool
1454 : REAL(kind=dp) :: exc
1455 :
1456 : CHARACTER(len=*), PARAMETER :: routineN = 'xc_exc_calc'
1457 :
1458 : INTEGER :: handle
1459 : REAL(dp) :: density_smooth_cut_range, rho_cutoff
1460 7775 : REAL(dp), DIMENSION(:, :, :), POINTER :: e_0
1461 : TYPE(xc_derivative_set_type) :: deriv_set
1462 : TYPE(xc_derivative_type), POINTER :: deriv
1463 : TYPE(xc_rho_set_type) :: rho_set
1464 :
1465 7775 : CALL timeset(routineN, handle)
1466 7775 : NULLIFY (deriv, e_0)
1467 7775 : exc = 0.0_dp
1468 :
1469 : ! this has to be consistent with what is done in xc_vxc_pw_create
1470 : CALL xc_rho_set_and_dset_create(rho_set=rho_set, &
1471 : deriv_set=deriv_set, deriv_order=0, &
1472 : rho_r=rho_r, rho_g=rho_g, tau=tau, xc_section=xc_section, &
1473 : pw_pool=pw_pool, &
1474 7775 : calc_potential=.FALSE.)
1475 7775 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::])
1476 :
1477 7775 : IF (ASSOCIATED(deriv)) THEN
1478 7775 : CALL xc_derivative_get(deriv, deriv_data=e_0)
1479 :
1480 : CALL section_vals_val_get(xc_section, "DENSITY_CUTOFF", &
1481 7775 : r_val=rho_cutoff)
1482 : CALL section_vals_val_get(xc_section, "DENSITY_SMOOTH_CUTOFF_RANGE", &
1483 7775 : r_val=density_smooth_cut_range)
1484 : CALL smooth_cutoff(pot=e_0, rho=rho_set%rho, &
1485 : rhoa=rho_set%rhoa, rhob=rho_set%rhob, &
1486 : rho_cutoff=rho_cutoff, &
1487 7775 : rho_smooth_cutoff_range=density_smooth_cut_range)
1488 :
1489 7775 : exc = accurate_sum(e_0)*rho_r(1)%pw_grid%dvol
1490 7775 : IF (rho_r(1)%pw_grid%para%mode == PW_MODE_DISTRIBUTED) THEN
1491 7650 : CALL rho_r(1)%pw_grid%para%group%sum(exc)
1492 : END IF
1493 :
1494 7775 : CALL xc_rho_set_release(rho_set, pw_pool=pw_pool)
1495 7775 : CALL xc_dset_release(deriv_set)
1496 : END IF
1497 7775 : CALL timestop(handle)
1498 147725 : END FUNCTION xc_exc_calc
1499 :
1500 : ! **************************************************************************************************
1501 : !> \brief Caller routine to calculate the second order potential in the direction of rho1_r
1502 : !> \param v_xc XC potential, will be allocated, to be integrated with the KS density
1503 : !> \param v_xc_tau ...
1504 : !> \param deriv_set XC derivatives from xc_prep_2nd_deriv
1505 : !> \param rho_set XC rho set from KS rho from xc_prep_2nd_deriv
1506 : !> \param rho1_r first-order density in r space
1507 : !> \param rho1_g first-order density in g space
1508 : !> \param tau1_r ...
1509 : !> \param pw_pool pw pool to create new grids
1510 : !> \param xc_section XC section to calculate the derivatives from
1511 : !> \param gapw whether to carry out GAPW (not possible with numerical derivatives)
1512 : !> \param vxg GAPW potential
1513 : !> \param lsd_singlets ...
1514 : !> \param do_excitations ...
1515 : !> \param do_triplet ...
1516 : !> \param do_tddft ...
1517 : !> \param compute_virial ...
1518 : !> \param virial_xc virial terms will be collected here
1519 : ! **************************************************************************************************
1520 13844 : SUBROUTINE xc_calc_2nd_deriv(v_xc, v_xc_tau, deriv_set, rho_set, rho1_r, rho1_g, tau1_r, &
1521 : pw_pool, xc_section, gapw, vxg, &
1522 : lsd_singlets, do_excitations, do_triplet, do_tddft, &
1523 : compute_virial, virial_xc)
1524 :
1525 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_xc, v_xc_tau
1526 : TYPE(xc_derivative_set_type) :: deriv_set
1527 : TYPE(xc_rho_set_type) :: rho_set
1528 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho1_r, tau1_r
1529 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho1_g
1530 : TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool
1531 : TYPE(section_vals_type), INTENT(IN), POINTER :: xc_section
1532 : LOGICAL, INTENT(IN) :: gapw
1533 : REAL(KIND=dp), DIMENSION(:, :, :, :), OPTIONAL, &
1534 : POINTER :: vxg
1535 : LOGICAL, INTENT(IN), OPTIONAL :: lsd_singlets, do_excitations, &
1536 : do_triplet, do_tddft, compute_virial
1537 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT), &
1538 : OPTIONAL :: virial_xc
1539 :
1540 : CHARACTER(len=*), PARAMETER :: routineN = 'xc_calc_2nd_deriv'
1541 :
1542 : INTEGER :: handle, ispin, nspins
1543 : INTEGER, DIMENSION(2, 3) :: bo
1544 : LOGICAL :: lsd, my_compute_virial, &
1545 : my_do_excitations, my_do_tddft, &
1546 : my_do_triplet, my_lsd_singlets
1547 : REAL(KIND=dp) :: fac
1548 : TYPE(section_vals_type), POINTER :: xc_fun_section
1549 : TYPE(xc_rho_cflags_type) :: needs
1550 : TYPE(xc_rho_set_type) :: rho1_set
1551 :
1552 13844 : CALL timeset(routineN, handle)
1553 :
1554 13844 : my_compute_virial = .FALSE.
1555 13844 : IF (PRESENT(compute_virial)) my_compute_virial = compute_virial
1556 :
1557 13844 : my_do_tddft = .FALSE.
1558 13844 : IF (PRESENT(do_tddft)) my_do_tddft = do_tddft
1559 :
1560 13844 : my_do_excitations = .FALSE.
1561 13844 : IF (PRESENT(do_excitations)) my_do_excitations = do_excitations
1562 :
1563 13844 : my_lsd_singlets = .FALSE.
1564 13844 : IF (PRESENT(lsd_singlets)) my_lsd_singlets = lsd_singlets
1565 :
1566 13844 : my_do_triplet = .FALSE.
1567 13844 : IF (PRESENT(do_triplet)) my_do_triplet = do_triplet
1568 :
1569 13844 : nspins = SIZE(rho1_r)
1570 13844 : lsd = (nspins == 2)
1571 13844 : IF (nspins == 1 .AND. my_do_excitations .AND. (my_lsd_singlets .OR. my_do_triplet)) THEN
1572 0 : nspins = 2
1573 0 : lsd = .TRUE.
1574 : END IF
1575 :
1576 13844 : NULLIFY (v_xc, v_xc_tau)
1577 57318 : ALLOCATE (v_xc(nspins))
1578 29630 : DO ispin = 1, nspins
1579 15786 : CALL pw_pool%create_pw(v_xc(ispin))
1580 29630 : CALL pw_zero(v_xc(ispin))
1581 : END DO
1582 :
1583 13844 : xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
1584 13844 : needs = xc_functionals_get_needs(xc_fun_section, lsd, .TRUE.)
1585 :
1586 13844 : IF (needs%tau .OR. needs%tau_spin) THEN
1587 514 : IF (.NOT. ASSOCIATED(tau1_r)) &
1588 0 : CPABORT("Tau-dependent functionals requires allocated kinetic energy density grid")
1589 2172 : ALLOCATE (v_xc_tau(nspins))
1590 1144 : DO ispin = 1, nspins
1591 630 : CALL pw_pool%create_pw(v_xc_tau(ispin))
1592 14474 : CALL pw_zero(v_xc_tau(ispin))
1593 : END DO
1594 : END IF
1595 :
1596 13844 : IF (section_get_lval(xc_section, "2ND_DERIV_ANALYTICAL") .AND. .NOT. my_do_tddft) THEN
1597 : !------!
1598 : ! rho1 !
1599 : !------!
1600 130640 : bo = rho1_r(1)%pw_grid%bounds_local
1601 : ! create the place where to store the argument for the functionals
1602 : CALL xc_rho_set_create(rho1_set, bo, &
1603 : rho_cutoff=section_get_rval(xc_section, "DENSITY_CUTOFF"), &
1604 : drho_cutoff=section_get_rval(xc_section, "GRADIENT_CUTOFF"), &
1605 13064 : tau_cutoff=section_get_rval(xc_section, "TAU_CUTOFF"))
1606 :
1607 : ! calculate the arguments needed by the functionals
1608 : CALL xc_rho_set_update(rho1_set, rho1_r, rho1_g, tau1_r, needs, &
1609 : section_get_ival(xc_section, "XC_GRID%XC_DERIV"), &
1610 : section_get_ival(xc_section, "XC_GRID%XC_SMOOTH_RHO"), &
1611 13064 : pw_pool)
1612 :
1613 13064 : fac = 0._dp
1614 13064 : IF (nspins == 1 .AND. my_do_excitations) THEN
1615 1064 : IF (my_lsd_singlets) fac = 1.0_dp
1616 1064 : IF (my_do_triplet) fac = -1.0_dp
1617 : END IF
1618 :
1619 : CALL xc_calc_2nd_deriv_analytical(v_xc, v_xc_tau, deriv_set, rho_set, &
1620 : rho1_set, pw_pool, xc_section, gapw, vxg=vxg, &
1621 26128 : tddfpt_fac=fac, compute_virial=compute_virial, virial_xc=virial_xc)
1622 :
1623 13064 : CALL xc_rho_set_release(rho1_set)
1624 :
1625 : ELSE
1626 780 : IF (gapw) CPABORT("Numerical 2nd derivatives not implemented with GAPW")
1627 :
1628 : CALL xc_calc_2nd_deriv_numerical(v_xc, v_xc_tau, rho_set, rho1_r, rho1_g, tau1_r, &
1629 : pw_pool, xc_section, &
1630 : my_do_excitations .AND. my_do_triplet, &
1631 780 : compute_virial, virial_xc, deriv_set)
1632 : END IF
1633 :
1634 13844 : CALL timestop(handle)
1635 :
1636 304568 : END SUBROUTINE xc_calc_2nd_deriv
1637 :
1638 : ! **************************************************************************************************
1639 : !> \brief calculates 2nd derivative numerically
1640 : !> \param v_xc potential to be calculated (has to be allocated already)
1641 : !> \param v_tau tau-part of the potential to be calculated (has to be allocated already)
1642 : !> \param rho_set KS density from xc_prep_2nd_deriv
1643 : !> \param rho1_r first-order density in r-space
1644 : !> \param rho1_g first-order density in g-space
1645 : !> \param tau1_r first-order kinetic-energy density in r-space
1646 : !> \param pw_pool pw pool for new grids
1647 : !> \param xc_section XC section to calculate the derivatives from
1648 : !> \param do_triplet ...
1649 : !> \param calc_virial whether to calculate virial terms
1650 : !> \param virial_xc collects stress tensor components (no metaGGAs!)
1651 : !> \param deriv_set deriv set from xc_prep_2nd_deriv (only for virials)
1652 : ! **************************************************************************************************
1653 832 : SUBROUTINE xc_calc_2nd_deriv_numerical(v_xc, v_tau, rho_set, rho1_r, rho1_g, tau1_r, &
1654 : pw_pool, xc_section, &
1655 : do_triplet, calc_virial, virial_xc, deriv_set)
1656 :
1657 : TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(IN), POINTER :: v_xc, v_tau
1658 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
1659 : TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(IN), POINTER :: rho1_r, tau1_r
1660 : TYPE(pw_c1d_gs_type), DIMENSION(:), INTENT(IN), POINTER :: rho1_g
1661 : TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool
1662 : TYPE(section_vals_type), INTENT(IN), POINTER :: xc_section
1663 : LOGICAL, INTENT(IN) :: do_triplet
1664 : LOGICAL, INTENT(IN), OPTIONAL :: calc_virial
1665 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT), &
1666 : OPTIONAL :: virial_xc
1667 : TYPE(xc_derivative_set_type), OPTIONAL :: deriv_set
1668 :
1669 : CHARACTER(len=*), PARAMETER :: routineN = 'xc_calc_2nd_deriv_numerical'
1670 : REAL(KIND=dp), DIMENSION(-4:4, 4), PARAMETER :: &
1671 : 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, &
1672 : 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, &
1673 : 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, &
1674 : 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])
1675 :
1676 : INTEGER :: handle, idir, ispin, nspins, istep, nsteps
1677 : INTEGER, DIMENSION(2, 3) :: bo
1678 : LOGICAL :: gradient_f, lsd, my_calc_virial, tau_f, laplace_f, rho_f
1679 : REAL(KIND=dp) :: exc, gradient_cut, h, weight, step, rho_cutoff
1680 832 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: dr1dr, dra1dra, drb1drb
1681 : REAL(KIND=dp), DIMENSION(3, 3) :: virial_dummy
1682 832 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: norm_drho, norm_drho2, norm_drho2a, &
1683 832 : norm_drho2b, norm_drhoa, norm_drhob, &
1684 2496 : rho, rho1, rho1a, rho1b, rhoa, rhob, &
1685 1664 : tau_a, tau_b, tau, tau1, tau1a, tau1b, laplace, laplace1, &
1686 832 : laplacea, laplaceb, laplace1a, laplace1b, &
1687 1664 : laplace2, laplace2a, laplace2b, deriv_data
1688 19968 : TYPE(cp_3d_r_cp_type), DIMENSION(3) :: drho, drho1, drho1a, drho1b, drhoa, drhob
1689 : TYPE(pw_r3d_rs_type) :: v_drho, v_drhoa, v_drhob
1690 832 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: vxc_rho, vxc_tau
1691 832 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
1692 832 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, tau_r
1693 : TYPE(pw_r3d_rs_type) :: virial_pw, v_laplace, v_laplacea, v_laplaceb
1694 : TYPE(section_vals_type), POINTER :: xc_fun_section
1695 : TYPE(xc_derivative_set_type) :: deriv_set1
1696 : TYPE(xc_rho_cflags_type) :: needs
1697 : TYPE(xc_rho_set_type) :: rho1_set, rho2_set
1698 :
1699 832 : CALL timeset(routineN, handle)
1700 :
1701 832 : my_calc_virial = .FALSE.
1702 832 : IF (PRESENT(calc_virial) .AND. PRESENT(virial_xc)) my_calc_virial = calc_virial
1703 :
1704 832 : nspins = SIZE(v_xc)
1705 :
1706 832 : NULLIFY (tau, tau_r, tau_a, tau_b)
1707 :
1708 832 : h = section_get_rval(xc_section, "STEP_SIZE")
1709 832 : nsteps = section_get_ival(xc_section, "NSTEPS")
1710 832 : IF (nsteps < LBOUND(weights, 2) .OR. nspins > UBOUND(weights, 2)) THEN
1711 0 : CPABORT("The number of steps must be a value from 1 to 4.")
1712 : END IF
1713 :
1714 832 : IF (nspins == 2) THEN
1715 324 : NULLIFY (vxc_rho, rho_g, vxc_tau)
1716 972 : ALLOCATE (rho_r(2))
1717 972 : DO ispin = 1, nspins
1718 972 : CALL pw_pool%create_pw(rho_r(ispin))
1719 : END DO
1720 324 : IF (ASSOCIATED(tau1_r) .AND. ASSOCIATED(v_tau)) THEN
1721 162 : ALLOCATE (tau_r(2))
1722 162 : DO ispin = 1, nspins
1723 162 : CALL pw_pool%create_pw(tau_r(ispin))
1724 : END DO
1725 : END IF
1726 324 : CALL xc_rho_set_get(rho_set, can_return_null=.TRUE., rhoa=rhoa, rhob=rhob, tau_a=tau_a, tau_b=tau_b)
1727 2496 : DO istep = -nsteps, nsteps
1728 2172 : IF (istep == 0) CYCLE
1729 1848 : weight = weights(istep, nsteps)/h
1730 1848 : step = REAL(istep, dp)*h
1731 : CALL calc_resp_potential_numer_ab(rho_r, rho_g, rho1_r, rhoa, rhob, vxc_rho, &
1732 1848 : tau_r, tau1_r, tau_a, tau_b, vxc_tau, xc_section, pw_pool, step)
1733 5544 : DO ispin = 1, nspins
1734 3696 : CALL pw_axpy(vxc_rho(ispin), v_xc(ispin), weight)
1735 5544 : IF (ASSOCIATED(vxc_tau) .AND. ASSOCIATED(v_tau)) THEN
1736 456 : CALL pw_axpy(vxc_tau(ispin), v_tau(ispin), weight)
1737 : END IF
1738 : END DO
1739 5544 : DO ispin = 1, nspins
1740 5544 : CALL vxc_rho(ispin)%release()
1741 : END DO
1742 1848 : DEALLOCATE (vxc_rho)
1743 2172 : IF (ASSOCIATED(vxc_tau)) THEN
1744 684 : DO ispin = 1, nspins
1745 684 : CALL vxc_tau(ispin)%release()
1746 : END DO
1747 228 : DEALLOCATE (vxc_tau)
1748 : END IF
1749 : END DO
1750 508 : ELSE IF (nspins == 1 .AND. do_triplet) THEN
1751 24 : NULLIFY (vxc_rho, vxc_tau, rho_g)
1752 72 : ALLOCATE (rho_r(2))
1753 72 : DO ispin = 1, 2
1754 72 : CALL pw_pool%create_pw(rho_r(ispin))
1755 : END DO
1756 24 : IF (ASSOCIATED(tau1_r) .AND. ASSOCIATED(v_tau)) THEN
1757 0 : ALLOCATE (tau_r(2))
1758 0 : DO ispin = 1, nspins
1759 0 : CALL pw_pool%create_pw(tau_r(ispin))
1760 : END DO
1761 : END IF
1762 24 : CALL xc_rho_set_get(rho_set, can_return_null=.TRUE., rhoa=rhoa, rhob=rhob, tau_a=tau_a, tau_b=tau_b)
1763 192 : DO istep = -nsteps, nsteps
1764 168 : IF (istep == 0) CYCLE
1765 144 : weight = weights(istep, nsteps)/h
1766 144 : step = REAL(istep, dp)*h
1767 : ! K(alpha,alpha)
1768 144 : !$OMP PARALLEL DEFAULT(NONE) SHARED(rho_r,rhoa,rhob,step,rho1_r,tau_r,tau_a,tau_b,tau1_r)
1769 : !$OMP WORKSHARE
1770 : rho_r(1)%array(:, :, :) = rhoa(:, :, :) + step*rho1_r(1)%array(:, :, :)
1771 : !$OMP END WORKSHARE NOWAIT
1772 : !$OMP WORKSHARE
1773 : rho_r(2)%array(:, :, :) = rhob(:, :, :)
1774 : !$OMP END WORKSHARE NOWAIT
1775 : IF (ASSOCIATED(tau1_r)) THEN
1776 : !$OMP WORKSHARE
1777 : tau_r(1)%array(:, :, :) = tau_a(:, :, :) + step*tau1_r(1)%array(:, :, :)
1778 : !$OMP END WORKSHARE NOWAIT
1779 : !$OMP WORKSHARE
1780 : tau_r(2)%array(:, :, :) = tau_b(:, :, :)
1781 : !$OMP END WORKSHARE NOWAIT
1782 : END IF
1783 : !$OMP END PARALLEL
1784 : CALL xc_vxc_pw_create(vxc_rho, vxc_tau, exc, rho_r, rho_g, tau_r, xc_section, &
1785 144 : pw_pool, .FALSE., virial_dummy)
1786 144 : CALL pw_axpy(vxc_rho(1), v_xc(1), weight)
1787 144 : IF (ASSOCIATED(vxc_tau) .AND. ASSOCIATED(v_tau)) THEN
1788 0 : CALL pw_axpy(vxc_tau(1), v_tau(1), weight)
1789 : END IF
1790 432 : DO ispin = 1, 2
1791 432 : CALL vxc_rho(ispin)%release()
1792 : END DO
1793 144 : DEALLOCATE (vxc_rho)
1794 144 : IF (ASSOCIATED(vxc_tau)) THEN
1795 0 : DO ispin = 1, 2
1796 0 : CALL vxc_tau(ispin)%release()
1797 : END DO
1798 0 : DEALLOCATE (vxc_tau)
1799 : END IF
1800 144 : !$OMP PARALLEL DEFAULT(NONE) SHARED(rho_r,rhoa,rhob,step,rho1_r,tau_r,tau_a,tau_b,tau1_r)
1801 : !$OMP WORKSHARE
1802 : ! K(alpha,beta)
1803 : rho_r(1)%array(:, :, :) = rhoa(:, :, :)
1804 : !$OMP END WORKSHARE NOWAIT
1805 : !$OMP WORKSHARE
1806 : rho_r(2)%array(:, :, :) = rhob(:, :, :) + step*rho1_r(1)%array(:, :, :)
1807 : !$OMP END WORKSHARE NOWAIT
1808 : IF (ASSOCIATED(tau1_r)) THEN
1809 : !$OMP WORKSHARE
1810 : tau_r(1)%array(:, :, :) = tau_a(:, :, :)
1811 : !$OMP END WORKSHARE NOWAIT
1812 : !$OMP WORKSHARE
1813 : tau_r(2)%array(:, :, :) = tau_b(:, :, :) + step*tau1_r(1)%array(:, :, :)
1814 : !$OMP END WORKSHARE NOWAIT
1815 : END IF
1816 : !$OMP END PARALLEL
1817 : CALL xc_vxc_pw_create(vxc_rho, vxc_tau, exc, rho_r, rho_g, tau_r, xc_section, &
1818 144 : pw_pool, .FALSE., virial_dummy)
1819 144 : CALL pw_axpy(vxc_rho(1), v_xc(1), weight)
1820 144 : IF (ASSOCIATED(vxc_tau) .AND. ASSOCIATED(v_tau)) THEN
1821 0 : CALL pw_axpy(vxc_tau(1), v_tau(1), weight)
1822 : END IF
1823 432 : DO ispin = 1, 2
1824 432 : CALL vxc_rho(ispin)%release()
1825 : END DO
1826 144 : DEALLOCATE (vxc_rho)
1827 168 : IF (ASSOCIATED(vxc_tau)) THEN
1828 0 : DO ispin = 1, 2
1829 0 : CALL vxc_tau(ispin)%release()
1830 : END DO
1831 0 : DEALLOCATE (vxc_tau)
1832 : END IF
1833 : END DO
1834 : ELSE
1835 484 : NULLIFY (vxc_rho, rho_r, rho_g, vxc_tau, tau_r, tau)
1836 968 : ALLOCATE (rho_r(1))
1837 484 : CALL pw_pool%create_pw(rho_r(1))
1838 484 : IF (ASSOCIATED(tau1_r) .AND. ASSOCIATED(v_tau)) THEN
1839 92 : ALLOCATE (tau_r(1))
1840 46 : CALL pw_pool%create_pw(tau_r(1))
1841 : END IF
1842 484 : CALL xc_rho_set_get(rho_set, can_return_null=.TRUE., rho=rho, tau=tau)
1843 3776 : DO istep = -nsteps, nsteps
1844 3292 : IF (istep == 0) CYCLE
1845 2808 : weight = weights(istep, nsteps)/h
1846 2808 : step = REAL(istep, dp)*h
1847 2808 : !$OMP PARALLEL DEFAULT(NONE) SHARED(rho_r,rho,step,rho1_r,tau1_r,tau,tau_r)
1848 : !$OMP WORKSHARE
1849 : rho_r(1)%array(:, :, :) = rho(:, :, :) + step*rho1_r(1)%array(:, :, :)
1850 : !$OMP END WORKSHARE NOWAIT
1851 : IF (ASSOCIATED(tau1_r) .AND. ASSOCIATED(tau) .AND. ASSOCIATED(tau1_r)) THEN
1852 : !$OMP WORKSHARE
1853 : tau_r(1)%array(:, :, :) = tau(:, :, :) + step*tau1_r(1)%array(:, :, :)
1854 : !$OMP END WORKSHARE NOWAIT
1855 : END IF
1856 : !$OMP END PARALLEL
1857 : CALL xc_vxc_pw_create(vxc_rho, vxc_tau, exc, rho_r, rho_g, tau_r, xc_section, &
1858 2808 : pw_pool, .FALSE., virial_dummy)
1859 2808 : CALL pw_axpy(vxc_rho(1), v_xc(1), weight)
1860 2808 : IF (ASSOCIATED(vxc_tau) .AND. ASSOCIATED(v_tau)) THEN
1861 276 : CALL pw_axpy(vxc_tau(1), v_tau(1), weight)
1862 : END IF
1863 2808 : CALL vxc_rho(1)%release()
1864 2808 : DEALLOCATE (vxc_rho)
1865 3292 : IF (ASSOCIATED(vxc_tau)) THEN
1866 276 : CALL vxc_tau(1)%release()
1867 276 : DEALLOCATE (vxc_tau)
1868 : END IF
1869 : END DO
1870 : END IF
1871 :
1872 832 : IF (my_calc_virial) THEN
1873 36 : lsd = (nspins == 2)
1874 36 : IF (nspins == 1 .AND. do_triplet) THEN
1875 0 : lsd = .TRUE.
1876 : END IF
1877 :
1878 36 : CALL check_for_derivatives(deriv_set, (nspins == 2), rho_f, gradient_f, tau_f, laplace_f)
1879 :
1880 : ! Calculate the virial terms
1881 : ! Those arising from the first derivatives are treated like in xc_calc_2nd_deriv_analytical
1882 : ! Those arising from the second derivatives are calculated numerically
1883 : ! We assume that all metaGGA functionals require the gradient
1884 36 : IF (gradient_f) THEN
1885 360 : bo = rho_set%local_bounds
1886 :
1887 : ! Create the work grid for the virial terms
1888 36 : CALL allocate_pw(virial_pw, pw_pool, bo)
1889 :
1890 36 : gradient_cut = section_get_rval(xc_section, "GRADIENT_CUTOFF")
1891 :
1892 : ! create the container to store the argument of the functionals
1893 : CALL xc_rho_set_create(rho1_set, bo, &
1894 : rho_cutoff=section_get_rval(xc_section, "DENSITY_CUTOFF"), &
1895 : drho_cutoff=gradient_cut, &
1896 36 : tau_cutoff=section_get_rval(xc_section, "TAU_CUTOFF"))
1897 :
1898 36 : xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
1899 36 : needs = xc_functionals_get_needs(xc_fun_section, lsd, .TRUE.)
1900 :
1901 : ! calculate the arguments needed by the functionals
1902 : CALL xc_rho_set_update(rho1_set, rho1_r, rho1_g, tau1_r, needs, &
1903 : section_get_ival(xc_section, "XC_GRID%XC_DERIV"), &
1904 : section_get_ival(xc_section, "XC_GRID%XC_SMOOTH_RHO"), &
1905 36 : pw_pool)
1906 :
1907 36 : IF (lsd) THEN
1908 : CALL xc_rho_set_get(rho_set, drhoa=drhoa, drhob=drhob, norm_drho=norm_drho, &
1909 : norm_drhoa=norm_drhoa, norm_drhob=norm_drhob, tau_a=tau_a, tau_b=tau_b, &
1910 10 : laplace_rhoa=laplacea, laplace_rhob=laplaceb, can_return_null=.TRUE.)
1911 : CALL xc_rho_set_get(rho1_set, rhoa=rho1a, rhob=rho1b, drhoa=drho1a, drhob=drho1b, laplace_rhoa=laplace1a, &
1912 10 : laplace_rhob=laplace1b, can_return_null=.TRUE.)
1913 :
1914 10 : CALL calc_drho_from_ab(drho, drhoa, drhob)
1915 10 : CALL calc_drho_from_ab(drho1, drho1a, drho1b)
1916 : ELSE
1917 26 : CALL xc_rho_set_get(rho_set, drho=drho, norm_drho=norm_drho, tau=tau, laplace_rho=laplace, can_return_null=.TRUE.)
1918 26 : CALL xc_rho_set_get(rho1_set, rho=rho1, drho=drho1, laplace_rho=laplace1, can_return_null=.TRUE.)
1919 : END IF
1920 :
1921 36 : CALL prepare_dr1dr(dr1dr, drho, drho1)
1922 :
1923 36 : IF (lsd) THEN
1924 10 : CALL prepare_dr1dr(dra1dra, drhoa, drho1a)
1925 10 : CALL prepare_dr1dr(drb1drb, drhob, drho1b)
1926 :
1927 10 : CALL allocate_pw(v_drho, pw_pool, bo)
1928 10 : CALL allocate_pw(v_drhoa, pw_pool, bo)
1929 10 : CALL allocate_pw(v_drhob, pw_pool, bo)
1930 :
1931 10 : IF (ASSOCIATED(norm_drhoa)) CALL apply_drho(deriv_set, [deriv_norm_drhoa], virial_pw, drhoa, drho1a, virial_xc, &
1932 10 : norm_drhoa, gradient_cut, dra1dra, v_drhoa%array)
1933 10 : IF (ASSOCIATED(norm_drhob)) CALL apply_drho(deriv_set, [deriv_norm_drhob], virial_pw, drhob, drho1b, virial_xc, &
1934 10 : norm_drhob, gradient_cut, drb1drb, v_drhob%array)
1935 10 : IF (ASSOCIATED(norm_drho)) CALL apply_drho(deriv_set, [deriv_norm_drho], virial_pw, drho, drho1, virial_xc, &
1936 6 : norm_drho, gradient_cut, dr1dr, v_drho%array)
1937 10 : IF (laplace_f) THEN
1938 2 : CALL xc_derivative_get(xc_dset_get_derivative(deriv_set, [deriv_laplace_rhoa]), deriv_data=deriv_data)
1939 2 : CPASSERT(ASSOCIATED(deriv_data))
1940 15026 : virial_pw%array(:, :, :) = -rho1a(:, :, :)
1941 2 : CALL virial_laplace(virial_pw, pw_pool, virial_xc, deriv_data)
1942 :
1943 2 : CALL allocate_pw(v_laplacea, pw_pool, bo)
1944 :
1945 2 : CALL xc_derivative_get(xc_dset_get_derivative(deriv_set, [deriv_laplace_rhob]), deriv_data=deriv_data)
1946 2 : CPASSERT(ASSOCIATED(deriv_data))
1947 15026 : virial_pw%array(:, :, :) = -rho1b(:, :, :)
1948 2 : CALL virial_laplace(virial_pw, pw_pool, virial_xc, deriv_data)
1949 :
1950 2 : CALL allocate_pw(v_laplaceb, pw_pool, bo)
1951 : END IF
1952 :
1953 : ELSE
1954 :
1955 : ! Create the work grid for the potential of the gradient part
1956 26 : CALL allocate_pw(v_drho, pw_pool, bo)
1957 :
1958 : CALL apply_drho(deriv_set, [deriv_norm_drho], virial_pw, drho, drho1, virial_xc, &
1959 26 : norm_drho, gradient_cut, dr1dr, v_drho%array)
1960 26 : IF (laplace_f) THEN
1961 2 : CALL xc_derivative_get(xc_dset_get_derivative(deriv_set, [deriv_laplace_rho]), deriv_data=deriv_data)
1962 2 : CPASSERT(ASSOCIATED(deriv_data))
1963 28862 : virial_pw%array(:, :, :) = -rho1(:, :, :)
1964 2 : CALL virial_laplace(virial_pw, pw_pool, virial_xc, deriv_data)
1965 :
1966 2 : CALL allocate_pw(v_laplace, pw_pool, bo)
1967 : END IF
1968 :
1969 : END IF
1970 :
1971 36 : IF (lsd) THEN
1972 150260 : rho_r(1)%array = rhoa
1973 150260 : rho_r(2)%array = rhob
1974 : ELSE
1975 701552 : rho_r(1)%array = rho
1976 : END IF
1977 36 : IF (ASSOCIATED(tau1_r)) THEN
1978 8 : IF (lsd) THEN
1979 60104 : tau_r(1)%array = tau_a
1980 60104 : tau_r(2)%array = tau_b
1981 : ELSE
1982 115448 : tau_r(1)%array = tau
1983 : END IF
1984 : END IF
1985 :
1986 : ! Create deriv sets with same densities but different gradients
1987 36 : CALL xc_dset_create(deriv_set1, pw_pool)
1988 :
1989 36 : rho_cutoff = section_get_rval(xc_section, "DENSITY_CUTOFF")
1990 :
1991 : ! create the place where to store the argument for the functionals
1992 : CALL xc_rho_set_create(rho2_set, bo, &
1993 : rho_cutoff=rho_cutoff, &
1994 : drho_cutoff=section_get_rval(xc_section, "GRADIENT_CUTOFF"), &
1995 36 : tau_cutoff=section_get_rval(xc_section, "TAU_CUTOFF"))
1996 :
1997 : ! calculate the arguments needed by the functionals
1998 : CALL xc_rho_set_update(rho2_set, rho_r, rho_g, tau_r, needs, &
1999 : section_get_ival(xc_section, "XC_GRID%XC_DERIV"), &
2000 : section_get_ival(xc_section, "XC_GRID%XC_SMOOTH_RHO"), &
2001 36 : pw_pool)
2002 :
2003 36 : IF (lsd) THEN
2004 : CALL xc_rho_set_get(rho1_set, rhoa=rho1a, rhob=rho1b, tau_a=tau1a, tau_b=tau1b, &
2005 10 : laplace_rhoa=laplace1a, laplace_rhob=laplace1b, can_return_null=.TRUE.)
2006 : CALL xc_rho_set_get(rho2_set, norm_drhoa=norm_drho2a, norm_drhob=norm_drho2b, &
2007 10 : norm_drho=norm_drho2, laplace_rhoa=laplace2a, laplace_rhob=laplace2b, can_return_null=.TRUE.)
2008 :
2009 64 : DO istep = -nsteps, nsteps
2010 54 : IF (istep == 0) CYCLE
2011 44 : weight = weights(istep, nsteps)/h
2012 44 : step = REAL(istep, dp)*h
2013 44 : IF (ASSOCIATED(norm_drhoa)) THEN
2014 44 : CALL get_derivs_rho(norm_drho2a, norm_drhoa, step, xc_fun_section, lsd, rho2_set, deriv_set1)
2015 : CALL update_deriv_rho(deriv_set1, [deriv_rhoa], bo, &
2016 44 : norm_drhoa, gradient_cut, weight, rho1a, v_drhoa%array)
2017 : CALL update_deriv_rho(deriv_set1, [deriv_rhob], bo, &
2018 44 : norm_drhoa, gradient_cut, weight, rho1b, v_drhoa%array)
2019 : CALL update_deriv_rho(deriv_set1, [deriv_norm_drhoa], bo, &
2020 44 : norm_drhoa, gradient_cut, weight, dra1dra, v_drhoa%array)
2021 : CALL update_deriv_drho_ab(deriv_set1, [deriv_norm_drhob], bo, &
2022 44 : norm_drhoa, gradient_cut, weight, dra1dra, drb1drb, v_drhoa%array, v_drhob%array)
2023 : CALL update_deriv_drho_ab(deriv_set1, [deriv_norm_drho], bo, &
2024 44 : norm_drhoa, gradient_cut, weight, dra1dra, dr1dr, v_drhoa%array, v_drho%array)
2025 44 : IF (tau_f) THEN
2026 : CALL update_deriv_rho(deriv_set1, [deriv_tau_a], bo, &
2027 8 : norm_drhoa, gradient_cut, weight, tau1a, v_drhoa%array)
2028 : CALL update_deriv_rho(deriv_set1, [deriv_tau_b], bo, &
2029 8 : norm_drhoa, gradient_cut, weight, tau1b, v_drhoa%array)
2030 : END IF
2031 44 : IF (laplace_f) THEN
2032 : CALL update_deriv_rho(deriv_set1, [deriv_laplace_rhoa], bo, &
2033 4 : norm_drhoa, gradient_cut, weight, laplace1a, v_drhoa%array)
2034 : CALL update_deriv_rho(deriv_set1, [deriv_laplace_rhob], bo, &
2035 4 : norm_drhoa, gradient_cut, weight, laplace1b, v_drhoa%array)
2036 : END IF
2037 : END IF
2038 :
2039 44 : IF (ASSOCIATED(norm_drhob)) THEN
2040 44 : CALL get_derivs_rho(norm_drho2b, norm_drhob, step, xc_fun_section, lsd, rho2_set, deriv_set1)
2041 : CALL update_deriv_rho(deriv_set1, [deriv_rhoa], bo, &
2042 44 : norm_drhob, gradient_cut, weight, rho1a, v_drhob%array)
2043 : CALL update_deriv_rho(deriv_set1, [deriv_rhob], bo, &
2044 44 : norm_drhob, gradient_cut, weight, rho1b, v_drhob%array)
2045 : CALL update_deriv_rho(deriv_set1, [deriv_norm_drhob], bo, &
2046 44 : norm_drhob, gradient_cut, weight, drb1drb, v_drhob%array)
2047 : CALL update_deriv_drho_ab(deriv_set1, [deriv_norm_drhoa], bo, &
2048 44 : norm_drhob, gradient_cut, weight, drb1drb, dra1dra, v_drhob%array, v_drhoa%array)
2049 : CALL update_deriv_drho_ab(deriv_set1, [deriv_norm_drho], bo, &
2050 44 : norm_drhob, gradient_cut, weight, drb1drb, dr1dr, v_drhob%array, v_drho%array)
2051 44 : IF (tau_f) THEN
2052 : CALL update_deriv_rho(deriv_set1, [deriv_tau_a], bo, &
2053 8 : norm_drhob, gradient_cut, weight, tau1a, v_drhob%array)
2054 : CALL update_deriv_rho(deriv_set1, [deriv_tau_b], bo, &
2055 8 : norm_drhob, gradient_cut, weight, tau1b, v_drhob%array)
2056 : END IF
2057 44 : IF (laplace_f) THEN
2058 : CALL update_deriv_rho(deriv_set1, [deriv_laplace_rhoa], bo, &
2059 4 : norm_drhob, gradient_cut, weight, laplace1a, v_drhob%array)
2060 : CALL update_deriv_rho(deriv_set1, [deriv_laplace_rhob], bo, &
2061 4 : norm_drhob, gradient_cut, weight, laplace1b, v_drhob%array)
2062 : END IF
2063 : END IF
2064 :
2065 44 : IF (ASSOCIATED(norm_drho)) THEN
2066 20 : CALL get_derivs_rho(norm_drho2, norm_drho, step, xc_fun_section, lsd, rho2_set, deriv_set1)
2067 : CALL update_deriv_rho(deriv_set1, [deriv_rhoa], bo, &
2068 20 : norm_drho, gradient_cut, weight, rho1a, v_drho%array)
2069 : CALL update_deriv_rho(deriv_set1, [deriv_rhob], bo, &
2070 20 : norm_drho, gradient_cut, weight, rho1b, v_drho%array)
2071 : CALL update_deriv_rho(deriv_set1, [deriv_norm_drho], bo, &
2072 20 : norm_drho, gradient_cut, weight, dr1dr, v_drho%array)
2073 : CALL update_deriv_drho_ab(deriv_set1, [deriv_norm_drhoa], bo, &
2074 20 : norm_drho, gradient_cut, weight, dr1dr, dra1dra, v_drho%array, v_drhoa%array)
2075 : CALL update_deriv_drho_ab(deriv_set1, [deriv_norm_drhob], bo, &
2076 20 : norm_drho, gradient_cut, weight, dr1dr, drb1drb, v_drho%array, v_drhob%array)
2077 20 : IF (tau_f) THEN
2078 : CALL update_deriv_rho(deriv_set1, [deriv_tau_a], bo, &
2079 8 : norm_drho, gradient_cut, weight, tau1a, v_drho%array)
2080 : CALL update_deriv_rho(deriv_set1, [deriv_tau_b], bo, &
2081 8 : norm_drho, gradient_cut, weight, tau1b, v_drho%array)
2082 : END IF
2083 20 : IF (laplace_f) THEN
2084 : CALL update_deriv_rho(deriv_set1, [deriv_laplace_rhoa], bo, &
2085 4 : norm_drho, gradient_cut, weight, laplace1a, v_drho%array)
2086 : CALL update_deriv_rho(deriv_set1, [deriv_laplace_rhob], bo, &
2087 4 : norm_drho, gradient_cut, weight, laplace1b, v_drho%array)
2088 : END IF
2089 : END IF
2090 :
2091 54 : IF (laplace_f) THEN
2092 :
2093 4 : CALL get_derivs_rho(laplace2a, laplacea, step, xc_fun_section, lsd, rho2_set, deriv_set1)
2094 :
2095 : ! Obtain the numerical 2nd derivatives w.r.t. to drho and collect the potential
2096 : CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [deriv_rhoa], bo, &
2097 4 : weight, rho1a, v_laplacea%array)
2098 : CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [deriv_rhob], bo, &
2099 4 : weight, rho1b, v_laplacea%array)
2100 4 : IF (ASSOCIATED(norm_drho)) THEN
2101 : CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [deriv_norm_drho], bo, &
2102 4 : weight, dr1dr, v_laplacea%array)
2103 : END IF
2104 4 : IF (ASSOCIATED(norm_drhoa)) THEN
2105 : CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [deriv_norm_drhoa], bo, &
2106 4 : weight, dra1dra, v_laplacea%array)
2107 : END IF
2108 4 : IF (ASSOCIATED(norm_drhob)) THEN
2109 : CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [deriv_norm_drhob], bo, &
2110 4 : weight, drb1drb, v_laplacea%array)
2111 : END IF
2112 :
2113 4 : IF (ASSOCIATED(tau1a)) THEN
2114 : CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [deriv_tau_a], bo, &
2115 4 : weight, tau1a, v_laplacea%array)
2116 : END IF
2117 4 : IF (ASSOCIATED(tau1b)) THEN
2118 : CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [deriv_tau_b], bo, &
2119 4 : weight, tau1b, v_laplacea%array)
2120 : END IF
2121 :
2122 : CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [deriv_laplace_rhoa], bo, &
2123 4 : weight, laplace1a, v_laplacea%array)
2124 :
2125 : CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [deriv_laplace_rhob], bo, &
2126 4 : weight, laplace1b, v_laplacea%array)
2127 :
2128 : ! The same for the beta spin
2129 4 : CALL get_derivs_rho(laplace2b, laplaceb, step, xc_fun_section, lsd, rho2_set, deriv_set1)
2130 :
2131 : ! Obtain the numerical 2nd derivatives w.r.t. to drho and collect the potential
2132 : CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [deriv_rhoa], bo, &
2133 4 : weight, rho1a, v_laplaceb%array)
2134 : CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [deriv_rhob], bo, &
2135 4 : weight, rho1b, v_laplaceb%array)
2136 4 : IF (ASSOCIATED(norm_drho)) THEN
2137 : CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [deriv_norm_drho], bo, &
2138 4 : weight, dr1dr, v_laplaceb%array)
2139 : END IF
2140 4 : IF (ASSOCIATED(norm_drhoa)) THEN
2141 : CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [deriv_norm_drhoa], bo, &
2142 4 : weight, dra1dra, v_laplaceb%array)
2143 : END IF
2144 4 : IF (ASSOCIATED(norm_drhob)) THEN
2145 : CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [deriv_norm_drhob], bo, &
2146 4 : weight, drb1drb, v_laplaceb%array)
2147 : END IF
2148 :
2149 4 : IF (tau_f) THEN
2150 : CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [deriv_tau_a], bo, &
2151 4 : weight, tau1a, v_laplaceb%array)
2152 : CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [deriv_tau_b], bo, &
2153 4 : weight, tau1b, v_laplaceb%array)
2154 : END IF
2155 :
2156 : CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [deriv_laplace_rhoa], bo, &
2157 4 : weight, laplace1a, v_laplaceb%array)
2158 :
2159 : CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [deriv_laplace_rhob], bo, &
2160 4 : weight, laplace1b, v_laplaceb%array)
2161 : END IF
2162 : END DO
2163 :
2164 10 : CALL virial_drho_drho(virial_pw, drhoa, v_drhoa, virial_xc)
2165 10 : CALL virial_drho_drho(virial_pw, drhob, v_drhob, virial_xc)
2166 10 : CALL virial_drho_drho(virial_pw, drho, v_drho, virial_xc)
2167 :
2168 10 : CALL deallocate_pw(v_drho, pw_pool)
2169 10 : CALL deallocate_pw(v_drhoa, pw_pool)
2170 10 : CALL deallocate_pw(v_drhob, pw_pool)
2171 :
2172 10 : IF (laplace_f) THEN
2173 15026 : virial_pw%array(:, :, :) = -rhoa(:, :, :)
2174 2 : CALL virial_laplace(virial_pw, pw_pool, virial_xc, v_laplacea%array)
2175 2 : CALL deallocate_pw(v_laplacea, pw_pool)
2176 :
2177 15026 : virial_pw%array(:, :, :) = -rhob(:, :, :)
2178 2 : CALL virial_laplace(virial_pw, pw_pool, virial_xc, v_laplaceb%array)
2179 2 : CALL deallocate_pw(v_laplaceb, pw_pool)
2180 : END IF
2181 :
2182 10 : CALL deallocate_pw(virial_pw, pw_pool)
2183 :
2184 40 : DO idir = 1, 3
2185 30 : DEALLOCATE (drho(idir)%array)
2186 40 : DEALLOCATE (drho1(idir)%array)
2187 : END DO
2188 10 : DEALLOCATE (dra1dra, drb1drb)
2189 :
2190 : ELSE
2191 26 : CALL xc_rho_set_get(rho1_set, rho=rho1, tau=tau1, laplace_rho=laplace1, can_return_null=.TRUE.)
2192 26 : CALL xc_rho_set_get(rho2_set, norm_drho=norm_drho2, laplace_rho=laplace2, can_return_null=.TRUE.)
2193 :
2194 200 : DO istep = -nsteps, nsteps
2195 174 : IF (istep == 0) CYCLE
2196 148 : weight = weights(istep, nsteps)/h
2197 148 : step = REAL(istep, dp)*h
2198 148 : CALL get_derivs_rho(norm_drho2, norm_drho, step, xc_fun_section, lsd, rho2_set, deriv_set1)
2199 :
2200 : ! Obtain the numerical 2nd derivatives w.r.t. to drho and collect the potential
2201 : CALL update_deriv_rho(deriv_set1, [deriv_rho], bo, &
2202 148 : norm_drho, gradient_cut, weight, rho1, v_drho%array)
2203 : CALL update_deriv_rho(deriv_set1, [deriv_norm_drho], bo, &
2204 148 : norm_drho, gradient_cut, weight, dr1dr, v_drho%array)
2205 :
2206 148 : IF (tau_f) THEN
2207 : CALL update_deriv_rho(deriv_set1, [deriv_tau], bo, &
2208 24 : norm_drho, gradient_cut, weight, tau1, v_drho%array)
2209 : END IF
2210 174 : IF (laplace_f) THEN
2211 : CALL update_deriv_rho(deriv_set1, [deriv_laplace_rho], bo, &
2212 12 : norm_drho, gradient_cut, weight, laplace1, v_drho%array)
2213 :
2214 12 : CALL get_derivs_rho(laplace2, laplace, step, xc_fun_section, lsd, rho2_set, deriv_set1)
2215 :
2216 : ! Obtain the numerical 2nd derivatives w.r.t. to drho and collect the potential
2217 : CALL update_deriv(deriv_set1, laplace, rho_cutoff, [deriv_rho], bo, &
2218 12 : weight, rho1, v_laplace%array)
2219 : CALL update_deriv(deriv_set1, laplace, rho_cutoff, [deriv_norm_drho], bo, &
2220 12 : weight, dr1dr, v_laplace%array)
2221 :
2222 12 : IF (tau_f) THEN
2223 : CALL update_deriv(deriv_set1, laplace, rho_cutoff, [deriv_tau], bo, &
2224 12 : weight, tau1, v_laplace%array)
2225 : END IF
2226 :
2227 : CALL update_deriv(deriv_set1, laplace, rho_cutoff, [deriv_laplace_rho], bo, &
2228 12 : weight, laplace1, v_laplace%array)
2229 : END IF
2230 : END DO
2231 :
2232 : ! Calculate the virial contribution from the potential
2233 26 : CALL virial_drho_drho(virial_pw, drho, v_drho, virial_xc)
2234 :
2235 26 : CALL deallocate_pw(v_drho, pw_pool)
2236 :
2237 26 : IF (laplace_f) THEN
2238 28862 : virial_pw%array(:, :, :) = -rho(:, :, :)
2239 2 : CALL virial_laplace(virial_pw, pw_pool, virial_xc, v_laplace%array)
2240 2 : CALL deallocate_pw(v_laplace, pw_pool)
2241 : END IF
2242 :
2243 26 : CALL deallocate_pw(virial_pw, pw_pool)
2244 : END IF
2245 :
2246 : END IF
2247 :
2248 36 : CALL xc_dset_release(deriv_set1)
2249 :
2250 36 : DEALLOCATE (dr1dr)
2251 :
2252 36 : CALL xc_rho_set_release(rho1_set)
2253 36 : CALL xc_rho_set_release(rho2_set)
2254 : END IF
2255 :
2256 2012 : DO ispin = 1, SIZE(rho_r)
2257 2012 : CALL pw_pool%give_back_pw(rho_r(ispin))
2258 : END DO
2259 832 : DEALLOCATE (rho_r)
2260 :
2261 832 : IF (ASSOCIATED(tau_r)) THEN
2262 254 : DO ispin = 1, SIZE(tau_r)
2263 254 : CALL pw_pool%give_back_pw(tau_r(ispin))
2264 : END DO
2265 100 : DEALLOCATE (tau_r)
2266 : END IF
2267 :
2268 832 : CALL timestop(handle)
2269 :
2270 37440 : END SUBROUTINE xc_calc_2nd_deriv_numerical
2271 :
2272 : ! **************************************************************************************************
2273 : !> \brief ...
2274 : !> \param rho_r ...
2275 : !> \param rho_g ...
2276 : !> \param rho1_r ...
2277 : !> \param rhoa ...
2278 : !> \param rhob ...
2279 : !> \param vxc_rho ...
2280 : !> \param tau_r ...
2281 : !> \param tau1_r ...
2282 : !> \param tau_a ...
2283 : !> \param tau_b ...
2284 : !> \param vxc_tau ...
2285 : !> \param xc_section ...
2286 : !> \param pw_pool ...
2287 : !> \param step ...
2288 : ! **************************************************************************************************
2289 1848 : SUBROUTINE calc_resp_potential_numer_ab(rho_r, rho_g, rho1_r, rhoa, rhob, vxc_rho, &
2290 : tau_r, tau1_r, tau_a, tau_b, vxc_tau, &
2291 : xc_section, pw_pool, step)
2292 :
2293 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER, INTENT(IN) :: vxc_rho, vxc_tau
2294 : TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(IN) :: rho1_r
2295 : TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(IN), POINTER :: tau1_r
2296 : TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool
2297 : TYPE(section_vals_type), INTENT(IN), POINTER :: xc_section
2298 : REAL(KIND=dp), INTENT(IN) :: step
2299 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER, INTENT(IN) :: rhoa, rhob, tau_a, tau_b
2300 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER, INTENT(IN) :: rho_r
2301 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
2302 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: tau_r
2303 :
2304 : CHARACTER(len=*), PARAMETER :: routineN = 'calc_resp_potential_numer_ab'
2305 :
2306 : INTEGER :: handle
2307 : REAL(KIND=dp) :: exc
2308 : REAL(KIND=dp), DIMENSION(3, 3) :: virial_dummy
2309 :
2310 1848 : CALL timeset(routineN, handle)
2311 :
2312 1848 : !$OMP PARALLEL DEFAULT(NONE) SHARED(rho_r,rhoa,rhob,step,rho1_r,tau_r,tau_a,tau_b,tau1_r)
2313 : !$OMP WORKSHARE
2314 : rho_r(1)%array(:, :, :) = rhoa(:, :, :) + step*rho1_r(1)%array(:, :, :)
2315 : !$OMP END WORKSHARE NOWAIT
2316 : !$OMP WORKSHARE
2317 : rho_r(2)%array(:, :, :) = rhob(:, :, :) + step*rho1_r(2)%array(:, :, :)
2318 : !$OMP END WORKSHARE NOWAIT
2319 : IF (ASSOCIATED(tau1_r) .AND. ASSOCIATED(tau_r) .AND. ASSOCIATED(tau_a) .AND. ASSOCIATED(tau_b)) THEN
2320 : !$OMP WORKSHARE
2321 : tau_r(1)%array(:, :, :) = tau_a(:, :, :) + step*tau1_r(1)%array(:, :, :)
2322 : !$OMP END WORKSHARE NOWAIT
2323 : !$OMP WORKSHARE
2324 : tau_r(2)%array(:, :, :) = tau_b(:, :, :) + step*tau1_r(2)%array(:, :, :)
2325 : !$OMP END WORKSHARE NOWAIT
2326 : END IF
2327 : !$OMP END PARALLEL
2328 : CALL xc_vxc_pw_create(vxc_rho, vxc_tau, exc, rho_r, rho_g, tau_r, xc_section, &
2329 1848 : pw_pool, .FALSE., virial_dummy)
2330 :
2331 1848 : CALL timestop(handle)
2332 :
2333 1848 : END SUBROUTINE calc_resp_potential_numer_ab
2334 :
2335 : ! **************************************************************************************************
2336 : !> \brief calculates stress tensor and potential contributions from the first derivative
2337 : !> \param deriv_set ...
2338 : !> \param description ...
2339 : !> \param virial_pw ...
2340 : !> \param drho ...
2341 : !> \param drho1 ...
2342 : !> \param virial_xc ...
2343 : !> \param norm_drho ...
2344 : !> \param gradient_cut ...
2345 : !> \param dr1dr ...
2346 : !> \param v_drho ...
2347 : ! **************************************************************************************************
2348 52 : SUBROUTINE apply_drho(deriv_set, description, virial_pw, drho, drho1, virial_xc, norm_drho, gradient_cut, dr1dr, v_drho)
2349 :
2350 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
2351 : INTEGER, DIMENSION(:), INTENT(in) :: description
2352 : TYPE(pw_r3d_rs_type), INTENT(IN) :: virial_pw
2353 : TYPE(cp_3d_r_cp_type), DIMENSION(3), INTENT(IN) :: drho, drho1
2354 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT) :: virial_xc
2355 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: norm_drho
2356 : REAL(KIND=dp), INTENT(IN) :: gradient_cut
2357 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: dr1dr
2358 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: v_drho
2359 :
2360 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_drho'
2361 :
2362 : INTEGER :: handle
2363 52 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: deriv_data
2364 : TYPE(xc_derivative_type), POINTER :: deriv_att
2365 :
2366 52 : CALL timeset(routineN, handle)
2367 :
2368 52 : deriv_att => xc_dset_get_derivative(deriv_set, description)
2369 52 : IF (ASSOCIATED(deriv_att)) THEN
2370 52 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
2371 52 : CALL virial_drho_drho1(virial_pw, drho, drho1, deriv_data, virial_xc)
2372 :
2373 52 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(dr1dr,gradient_cut,norm_drho,v_drho,deriv_data)
2374 : v_drho(:, :, :) = v_drho(:, :, :) + &
2375 : deriv_data(:, :, :)*dr1dr(:, :, :)/MAX(gradient_cut, norm_drho(:, :, :))**2
2376 : !$OMP END PARALLEL WORKSHARE
2377 : END IF
2378 :
2379 52 : CALL timestop(handle)
2380 :
2381 52 : END SUBROUTINE apply_drho
2382 :
2383 : ! **************************************************************************************************
2384 : !> \brief adds potential contributions from derivatives of rho or diagonal terms of norm_drho
2385 : !> \param deriv_set1 ...
2386 : !> \param description ...
2387 : !> \param bo ...
2388 : !> \param norm_drho norm_drho of which derivative is calculated
2389 : !> \param gradient_cut ...
2390 : !> \param h ...
2391 : !> \param rho1 function to contract the derivative with (rho1 for rho, dr1dr for norm_drho)
2392 : !> \param v_drho ...
2393 : ! **************************************************************************************************
2394 728 : SUBROUTINE update_deriv_rho(deriv_set1, description, bo, norm_drho, gradient_cut, weight, rho1, v_drho)
2395 :
2396 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set1
2397 : INTEGER, DIMENSION(:), INTENT(in) :: description
2398 : INTEGER, DIMENSION(2, 3), INTENT(IN) :: bo
2399 : REAL(KIND=dp), DIMENSION(bo(1, 1):bo(2, 1), bo(1, &
2400 : 2):bo(2, 2), bo(1, 3):bo(2, 3)), INTENT(IN) :: norm_drho
2401 : REAL(KIND=dp), INTENT(IN) :: gradient_cut, weight
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(IN) :: rho1
2404 : REAL(KIND=dp), DIMENSION(bo(1, 1):bo(2, 1), bo(1, &
2405 : 2):bo(2, 2), bo(1, 3):bo(2, 3)), INTENT(INOUT) :: v_drho
2406 :
2407 : CHARACTER(len=*), PARAMETER :: routineN = 'update_deriv_rho'
2408 :
2409 : INTEGER :: handle, i, j, k
2410 : REAL(KIND=dp) :: de
2411 728 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: deriv_data1
2412 : TYPE(xc_derivative_type), POINTER :: deriv_att1
2413 :
2414 728 : CALL timeset(routineN, handle)
2415 :
2416 : ! Obtain the numerical 2nd derivatives w.r.t. to drho and collect the potential
2417 728 : deriv_att1 => xc_dset_get_derivative(deriv_set1, description)
2418 728 : IF (ASSOCIATED(deriv_att1)) THEN
2419 728 : CALL xc_derivative_get(deriv_att1, deriv_data=deriv_data1)
2420 : !$OMP PARALLEL DO DEFAULT(NONE) SHARED(bo,deriv_data1,weight,norm_drho,v_drho,rho1,gradient_cut) &
2421 728 : !$OMP PRIVATE(i,j,k,de) 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) SHARED(bo,deriv_data1,weight,v,rho1,rho, rho_cutoff) &
2471 120 : !$OMP PRIVATE(i,j,k,de) COLLAPSE(3)
2472 : DO k = bo(1, 3), bo(2, 3)
2473 : DO j = bo(1, 2), bo(2, 2)
2474 : DO i = bo(1, 1), bo(2, 1)
2475 : ! We have to consider that the given density (mostly the Laplacian) may have positive and negative values
2476 : de = weight*deriv_data1(i, j, k)/SIGN(MAX(ABS(rho(i, j, k)), rho_cutoff), rho(i, j, k))
2477 : v(i, j, k) = v(i, j, k) + de*rho1(i, j, k)
2478 : END DO
2479 : END DO
2480 : END DO
2481 : !$OMP END PARALLEL DO
2482 : END IF
2483 :
2484 120 : CALL timestop(handle)
2485 :
2486 120 : END SUBROUTINE update_deriv
2487 :
2488 : ! **************************************************************************************************
2489 : !> \brief adds mixed derivatives of norm_drho
2490 : !> \param deriv_set1 ...
2491 : !> \param description ...
2492 : !> \param bo ...
2493 : !> \param norm_drhoa norm_drho of which derivatives is calculated
2494 : !> \param gradient_cut ...
2495 : !> \param h ...
2496 : !> \param dra1dra dr1dr corresponding to norm_drho
2497 : !> \param drb1drb ...
2498 : !> \param v_drhoa potential corresponding to norm_drho
2499 : !> \param v_drhob ...
2500 : ! **************************************************************************************************
2501 216 : SUBROUTINE update_deriv_drho_ab(deriv_set1, description, bo, &
2502 216 : norm_drhoa, gradient_cut, weight, dra1dra, drb1drb, v_drhoa, v_drhob)
2503 :
2504 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set1
2505 : INTEGER, DIMENSION(:), INTENT(in) :: description
2506 : INTEGER, DIMENSION(2, 3), INTENT(IN) :: bo
2507 : REAL(KIND=dp), DIMENSION(bo(1, 1):bo(2, 1), bo(1, &
2508 : 2):bo(2, 2), bo(1, 3):bo(2, 3)), INTENT(IN) :: norm_drhoa
2509 : REAL(KIND=dp), INTENT(IN) :: gradient_cut, weight
2510 : REAL(KIND=dp), DIMENSION(bo(1, 1):bo(2, 1), bo(1, &
2511 : 2):bo(2, 2), bo(1, 3):bo(2, 3)), INTENT(IN) :: dra1dra, drb1drb
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(INOUT) :: v_drhoa, v_drhob
2514 :
2515 : CHARACTER(len=*), PARAMETER :: routineN = 'update_deriv_drho_ab'
2516 :
2517 : INTEGER :: handle, i, j, k
2518 : REAL(KIND=dp) :: de
2519 216 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: deriv_data1
2520 : TYPE(xc_derivative_type), POINTER :: deriv_att1
2521 :
2522 216 : CALL timeset(routineN, handle)
2523 :
2524 216 : deriv_att1 => xc_dset_get_derivative(deriv_set1, description)
2525 216 : IF (ASSOCIATED(deriv_att1)) THEN
2526 168 : CALL xc_derivative_get(deriv_att1, deriv_data=deriv_data1)
2527 : !$OMP PARALLEL DO PRIVATE(k,j,i,de) DEFAULT(NONE) &
2528 168 : !$OMP SHARED(bo,drb1drb,dra1dra,deriv_data1,weight,gradient_cut,norm_drhoa,v_drhoa,v_drhob) COLLAPSE(3)
2529 : DO k = bo(1, 3), bo(2, 3)
2530 : DO j = bo(1, 2), bo(2, 2)
2531 : DO i = bo(1, 1), bo(2, 1)
2532 : ! We introduce a factor of two because we will average between both numerical derivatives
2533 : de = 0.5_dp*weight*deriv_data1(i, j, k)/MAX(gradient_cut, norm_drhoa(i, j, k))**2
2534 : v_drhoa(i, j, k) = v_drhoa(i, j, k) - de*drb1drb(i, j, k)
2535 : v_drhob(i, j, k) = v_drhob(i, j, k) - de*dra1dra(i, j, k)
2536 : END DO
2537 : END DO
2538 : END DO
2539 : END IF
2540 :
2541 216 : CALL timestop(handle)
2542 :
2543 216 : END SUBROUTINE update_deriv_drho_ab
2544 :
2545 : ! **************************************************************************************************
2546 : !> \brief calculate derivative sets for helper points
2547 : !> \param norm_drho2 norm_drho of new points
2548 : !> \param norm_drho norm_drho of KS density
2549 : !> \param h ...
2550 : !> \param xc_fun_section ...
2551 : !> \param lsd ...
2552 : !> \param rho2_set rho_set for new points
2553 : !> \param deriv_set1 will contain derivatives of the perturbed density
2554 : ! **************************************************************************************************
2555 276 : SUBROUTINE get_derivs_rho(norm_drho2, norm_drho, step, xc_fun_section, lsd, rho2_set, deriv_set1)
2556 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: norm_drho2
2557 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: norm_drho
2558 : REAL(KIND=dp), INTENT(IN) :: step
2559 : TYPE(section_vals_type), INTENT(IN), POINTER :: xc_fun_section
2560 : LOGICAL, INTENT(IN) :: lsd
2561 : TYPE(xc_rho_set_type), INTENT(INOUT) :: rho2_set
2562 : TYPE(xc_derivative_set_type) :: deriv_set1
2563 :
2564 : CHARACTER(len=*), PARAMETER :: routineN = 'get_derivs_rho'
2565 :
2566 : INTEGER :: handle
2567 :
2568 276 : CALL timeset(routineN, handle)
2569 :
2570 : ! Copy the densities, do one step into the direction of drho
2571 276 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(norm_drho,norm_drho2,step)
2572 : norm_drho2 = norm_drho*(1.0_dp + step)
2573 : !$OMP END PARALLEL WORKSHARE
2574 :
2575 276 : CALL xc_dset_zero_all(deriv_set1)
2576 :
2577 : ! Calculate the derivatives of the functional
2578 : CALL xc_functionals_eval(xc_fun_section, &
2579 : lsd=lsd, &
2580 : rho_set=rho2_set, &
2581 : deriv_set=deriv_set1, &
2582 276 : deriv_order=1)
2583 :
2584 : ! Return to the original values
2585 276 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(norm_drho,norm_drho2)
2586 : norm_drho2 = norm_drho
2587 : !$OMP END PARALLEL WORKSHARE
2588 :
2589 276 : CALL divide_by_norm_drho(deriv_set1, rho2_set, lsd)
2590 :
2591 276 : CALL timestop(handle)
2592 :
2593 276 : END SUBROUTINE get_derivs_rho
2594 :
2595 : ! **************************************************************************************************
2596 : !> \brief Calculates the second derivative of E_xc at rho in the direction
2597 : !> rho1 (if you see the second derivative as bilinear form)
2598 : !> partial_rho|_(rho=rho) partial_rho|_(rho=rho) E_xc drho(rho1)drho
2599 : !> The other direction is still undetermined, thus it returns
2600 : !> a potential (partial integration is performed to reduce it to
2601 : !> function of rho, removing the dependence from its partial derivs)
2602 : !> Has to be called after the setup by xc_prep_2nd_deriv.
2603 : !> \param v_xc exchange-correlation potential
2604 : !> \param v_xc_tau ...
2605 : !> \param deriv_set derivatives of the exchange-correlation potential
2606 : !> \param rho_set object containing the density at which the derivatives were calculated
2607 : !> \param rho1_set object containing the density with which to fold
2608 : !> \param pw_pool the pool for the grids
2609 : !> \param xc_section XC parameters
2610 : !> \param gapw Gaussian and augmented plane waves calculation
2611 : !> \param vxg ...
2612 : !> \param tddfpt_fac factor that multiplies the crossterms (tddfpt triplets
2613 : !> on a closed shell system it should be -1, defaults to 1)
2614 : !> \param compute_virial ...
2615 : !> \param virial_xc ...
2616 : !> \note
2617 : !> The old version of this routine was smarter: it handled split_desc(1)
2618 : !> and split_desc(2) separately, thus the code automatically handled all
2619 : !> possible cross terms (you only had to check if it was diagonal to avoid
2620 : !> double counting). I think that is the way to go if you want to add more
2621 : !> terms (tau,rho in LSD,...). The problem with the old code was that it
2622 : !> because of the old functional structure it sometime guessed wrongly
2623 : !> which derivative was where. There were probably still bugs with gradient
2624 : !> corrected functionals (never tested), and it didn't contain first
2625 : !> derivatives with respect to drho (that contribute also to the second
2626 : !> derivative wrt. rho).
2627 : !> The code was a little complex because it really tried to handle any
2628 : !> functional derivative in the most efficient way with the given contents of
2629 : !> rho_set.
2630 : !> Anyway I strongly encourage whoever wants to modify this code to give a
2631 : !> look to the old version. [fawzi]
2632 : ! **************************************************************************************************
2633 31304 : SUBROUTINE xc_calc_2nd_deriv_analytical(v_xc, v_xc_tau, deriv_set, rho_set, rho1_set, &
2634 : pw_pool, xc_section, gapw, vxg, tddfpt_fac, &
2635 : compute_virial, virial_xc)
2636 :
2637 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_xc, v_xc_tau
2638 : TYPE(xc_derivative_set_type) :: deriv_set
2639 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set, rho1_set
2640 : TYPE(pw_pool_type), POINTER :: pw_pool
2641 : TYPE(section_vals_type), POINTER :: xc_section
2642 : LOGICAL, INTENT(IN), OPTIONAL :: gapw
2643 : REAL(kind=dp), DIMENSION(:, :, :, :), OPTIONAL, &
2644 : POINTER :: vxg
2645 : REAL(kind=dp), INTENT(in), OPTIONAL :: tddfpt_fac
2646 : LOGICAL, INTENT(IN), OPTIONAL :: compute_virial
2647 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT), &
2648 : OPTIONAL :: virial_xc
2649 :
2650 : CHARACTER(len=*), PARAMETER :: routineN = 'xc_calc_2nd_deriv_analytical'
2651 :
2652 : INTEGER :: handle, i, ia, idir, ir, ispin, j, jdir, &
2653 : k, nspins, xc_deriv_method_id
2654 : INTEGER, DIMENSION(2, 3) :: bo
2655 : LOGICAL :: gradient_f, lsd, my_compute_virial, &
2656 : my_gapw, tau_f, laplace_f, rho_f
2657 : REAL(KIND=dp) :: fac, gradient_cut, tmp, factor2
2658 31304 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: dr1dr, dra1dra, drb1drb
2659 62608 : REAL(kind=dp), DIMENSION(:, :, :), POINTER :: deriv_data, e_drhoa, e_drhob, &
2660 31304 : e_drho, norm_drho, norm_drhoa, &
2661 31304 : norm_drhob, rho1, rho1a, rho1b, &
2662 31304 : tau1, tau1a, tau1b, laplace1, laplace1a, laplace1b, &
2663 31304 : rho, rhoa, rhob
2664 594776 : TYPE(cp_3d_r_cp_type), DIMENSION(3) :: drho, drho1, drho1a, drho1b, drhoa, drhob
2665 31304 : TYPE(pw_r3d_rs_type), DIMENSION(:), ALLOCATABLE :: v_drhoa, v_drhob, v_drho, v_laplace
2666 31304 : TYPE(pw_r3d_rs_type), DIMENSION(:, :), ALLOCATABLE :: v_drho_r
2667 : TYPE(pw_r3d_rs_type) :: virial_pw
2668 : TYPE(pw_c1d_gs_type) :: tmp_g, vxc_g
2669 : TYPE(xc_derivative_type), POINTER :: deriv_att
2670 :
2671 31304 : CALL timeset(routineN, handle)
2672 :
2673 31304 : NULLIFY (e_drhoa, e_drhob, e_drho)
2674 :
2675 31304 : my_gapw = .FALSE.
2676 31304 : IF (PRESENT(gapw)) my_gapw = gapw
2677 :
2678 31304 : my_compute_virial = .FALSE.
2679 31304 : IF (PRESENT(compute_virial)) my_compute_virial = compute_virial
2680 :
2681 31304 : CPASSERT(ASSOCIATED(v_xc))
2682 31304 : CPASSERT(ASSOCIATED(xc_section))
2683 31304 : IF (my_gapw) THEN
2684 10060 : CPASSERT(PRESENT(vxg))
2685 : END IF
2686 31304 : IF (my_compute_virial) THEN
2687 354 : CPASSERT(PRESENT(virial_xc))
2688 : END IF
2689 :
2690 : CALL section_vals_val_get(xc_section, "XC_GRID%XC_DERIV", &
2691 31304 : i_val=xc_deriv_method_id)
2692 31304 : CALL xc_rho_set_get(rho_set, drho_cutoff=gradient_cut)
2693 31304 : nspins = SIZE(v_xc)
2694 31304 : lsd = ASSOCIATED(rho_set%rhoa)
2695 31304 : fac = 0.0_dp
2696 31304 : factor2 = 1.0_dp
2697 31304 : IF (PRESENT(tddfpt_fac)) fac = tddfpt_fac
2698 31304 : IF (PRESENT(tddfpt_fac)) factor2 = tddfpt_fac
2699 :
2700 313040 : bo = rho_set%local_bounds
2701 :
2702 31304 : CALL check_for_derivatives(deriv_set, lsd, rho_f, gradient_f, tau_f, laplace_f)
2703 :
2704 31304 : IF (tau_f) THEN
2705 614 : CPASSERT(ASSOCIATED(v_xc_tau))
2706 : END IF
2707 :
2708 31304 : IF (gradient_f) THEN
2709 203810 : ALLOCATE (v_drho_r(3, nspins), v_drho(nspins))
2710 40762 : DO ispin = 1, nspins
2711 84136 : DO idir = 1, 3
2712 84136 : CALL allocate_pw(v_drho_r(idir, ispin), pw_pool, bo)
2713 : END DO
2714 40762 : CALL allocate_pw(v_drho(ispin), pw_pool, bo)
2715 : END DO
2716 :
2717 19728 : IF (xc_requires_tmp_g(xc_deriv_method_id) .AND. .NOT. my_gapw) THEN
2718 12034 : IF (ASSOCIATED(pw_pool)) THEN
2719 12034 : CALL pw_pool%create_pw(tmp_g)
2720 12034 : CALL pw_pool%create_pw(vxc_g)
2721 : ELSE
2722 : ! remember to refix for gapw
2723 0 : CPABORT("XC_DERIV method is not implemented in GAPW")
2724 : END IF
2725 : END IF
2726 : END IF
2727 :
2728 65826 : DO ispin = 1, nspins
2729 965751470 : v_xc(ispin)%array = 0.0_dp
2730 : END DO
2731 :
2732 31304 : IF (tau_f) THEN
2733 1290 : DO ispin = 1, nspins
2734 14079874 : v_xc_tau(ispin)%array = 0.0_dp
2735 : END DO
2736 : END IF
2737 :
2738 31304 : IF (laplace_f .AND. my_gapw) &
2739 0 : CPABORT("Laplace-dependent functional not implemented with GAPW!")
2740 :
2741 31304 : IF (my_compute_virial .AND. (gradient_f .OR. laplace_f)) CALL allocate_pw(virial_pw, pw_pool, bo)
2742 :
2743 31304 : IF (lsd) THEN
2744 :
2745 : !-------------------!
2746 : ! UNrestricted case !
2747 : !-------------------!
2748 :
2749 4512 : CALL xc_rho_set_get(rho1_set, rhoa=rho1a, rhob=rho1b)
2750 :
2751 4512 : IF (gradient_f) THEN
2752 : CALL xc_rho_set_get(rho_set, drhoa=drhoa, drhob=drhob, &
2753 2176 : norm_drho=norm_drho, norm_drhoa=norm_drhoa, norm_drhob=norm_drhob)
2754 2176 : CALL xc_rho_set_get(rho1_set, drhoa=drho1a, drhob=drho1b)
2755 :
2756 2176 : CALL calc_drho_from_ab(drho, drhoa, drhob)
2757 2176 : CALL calc_drho_from_ab(drho1, drho1a, drho1b)
2758 :
2759 2176 : CALL prepare_dr1dr(dra1dra, drhoa, drho1a)
2760 2176 : IF (nspins /= 1) THEN
2761 1306 : CALL prepare_dr1dr(drb1drb, drhob, drho1b)
2762 1306 : CALL prepare_dr1dr(dr1dr, drho, drho1)
2763 : ELSE
2764 870 : CALL prepare_dr1dr(drb1drb, drhob, drho1b)
2765 870 : CALL prepare_dr1dr_ab(dr1dr, drhoa, drhob, drho1a, drho1b, fac)
2766 : END IF
2767 :
2768 17844 : ALLOCATE (v_drhoa(nspins), v_drhob(nspins))
2769 5658 : DO ispin = 1, nspins
2770 3482 : CALL allocate_pw(v_drhoa(ispin), pw_pool, bo)
2771 5658 : CALL allocate_pw(v_drhob(ispin), pw_pool, bo)
2772 : END DO
2773 :
2774 : END IF
2775 :
2776 4512 : IF (laplace_f) THEN
2777 26 : CALL xc_rho_set_get(rho1_set, laplace_rhoa=laplace1a, laplace_rhob=laplace1b)
2778 :
2779 130 : ALLOCATE (v_laplace(nspins))
2780 78 : DO ispin = 1, nspins
2781 78 : CALL allocate_pw(v_laplace(ispin), pw_pool, bo)
2782 : END DO
2783 :
2784 26 : IF (my_compute_virial) CALL xc_rho_set_get(rho_set, rhoa=rhoa, rhob=rhob)
2785 : END IF
2786 :
2787 4512 : IF (tau_f) THEN
2788 62 : CALL xc_rho_set_get(rho1_set, tau_a=tau1a, tau_b=tau1b)
2789 : END IF
2790 :
2791 4512 : IF (nspins /= 1) THEN
2792 :
2793 3218 : $:add_2nd_derivative_terms(arguments_openshell)
2794 :
2795 : ELSE
2796 :
2797 1294 : $:add_2nd_derivative_terms(arguments_triplet_outer, arguments_triplet_inner)
2798 :
2799 : END IF
2800 :
2801 4512 : IF (gradient_f) THEN
2802 :
2803 2176 : IF (my_compute_virial) THEN
2804 8 : CALL virial_drho_drho(virial_pw, drhoa, v_drhoa(1), virial_xc)
2805 8 : CALL virial_drho_drho(virial_pw, drhob, v_drhob(2), virial_xc)
2806 32 : DO idir = 1, 3
2807 24 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(drho,idir,v_drho,virial_pw)
2808 : virial_pw%array(:, :, :) = drho(idir)%array(:, :, :)*(v_drho(1)%array(:, :, :) + v_drho(2)%array(:, :, :))
2809 : !$OMP END PARALLEL WORKSHARE
2810 80 : DO jdir = 1, idir
2811 : tmp = -0.5_dp*virial_pw%pw_grid%dvol*accurate_dot_product(virial_pw%array(:, :, :), &
2812 48 : drho(jdir)%array(:, :, :))
2813 48 : virial_xc(jdir, idir) = virial_xc(jdir, idir) + tmp
2814 72 : virial_xc(idir, jdir) = virial_xc(jdir, idir)
2815 : END DO
2816 : END DO
2817 : END IF ! my_compute_virial
2818 :
2819 2176 : IF (my_gapw) THEN
2820 : !$OMP PARALLEL DO PRIVATE(ia,idir,ispin,ir) DEFAULT(NONE) &
2821 : !$OMP SHARED(bo,nspins,vxg,drhoa,drhob,v_drhoa,v_drhob,v_drho, &
2822 602 : !$OMP e_drhoa,e_drhob,e_drho,drho1a,drho1b,fac,drho,drho1) COLLAPSE(3)
2823 : DO ir = bo(1, 2), bo(2, 2)
2824 : DO ia = bo(1, 1), bo(2, 1)
2825 : DO idir = 1, 3
2826 : DO ispin = 1, nspins
2827 : vxg(idir, ia, ir, ispin) = &
2828 : -(v_drhoa(ispin)%array(ia, ir, 1)*drhoa(idir)%array(ia, ir, 1) + &
2829 : v_drhob(ispin)%array(ia, ir, 1)*drhob(idir)%array(ia, ir, 1) + &
2830 : v_drho(ispin)%array(ia, ir, 1)*drho(idir)%array(ia, ir, 1))
2831 : END DO
2832 : IF (ASSOCIATED(e_drhoa)) THEN
2833 : vxg(idir, ia, ir, 1) = vxg(idir, ia, ir, 1) + &
2834 : e_drhoa(ia, ir, 1)*drho1a(idir)%array(ia, ir, 1)
2835 : END IF
2836 : IF (nspins /= 1 .AND. ASSOCIATED(e_drhob)) THEN
2837 : vxg(idir, ia, ir, 2) = vxg(idir, ia, ir, 2) + &
2838 : e_drhob(ia, ir, 1)*drho1b(idir)%array(ia, ir, 1)
2839 : END IF
2840 : IF (ASSOCIATED(e_drho)) THEN
2841 : IF (nspins /= 1) THEN
2842 : vxg(idir, ia, ir, 1) = vxg(idir, ia, ir, 1) + &
2843 : e_drho(ia, ir, 1)*drho1(idir)%array(ia, ir, 1)
2844 : vxg(idir, ia, ir, 2) = vxg(idir, ia, ir, 2) + &
2845 : e_drho(ia, ir, 1)*drho1(idir)%array(ia, ir, 1)
2846 : ELSE
2847 : vxg(idir, ia, ir, 1) = vxg(idir, ia, ir, 1) + &
2848 : e_drho(ia, ir, 1)*(drho1a(idir)%array(ia, ir, 1) + &
2849 : fac*drho1b(idir)%array(ia, ir, 1))
2850 : END IF
2851 : END IF
2852 : END DO
2853 : END DO
2854 : END DO
2855 : ELSE
2856 :
2857 : ! partial integration
2858 6296 : DO idir = 1, 3
2859 :
2860 12978 : DO ispin = 1, nspins
2861 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) &
2862 12978 : !$OMP SHARED(v_drho_r,v_drhoa,v_drhob,v_drho,drhoa,drhob,drho,ispin,idir)
2863 : v_drho_r(idir, ispin)%array(:, :, :) = &
2864 : v_drhoa(ispin)%array(:, :, :)*drhoa(idir)%array(:, :, :) + &
2865 : v_drhob(ispin)%array(:, :, :)*drhob(idir)%array(:, :, :) + &
2866 : v_drho(ispin)%array(:, :, :)*drho(idir)%array(:, :, :)
2867 : !$OMP END PARALLEL WORKSHARE
2868 : END DO
2869 4722 : IF (ASSOCIATED(e_drhoa)) THEN
2870 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) &
2871 4722 : !$OMP SHARED(v_drho_r,e_drhoa,drho1a,idir)
2872 : v_drho_r(idir, 1)%array(:, :, :) = v_drho_r(idir, 1)%array(:, :, :) - &
2873 : e_drhoa(:, :, :)*drho1a(idir)%array(:, :, :)
2874 : !$OMP END PARALLEL WORKSHARE
2875 : END IF
2876 4722 : IF (nspins /= 1 .AND. ASSOCIATED(e_drhob)) THEN
2877 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE)&
2878 3534 : !$OMP SHARED(v_drho_r,e_drhob,drho1b,idir)
2879 : v_drho_r(idir, 2)%array(:, :, :) = v_drho_r(idir, 2)%array(:, :, :) - &
2880 : e_drhob(:, :, :)*drho1b(idir)%array(:, :, :)
2881 : !$OMP END PARALLEL WORKSHARE
2882 : END IF
2883 6296 : IF (ASSOCIATED(e_drho)) THEN
2884 : !$OMP PARALLEL DO PRIVATE(k,j,i) DEFAULT(NONE)&
2885 4722 : !$OMP SHARED(bo,v_drho_r,e_drho,drho1a,drho1b,drho1,fac,idir,nspins) COLLAPSE(3)
2886 : DO k = bo(1, 3), bo(2, 3)
2887 : DO j = bo(1, 2), bo(2, 2)
2888 : DO i = bo(1, 1), bo(2, 1)
2889 : IF (nspins /= 1) THEN
2890 : v_drho_r(idir, 1)%array(i, j, k) = v_drho_r(idir, 1)%array(i, j, k) - &
2891 : e_drho(i, j, k)*drho1(idir)%array(i, j, k)
2892 : v_drho_r(idir, 2)%array(i, j, k) = v_drho_r(idir, 2)%array(i, j, k) - &
2893 : e_drho(i, j, k)*drho1(idir)%array(i, j, k)
2894 : ELSE
2895 : v_drho_r(idir, 1)%array(i, j, k) = v_drho_r(idir, 1)%array(i, j, k) - &
2896 : e_drho(i, j, k)*(drho1a(idir)%array(i, j, k) + &
2897 : fac*drho1b(idir)%array(i, j, k))
2898 : END IF
2899 : END DO
2900 : END DO
2901 : END DO
2902 : END IF
2903 : END DO
2904 :
2905 4326 : DO ispin = 1, nspins
2906 : ! partial integration
2907 4326 : CALL xc_pw_divergence(xc_deriv_method_id, v_drho_r(:, ispin), tmp_g, vxc_g, v_xc(ispin))
2908 : END DO ! ispin
2909 :
2910 : END IF
2911 :
2912 8704 : DO idir = 1, 3
2913 6528 : DEALLOCATE (drho(idir)%array)
2914 8704 : DEALLOCATE (drho1(idir)%array)
2915 : END DO
2916 :
2917 5658 : DO ispin = 1, nspins
2918 3482 : CALL deallocate_pw(v_drhoa(ispin), pw_pool)
2919 5658 : CALL deallocate_pw(v_drhob(ispin), pw_pool)
2920 : END DO
2921 :
2922 2176 : DEALLOCATE (v_drhoa, v_drhob)
2923 :
2924 : END IF ! gradient_f
2925 :
2926 4512 : IF (laplace_f .AND. my_compute_virial) THEN
2927 0 : virial_pw%array(:, :, :) = -rhoa(:, :, :)
2928 0 : CALL virial_laplace(virial_pw, pw_pool, virial_xc, v_laplace(1)%array)
2929 0 : virial_pw%array(:, :, :) = -rhob(:, :, :)
2930 0 : CALL virial_laplace(virial_pw, pw_pool, virial_xc, v_laplace(2)%array)
2931 : END IF
2932 :
2933 : ELSE
2934 :
2935 : !-----------------!
2936 : ! restricted case !
2937 : !-----------------!
2938 :
2939 26792 : CALL xc_rho_set_get(rho1_set, rho=rho1)
2940 :
2941 26792 : IF (gradient_f) THEN
2942 17552 : CALL xc_rho_set_get(rho_set, drho=drho, norm_drho=norm_drho)
2943 17552 : CALL xc_rho_set_get(rho1_set, drho=drho1)
2944 17552 : CALL prepare_dr1dr(dr1dr, drho, drho1)
2945 : END IF
2946 :
2947 26792 : IF (laplace_f) THEN
2948 126 : CALL xc_rho_set_get(rho1_set, laplace_rho=laplace1)
2949 :
2950 504 : ALLOCATE (v_laplace(nspins))
2951 252 : DO ispin = 1, nspins
2952 252 : CALL allocate_pw(v_laplace(ispin), pw_pool, bo)
2953 : END DO
2954 :
2955 126 : IF (my_compute_virial) CALL xc_rho_set_get(rho_set, rho=rho)
2956 : END IF
2957 :
2958 26792 : IF (tau_f) THEN
2959 552 : CALL xc_rho_set_get(rho1_set, tau=tau1)
2960 : END IF
2961 :
2962 292448 : $:add_2nd_derivative_terms(arguments_closedshell)
2963 :
2964 26792 : IF (gradient_f) THEN
2965 :
2966 17552 : IF (my_compute_virial) THEN
2967 220 : CALL virial_drho_drho(virial_pw, drho, v_drho(1), virial_xc)
2968 : END IF ! my_compute_virial
2969 :
2970 17552 : IF (my_gapw) THEN
2971 :
2972 26640 : DO idir = 1, 3
2973 : !$OMP PARALLEL DO PRIVATE(ia,ir) DEFAULT(NONE) &
2974 26640 : !$OMP SHARED(bo,vxg,drho,v_drho,e_drho,drho1,idir,factor2) COLLAPSE(2)
2975 : DO ia = bo(1, 1), bo(2, 1)
2976 : DO ir = bo(1, 2), bo(2, 2)
2977 : vxg(idir, ia, ir, 1) = -drho(idir)%array(ia, ir, 1)*v_drho(1)%array(ia, ir, 1)
2978 : IF (ASSOCIATED(e_drho)) THEN
2979 : vxg(idir, ia, ir, 1) = vxg(idir, ia, ir, 1) + factor2*drho1(idir)%array(ia, ir, 1)*e_drho(ia, ir, 1)
2980 : END IF
2981 : END DO
2982 : END DO
2983 : END DO
2984 :
2985 : ELSE
2986 : ! partial integration
2987 43568 : DO idir = 1, 3
2988 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE)&
2989 43568 : !$OMP SHARED(v_drho_r,drho,v_drho,drho1,e_drho,idir)
2990 : v_drho_r(idir, 1)%array(:, :, :) = drho(idir)%array(:, :, :)*v_drho(1)%array(:, :, :) - &
2991 : drho1(idir)%array(:, :, :)*e_drho(:, :, :)
2992 : !$OMP END PARALLEL WORKSHARE
2993 : END DO
2994 :
2995 10892 : CALL xc_pw_divergence(xc_deriv_method_id, v_drho_r(:, 1), tmp_g, vxc_g, v_xc(1))
2996 : END IF
2997 :
2998 : END IF
2999 :
3000 26792 : IF (laplace_f .AND. my_compute_virial) THEN
3001 265668 : virial_pw%array(:, :, :) = -rho(:, :, :)
3002 12 : CALL virial_laplace(virial_pw, pw_pool, virial_xc, v_laplace(1)%array)
3003 : END IF
3004 :
3005 : END IF
3006 :
3007 31304 : IF (laplace_f) THEN
3008 330 : DO ispin = 1, nspins
3009 178 : CALL xc_pw_laplace(v_laplace(ispin), pw_pool, xc_deriv_method_id)
3010 330 : CALL pw_axpy(v_laplace(ispin), v_xc(ispin))
3011 : END DO
3012 : END IF
3013 :
3014 31304 : IF (gradient_f) THEN
3015 :
3016 40762 : DO ispin = 1, nspins
3017 21034 : CALL deallocate_pw(v_drho(ispin), pw_pool)
3018 103864 : DO idir = 1, 3
3019 84136 : CALL deallocate_pw(v_drho_r(idir, ispin), pw_pool)
3020 : END DO
3021 : END DO
3022 19728 : DEALLOCATE (v_drho, v_drho_r)
3023 :
3024 : END IF
3025 :
3026 31304 : IF (laplace_f) THEN
3027 330 : DO ispin = 1, nspins
3028 330 : CALL deallocate_pw(v_laplace(ispin), pw_pool)
3029 : END DO
3030 152 : DEALLOCATE (v_laplace)
3031 : END IF
3032 :
3033 31304 : IF (ASSOCIATED(tmp_g%pw_grid) .AND. ASSOCIATED(pw_pool)) THEN
3034 12034 : CALL pw_pool%give_back_pw(tmp_g)
3035 : END IF
3036 :
3037 31304 : IF (ASSOCIATED(vxc_g%pw_grid) .AND. ASSOCIATED(pw_pool)) THEN
3038 12034 : CALL pw_pool%give_back_pw(vxc_g)
3039 : END IF
3040 :
3041 31304 : IF (my_compute_virial .AND. (gradient_f .OR. laplace_f)) THEN
3042 228 : CALL deallocate_pw(virial_pw, pw_pool)
3043 : END IF
3044 :
3045 31304 : CALL timestop(handle)
3046 93912 : END SUBROUTINE xc_calc_2nd_deriv_analytical
3047 :
3048 : ! **************************************************************************************************
3049 : !> \brief allocates grids using pw_pool (if associated) or with bounds
3050 : !> \param pw ...
3051 : !> \param pw_pool ...
3052 : !> \param bo ...
3053 : ! **************************************************************************************************
3054 91604 : SUBROUTINE allocate_pw(pw, pw_pool, bo)
3055 : TYPE(pw_r3d_rs_type), INTENT(OUT) :: pw
3056 : TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool
3057 : INTEGER, DIMENSION(2, 3), INTENT(IN) :: bo
3058 :
3059 91604 : IF (ASSOCIATED(pw_pool)) THEN
3060 60584 : CALL pw_pool%create_pw(pw)
3061 60584 : CALL pw_zero(pw)
3062 : ELSE
3063 155100 : ALLOCATE (pw%array(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
3064 79163040 : pw%array = 0.0_dp
3065 : END IF
3066 :
3067 91604 : END SUBROUTINE allocate_pw
3068 :
3069 : ! **************************************************************************************************
3070 : !> \brief deallocates grid allocated with allocate_pw
3071 : !> \param pw ...
3072 : !> \param pw_pool ...
3073 : ! **************************************************************************************************
3074 91604 : SUBROUTINE deallocate_pw(pw, pw_pool)
3075 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: pw
3076 : TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool
3077 :
3078 91604 : IF (ASSOCIATED(pw_pool)) THEN
3079 60584 : CALL pw_pool%give_back_pw(pw)
3080 : ELSE
3081 31020 : CALL pw%release()
3082 : END IF
3083 :
3084 91604 : END SUBROUTINE deallocate_pw
3085 :
3086 : ! **************************************************************************************************
3087 : !> \brief updates virial from first derivative w.r.t. norm_drho
3088 : !> \param virial_pw ...
3089 : !> \param drho ...
3090 : !> \param drho1 ...
3091 : !> \param deriv_data ...
3092 : !> \param virial_xc ...
3093 : ! **************************************************************************************************
3094 296 : SUBROUTINE virial_drho_drho1(virial_pw, drho, drho1, deriv_data, virial_xc)
3095 : TYPE(pw_r3d_rs_type), INTENT(IN) :: virial_pw
3096 : TYPE(cp_3d_r_cp_type), DIMENSION(3), INTENT(IN) :: drho, drho1
3097 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: deriv_data
3098 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT) :: virial_xc
3099 :
3100 : INTEGER :: idir, jdir
3101 : REAL(KIND=dp) :: tmp
3102 :
3103 1184 : DO idir = 1, 3
3104 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE)&
3105 888 : !$OMP SHARED(drho,idir,virial_pw,deriv_data)
3106 : virial_pw%array(:, :, :) = drho(idir)%array(:, :, :)*deriv_data(:, :, :)
3107 : !$OMP END PARALLEL WORKSHARE
3108 3848 : DO jdir = 1, 3
3109 : tmp = virial_pw%pw_grid%dvol*accurate_dot_product(virial_pw%array(:, :, :), &
3110 2664 : drho1(jdir)%array(:, :, :))
3111 2664 : virial_xc(jdir, idir) = virial_xc(jdir, idir) + tmp
3112 3552 : virial_xc(idir, jdir) = virial_xc(idir, jdir) + tmp
3113 : END DO
3114 : END DO
3115 :
3116 296 : END SUBROUTINE virial_drho_drho1
3117 :
3118 : ! **************************************************************************************************
3119 : !> \brief Adds virial contribution from second order potential parts
3120 : !> \param virial_pw ...
3121 : !> \param drho ...
3122 : !> \param v_drho ...
3123 : !> \param virial_xc ...
3124 : ! **************************************************************************************************
3125 292 : SUBROUTINE virial_drho_drho(virial_pw, drho, v_drho, virial_xc)
3126 : TYPE(pw_r3d_rs_type), INTENT(IN) :: virial_pw
3127 : TYPE(cp_3d_r_cp_type), DIMENSION(3), INTENT(IN) :: drho
3128 : TYPE(pw_r3d_rs_type), INTENT(IN) :: v_drho
3129 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT) :: virial_xc
3130 :
3131 : INTEGER :: idir, jdir
3132 : REAL(KIND=dp) :: tmp
3133 :
3134 1168 : DO idir = 1, 3
3135 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE)&
3136 876 : !$OMP SHARED(drho,idir,v_drho,virial_pw)
3137 : virial_pw%array(:, :, :) = drho(idir)%array(:, :, :)*v_drho%array(:, :, :)
3138 : !$OMP END PARALLEL WORKSHARE
3139 2920 : DO jdir = 1, idir
3140 : tmp = -virial_pw%pw_grid%dvol*accurate_dot_product(virial_pw%array(:, :, :), &
3141 1752 : drho(jdir)%array(:, :, :))
3142 1752 : virial_xc(jdir, idir) = virial_xc(jdir, idir) + tmp
3143 2628 : virial_xc(idir, jdir) = virial_xc(jdir, idir)
3144 : END DO
3145 : END DO
3146 :
3147 292 : END SUBROUTINE virial_drho_drho
3148 :
3149 : ! **************************************************************************************************
3150 : !> \brief ...
3151 : !> \param rho_r ...
3152 : !> \param pw_pool ...
3153 : !> \param virial_xc ...
3154 : !> \param deriv_data ...
3155 : ! **************************************************************************************************
3156 132 : SUBROUTINE virial_laplace(rho_r, pw_pool, virial_xc, deriv_data)
3157 : TYPE(pw_r3d_rs_type), TARGET :: rho_r
3158 : TYPE(pw_pool_type), POINTER, INTENT(IN) :: pw_pool
3159 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT) :: virial_xc
3160 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: deriv_data
3161 :
3162 : CHARACTER(len=*), PARAMETER :: routineN = 'virial_laplace'
3163 :
3164 : INTEGER :: handle, idir, jdir
3165 : TYPE(pw_r3d_rs_type), POINTER :: virial_pw
3166 : TYPE(pw_c1d_gs_type), POINTER :: tmp_g, rho_g
3167 : INTEGER, DIMENSION(3) :: my_deriv
3168 :
3169 132 : CALL timeset(routineN, handle)
3170 :
3171 132 : NULLIFY (virial_pw, tmp_g, rho_g)
3172 132 : ALLOCATE (virial_pw, tmp_g, rho_g)
3173 132 : CALL pw_pool%create_pw(virial_pw)
3174 132 : CALL pw_pool%create_pw(tmp_g)
3175 132 : CALL pw_pool%create_pw(rho_g)
3176 132 : CALL pw_zero(virial_pw)
3177 132 : CALL pw_transfer(rho_r, rho_g)
3178 528 : DO idir = 1, 3
3179 1320 : DO jdir = idir, 3
3180 792 : CALL pw_copy(rho_g, tmp_g)
3181 :
3182 792 : my_deriv = 0
3183 792 : my_deriv(idir) = 1
3184 792 : my_deriv(jdir) = my_deriv(jdir) + 1
3185 :
3186 792 : CALL pw_derive(tmp_g, my_deriv)
3187 792 : CALL pw_transfer(tmp_g, virial_pw)
3188 : virial_xc(idir, jdir) = virial_xc(idir, jdir) - 2.0_dp*virial_pw%pw_grid%dvol* &
3189 : accurate_dot_product(virial_pw%array(:, :, :), &
3190 792 : deriv_data(:, :, :))
3191 1188 : virial_xc(jdir, idir) = virial_xc(idir, jdir)
3192 : END DO
3193 : END DO
3194 132 : CALL pw_pool%give_back_pw(virial_pw)
3195 132 : CALL pw_pool%give_back_pw(tmp_g)
3196 132 : CALL pw_pool%give_back_pw(rho_g)
3197 132 : DEALLOCATE (virial_pw, tmp_g, rho_g)
3198 :
3199 132 : CALL timestop(handle)
3200 :
3201 132 : END SUBROUTINE virial_laplace
3202 :
3203 : ! **************************************************************************************************
3204 : !> \brief Prepare objects for the calculation of the 2nd derivatives of the density functional.
3205 : !> The calculation must then be performed with xc_calc_2nd_deriv.
3206 : !> \param deriv_set object containing the XC derivatives (out)
3207 : !> \param rho_set object that will contain the density at which the
3208 : !> derivatives were calculated
3209 : !> \param rho_r the place where you evaluate the derivative
3210 : !> \param pw_pool the pool for the grids
3211 : !> \param xc_section which functional should be used and how to calculate it
3212 : !> \param tau_r kinetic energy density in real space
3213 : ! **************************************************************************************************
3214 11646 : SUBROUTINE xc_prep_2nd_deriv(deriv_set, &
3215 : rho_set, rho_r, pw_pool, xc_section, tau_r)
3216 :
3217 : TYPE(xc_derivative_set_type) :: deriv_set
3218 : TYPE(xc_rho_set_type) :: rho_set
3219 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
3220 : TYPE(pw_pool_type), POINTER :: pw_pool
3221 : TYPE(section_vals_type), POINTER :: xc_section
3222 : TYPE(pw_r3d_rs_type), DIMENSION(:), OPTIONAL, POINTER :: tau_r
3223 :
3224 : CHARACTER(len=*), PARAMETER :: routineN = 'xc_prep_2nd_deriv'
3225 :
3226 : INTEGER :: handle, nspins
3227 : LOGICAL :: lsd
3228 11646 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
3229 11646 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: tau
3230 :
3231 11646 : CALL timeset(routineN, handle)
3232 :
3233 11646 : CPASSERT(ASSOCIATED(xc_section))
3234 11646 : CPASSERT(ASSOCIATED(pw_pool))
3235 :
3236 11646 : nspins = SIZE(rho_r)
3237 11646 : lsd = (nspins /= 1)
3238 :
3239 11646 : NULLIFY (rho_g, tau)
3240 11646 : IF (PRESENT(tau_r)) &
3241 11218 : tau => tau_r
3242 :
3243 11646 : IF (section_get_lval(xc_section, "2ND_DERIV_ANALYTICAL")) THEN
3244 : CALL xc_rho_set_and_dset_create(rho_set, deriv_set, 2, &
3245 : rho_r, rho_g, tau, xc_section, pw_pool, &
3246 11562 : calc_potential=.TRUE.)
3247 : ELSE
3248 : CALL xc_rho_set_and_dset_create(rho_set, deriv_set, 1, &
3249 : rho_r, rho_g, tau, xc_section, pw_pool, &
3250 84 : calc_potential=.TRUE.)
3251 : END IF
3252 :
3253 11646 : CALL timestop(handle)
3254 :
3255 11646 : END SUBROUTINE xc_prep_2nd_deriv
3256 :
3257 : ! **************************************************************************************************
3258 : !> \brief divides derivatives from deriv_set by norm_drho
3259 : !> \param deriv_set ...
3260 : !> \param rho_set ...
3261 : !> \param lsd ...
3262 : ! **************************************************************************************************
3263 138635 : SUBROUTINE divide_by_norm_drho(deriv_set, rho_set, lsd)
3264 :
3265 : TYPE(xc_derivative_set_type), INTENT(INOUT) :: deriv_set
3266 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
3267 : LOGICAL, INTENT(IN) :: lsd
3268 :
3269 138635 : INTEGER, DIMENSION(:), POINTER :: split_desc
3270 : INTEGER :: idesc
3271 : INTEGER, DIMENSION(2, 3) :: bo
3272 : REAL(KIND=dp) :: drho_cutoff
3273 138635 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: norm_drho, norm_drhoa, norm_drhob
3274 1663620 : TYPE(cp_3d_r_cp_type), DIMENSION(3) :: drho, drhoa, drhob
3275 : TYPE(cp_sll_xc_deriv_type), POINTER :: pos
3276 : TYPE(xc_derivative_type), POINTER :: deriv_att
3277 :
3278 : ! check for unknown derivatives and divide by norm_drho where necessary
3279 :
3280 1386350 : bo = rho_set%local_bounds
3281 : CALL xc_rho_set_get(rho_set, drho_cutoff=drho_cutoff, norm_drho=norm_drho, &
3282 : norm_drhoa=norm_drhoa, norm_drhob=norm_drhob, &
3283 138635 : drho=drho, drhoa=drhoa, drhob=drhob, can_return_null=.TRUE.)
3284 :
3285 138635 : pos => deriv_set%derivs
3286 594946 : DO WHILE (cp_sll_xc_deriv_next(pos, el_att=deriv_att))
3287 456311 : CALL xc_derivative_get(deriv_att, split_desc=split_desc)
3288 973622 : DO idesc = 1, SIZE(split_desc)
3289 456311 : SELECT CASE (split_desc(idesc))
3290 : CASE (deriv_norm_drho)
3291 114495 : IF (ASSOCIATED(norm_drho)) THEN
3292 114495 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(deriv_att,norm_drho,drho_cutoff)
3293 : deriv_att%deriv_data(:, :, :) = deriv_att%deriv_data(:, :, :)/ &
3294 : MAX(norm_drho(:, :, :), drho_cutoff)
3295 : !$OMP END PARALLEL WORKSHARE
3296 0 : ELSE IF (ASSOCIATED(drho(1)%array)) THEN
3297 0 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(deriv_att,drho,drho_cutoff)
3298 : deriv_att%deriv_data(:, :, :) = deriv_att%deriv_data(:, :, :)/ &
3299 : MAX(SQRT(drho(1)%array(:, :, :)**2 + &
3300 : drho(2)%array(:, :, :)**2 + &
3301 : drho(3)%array(:, :, :)**2), drho_cutoff)
3302 : !$OMP END PARALLEL WORKSHARE
3303 0 : ELSE IF (ASSOCIATED(drhoa(1)%array) .AND. ASSOCIATED(drhob(1)%array)) THEN
3304 0 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(deriv_att,drhoa,drhob,drho_cutoff)
3305 : deriv_att%deriv_data(:, :, :) = deriv_att%deriv_data(:, :, :)/ &
3306 : MAX(SQRT((drhoa(1)%array(:, :, :) + drhob(1)%array(:, :, :))**2 + &
3307 : (drhoa(2)%array(:, :, :) + drhob(2)%array(:, :, :))**2 + &
3308 : (drhoa(3)%array(:, :, :) + drhob(3)%array(:, :, :))**2), drho_cutoff)
3309 : !$OMP END PARALLEL WORKSHARE
3310 : ELSE
3311 0 : CPABORT("Normalization of derivative requires any of norm_drho, drho or drhoa+drhob!")
3312 : END IF
3313 : CASE (deriv_norm_drhoa)
3314 19298 : IF (ASSOCIATED(norm_drhoa)) THEN
3315 19298 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(deriv_att,norm_drhoa,drho_cutoff)
3316 : deriv_att%deriv_data(:, :, :) = deriv_att%deriv_data(:, :, :)/ &
3317 : MAX(norm_drhoa(:, :, :), drho_cutoff)
3318 : !$OMP END PARALLEL WORKSHARE
3319 0 : ELSE IF (ASSOCIATED(drhoa(1)%array)) THEN
3320 0 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(deriv_att,drhoa,drho_cutoff)
3321 : deriv_att%deriv_data(:, :, :) = deriv_att%deriv_data(:, :, :)/ &
3322 : MAX(SQRT(drhoa(1)%array(:, :, :)**2 + &
3323 : drhoa(2)%array(:, :, :)**2 + &
3324 : drhoa(3)%array(:, :, :)**2), drho_cutoff)
3325 : !$OMP END PARALLEL WORKSHARE
3326 : ELSE
3327 0 : CPABORT("Normalization of derivative requires any of norm_drhoa or drhoa!")
3328 : END IF
3329 : CASE (deriv_norm_drhob)
3330 19294 : IF (ASSOCIATED(norm_drhob)) THEN
3331 19294 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(deriv_att,norm_drhob,drho_cutoff)
3332 : deriv_att%deriv_data(:, :, :) = deriv_att%deriv_data(:, :, :)/ &
3333 : MAX(norm_drhob(:, :, :), drho_cutoff)
3334 : !$OMP END PARALLEL WORKSHARE
3335 0 : ELSE IF (ASSOCIATED(drhob(1)%array)) THEN
3336 0 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(deriv_att,drhob,drho_cutoff)
3337 : deriv_att%deriv_data(:, :, :) = deriv_att%deriv_data(:, :, :)/ &
3338 : MAX(SQRT(drhob(1)%array(:, :, :)**2 + &
3339 : drhob(2)%array(:, :, :)**2 + &
3340 : drhob(3)%array(:, :, :)**2), drho_cutoff)
3341 : !$OMP END PARALLEL WORKSHARE
3342 : ELSE
3343 0 : CPABORT("Normalization of derivative requires any of norm_drhob or drhob!")
3344 : END IF
3345 : CASE (deriv_rho, deriv_tau, deriv_laplace_rho)
3346 158321 : IF (lsd) &
3347 0 : CPABORT(TRIM(id_to_desc(split_desc(idesc)))//" not handled in lsd!'")
3348 : CASE (deriv_rhoa, deriv_rhob, deriv_tau_a, deriv_tau_b, deriv_laplace_rhoa, deriv_laplace_rhob)
3349 : CASE default
3350 378676 : CPABORT("Unknown derivative id")
3351 : END SELECT
3352 : END DO
3353 : END DO
3354 :
3355 138635 : END SUBROUTINE divide_by_norm_drho
3356 :
3357 : ! **************************************************************************************************
3358 : !> \brief allocates and calculates drho from given spin densities drhoa, drhob
3359 : !> \param drho ...
3360 : !> \param drhoa ...
3361 : !> \param drhob ...
3362 : ! **************************************************************************************************
3363 17488 : SUBROUTINE calc_drho_from_ab(drho, drhoa, drhob)
3364 : TYPE(cp_3d_r_cp_type), DIMENSION(3), INTENT(OUT) :: drho
3365 : TYPE(cp_3d_r_cp_type), DIMENSION(3), INTENT(IN) :: drhoa, drhob
3366 :
3367 : CHARACTER(len=*), PARAMETER :: routineN = 'calc_drho_from_ab'
3368 :
3369 : INTEGER :: handle, idir
3370 :
3371 4372 : CALL timeset(routineN, handle)
3372 :
3373 17488 : DO idir = 1, 3
3374 13116 : NULLIFY (drho(idir)%array)
3375 : ALLOCATE (drho(idir)%array(LBOUND(drhoa(1)%array, 1):UBOUND(drhoa(1)%array, 1), &
3376 : LBOUND(drhoa(1)%array, 2):UBOUND(drhoa(1)%array, 2), &
3377 65580 : LBOUND(drhoa(1)%array, 3):UBOUND(drhoa(1)%array, 3)))
3378 17488 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(drho,drhoa,drhob,idir)
3379 : drho(idir)%array(:, :, :) = drhoa(idir)%array(:, :, :) + drhob(idir)%array(:, :, :)
3380 : !$OMP END PARALLEL WORKSHARE
3381 : END DO
3382 :
3383 4372 : CALL timestop(handle)
3384 :
3385 4372 : END SUBROUTINE
3386 :
3387 : ! **************************************************************************************************
3388 : !> \brief allocates and calculates dot products of two density gradients
3389 : !> \param dr1dr ...
3390 : !> \param drho ...
3391 : !> \param drho1 ...
3392 : ! **************************************************************************************************
3393 23266 : SUBROUTINE prepare_dr1dr(dr1dr, drho, drho1)
3394 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
3395 : INTENT(OUT) :: dr1dr
3396 : TYPE(cp_3d_r_cp_type), DIMENSION(3), INTENT(IN) :: drho, drho1
3397 :
3398 : CHARACTER(len=*), PARAMETER :: routineN = 'prepare_dr1dr'
3399 :
3400 : INTEGER :: handle, idir
3401 :
3402 23266 : CALL timeset(routineN, handle)
3403 :
3404 0 : ALLOCATE (dr1dr(LBOUND(drho(1)%array, 1):UBOUND(drho(1)%array, 1), &
3405 : LBOUND(drho(1)%array, 2):UBOUND(drho(1)%array, 2), &
3406 116330 : LBOUND(drho(1)%array, 3):UBOUND(drho(1)%array, 3)))
3407 :
3408 23266 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(dr1dr,drho,drho1)
3409 : dr1dr(:, :, :) = drho(1)%array(:, :, :)*drho1(1)%array(:, :, :)
3410 : !$OMP END PARALLEL WORKSHARE
3411 69798 : DO idir = 2, 3
3412 69798 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(dr1dr,drho,drho1,idir)
3413 : dr1dr(:, :, :) = dr1dr(:, :, :) + drho(idir)%array(:, :, :)*drho1(idir)%array(:, :, :)
3414 : !$OMP END PARALLEL WORKSHARE
3415 : END DO
3416 :
3417 23266 : CALL timestop(handle)
3418 :
3419 23266 : END SUBROUTINE prepare_dr1dr
3420 :
3421 : ! **************************************************************************************************
3422 : !> \brief allocates and calculates dot product of two densities for triplets
3423 : !> \param dr1dr ...
3424 : !> \param drhoa ...
3425 : !> \param drhob ...
3426 : !> \param drho1a ...
3427 : !> \param drho1b ...
3428 : !> \param fac ...
3429 : ! **************************************************************************************************
3430 870 : SUBROUTINE prepare_dr1dr_ab(dr1dr, drhoa, drhob, drho1a, drho1b, fac)
3431 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
3432 : INTENT(OUT) :: dr1dr
3433 : TYPE(cp_3d_r_cp_type), DIMENSION(3), INTENT(IN) :: drhoa, drhob, drho1a, drho1b
3434 : REAL(KIND=dp), INTENT(IN) :: fac
3435 :
3436 : CHARACTER(len=*), PARAMETER :: routineN = 'prepare_dr1dr_ab'
3437 :
3438 : INTEGER :: handle, idir
3439 :
3440 870 : CALL timeset(routineN, handle)
3441 :
3442 0 : ALLOCATE (dr1dr(LBOUND(drhoa(1)%array, 1):UBOUND(drhoa(1)%array, 1), &
3443 : LBOUND(drhoa(1)%array, 2):UBOUND(drhoa(1)%array, 2), &
3444 4350 : LBOUND(drhoa(1)%array, 3):UBOUND(drhoa(1)%array, 3)))
3445 :
3446 870 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(fac,dr1dr,drho1a,drho1b,drhoa,drhob)
3447 : dr1dr(:, :, :) = drhoa(1)%array(:, :, :)*(drho1a(1)%array(:, :, :) + &
3448 : fac*drho1b(1)%array(:, :, :)) + &
3449 : drhob(1)%array(:, :, :)*(fac*drho1a(1)%array(:, :, :) + &
3450 : drho1b(1)%array(:, :, :))
3451 : !$OMP END PARALLEL WORKSHARE
3452 2610 : DO idir = 2, 3
3453 2610 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(fac,dr1dr,drho1a,drho1b,drhoa,drhob,idir)
3454 : dr1dr(:, :, :) = dr1dr(:, :, :) + &
3455 : drhoa(idir)%array(:, :, :)*(drho1a(idir)%array(:, :, :) + &
3456 : fac*drho1b(idir)%array(:, :, :)) + &
3457 : drhob(idir)%array(:, :, :)*(fac*drho1a(idir)%array(:, :, :) + &
3458 : drho1b(idir)%array(:, :, :))
3459 : !$OMP END PARALLEL WORKSHARE
3460 : END DO
3461 :
3462 870 : CALL timestop(handle)
3463 :
3464 870 : END SUBROUTINE prepare_dr1dr_ab
3465 :
3466 : ! **************************************************************************************************
3467 : !> \brief checks for gradients
3468 : !> \param deriv_set ...
3469 : !> \param lsd ...
3470 : !> \param gradient_f ...
3471 : !> \param tau_f ...
3472 : !> \param laplace_f ...
3473 : ! **************************************************************************************************
3474 140200 : SUBROUTINE check_for_derivatives(deriv_set, lsd, rho_f, gradient_f, tau_f, laplace_f)
3475 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
3476 : LOGICAL, INTENT(IN) :: lsd
3477 : LOGICAL, INTENT(OUT) :: rho_f, gradient_f, tau_f, laplace_f
3478 :
3479 : CHARACTER(len=*), PARAMETER :: routineN = 'check_for_derivatives'
3480 :
3481 : INTEGER :: handle, iorder, order
3482 140200 : INTEGER, DIMENSION(:), POINTER :: split_desc
3483 : TYPE(cp_sll_xc_deriv_type), POINTER :: pos
3484 : TYPE(xc_derivative_type), POINTER :: deriv_att
3485 :
3486 140200 : CALL timeset(routineN, handle)
3487 :
3488 140200 : rho_f = .FALSE.
3489 140200 : gradient_f = .FALSE.
3490 140200 : tau_f = .FALSE.
3491 140200 : laplace_f = .FALSE.
3492 : ! check for unknown derivatives
3493 140200 : pos => deriv_set%derivs
3494 639140 : DO WHILE (cp_sll_xc_deriv_next(pos, el_att=deriv_att))
3495 : CALL xc_derivative_get(deriv_att, order=order, &
3496 498940 : split_desc=split_desc)
3497 639140 : IF (lsd) THEN
3498 305685 : DO iorder = 1, size(split_desc)
3499 152117 : SELECT CASE (split_desc(iorder))
3500 : CASE (deriv_rhoa, deriv_rhob)
3501 80026 : rho_f = .TRUE.
3502 : CASE (deriv_norm_drho, deriv_norm_drhoa, deriv_norm_drhob)
3503 70238 : gradient_f = .TRUE.
3504 : CASE (deriv_tau_a, deriv_tau_b)
3505 2376 : tau_f = .TRUE.
3506 : CASE (deriv_laplace_rhoa, deriv_laplace_rhob)
3507 928 : laplace_f = .TRUE.
3508 : CASE (deriv_rho, deriv_tau, deriv_laplace_rho)
3509 0 : CPABORT("Derivative not handled in lsd!")
3510 : CASE default
3511 153568 : CPABORT("Unknown derivative id")
3512 : END SELECT
3513 : END DO
3514 : ELSE
3515 644385 : DO iorder = 1, size(split_desc)
3516 346823 : SELECT CASE (split_desc(iorder))
3517 : CASE (deriv_rho)
3518 176005 : rho_f = .TRUE.
3519 : CASE (deriv_tau)
3520 4712 : tau_f = .TRUE.
3521 : CASE (deriv_norm_drho)
3522 115543 : gradient_f = .TRUE.
3523 : CASE (deriv_laplace_rho)
3524 1302 : laplace_f = .TRUE.
3525 : CASE default
3526 297562 : CPABORT("Unknown derivative id")
3527 : END SELECT
3528 : END DO
3529 : END IF
3530 : END DO
3531 :
3532 140200 : CALL timestop(handle)
3533 :
3534 140200 : END SUBROUTINE check_for_derivatives
3535 :
3536 : END MODULE xc
3537 :
|
