Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Exchange and Correlation functional calculations
10 : !> \par History
11 : !> (13-Feb-2001) JGH, based on earlier version of apsi
12 : !> 02.2003 Many many changes [fawzi]
13 : !> 03.2004 new xc interface [fawzi]
14 : !> 04.2004 kinetic functionals [fawzi]
15 : !> \author fawzi
16 : ! **************************************************************************************************
17 : MODULE xc
18 : #:include 'xc.fypp'
19 : USE cp_array_utils, ONLY: cp_3d_r_cp_type
20 : USE cp_linked_list_xc_deriv, ONLY: cp_sll_xc_deriv_next, &
21 : cp_sll_xc_deriv_type
22 : USE cp_log_handling, ONLY: cp_get_default_logger, &
23 : cp_logger_get_default_unit_nr, &
24 : cp_logger_type, &
25 : cp_to_string
26 : USE input_section_types, ONLY: section_get_ival, &
27 : section_get_lval, &
28 : section_get_rval, &
29 : section_vals_get_subs_vals, &
30 : section_vals_type, &
31 : section_vals_val_get
32 : USE kahan_sum, ONLY: accurate_dot_product, &
33 : accurate_sum
34 : USE kinds, ONLY: default_path_length, &
35 : dp
36 : USE pw_grid_types, ONLY: PW_MODE_DISTRIBUTED, &
37 : pw_grid_type
38 : USE pw_methods, ONLY: pw_axpy, &
39 : pw_copy, &
40 : pw_derive, &
41 : pw_scale, &
42 : pw_transfer, &
43 : pw_zero, pw_integrate_function
44 : USE pw_pool_types, ONLY: &
45 : pw_pool_type
46 : USE pw_types, ONLY: &
47 : pw_c1d_gs_type, pw_r3d_rs_type
48 : USE xc_derivative_desc, ONLY: &
49 : deriv_rho, deriv_rhoa, deriv_rhob, &
50 : deriv_norm_drhoa, deriv_norm_drhob, deriv_norm_drho, deriv_tau_a, deriv_tau_b, deriv_tau, &
51 : deriv_laplace_rho, deriv_laplace_rhoa, deriv_laplace_rhob, id_to_desc
52 : USE xc_derivative_set_types, ONLY: xc_derivative_set_type, &
53 : xc_dset_create, &
54 : xc_dset_get_derivative, &
55 : xc_dset_release, &
56 : xc_dset_zero_all, xc_dset_recover_pw
57 : USE xc_derivative_types, ONLY: xc_derivative_get, &
58 : xc_derivative_type
59 : USE xc_derivatives, ONLY: xc_functionals_eval, &
60 : xc_functionals_get_needs
61 : USE xc_input_constants, ONLY: &
62 : xc_debug_new_routine, xc_default_f_routine, xc_test_lsd_f_routine
63 : USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
64 : USE xc_rho_set_types, ONLY: xc_rho_set_create, &
65 : xc_rho_set_get, &
66 : xc_rho_set_release, &
67 : xc_rho_set_type, &
68 : xc_rho_set_update, xc_rho_set_recover_pw
69 : USE xc_util, ONLY: xc_pw_smooth, xc_pw_laplace, xc_pw_divergence, xc_requires_tmp_g
70 : #include "../base/base_uses.f90"
71 :
72 : IMPLICIT NONE
73 : PRIVATE
74 : PUBLIC :: xc_vxc_pw_create1, xc_vxc_pw_create, &
75 : xc_exc_calc, xc_calc_2nd_deriv_analytical, xc_calc_2nd_deriv_numerical, xc_calc_3rd_deriv_analytical, &
76 : xc_calc_2nd_deriv, xc_prep_2nd_deriv, xc_prep_3rd_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 7696 : 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 15392 : calc_potential=.FALSE.)
101 7696 : res = (needs%tau_spin .OR. needs%tau)
102 :
103 7696 : END FUNCTION xc_uses_kinetic_energy_density
104 :
105 : ! **************************************************************************************************
106 : !> \brief ...
107 : !> \param xc_fun_section ...
108 : !> \param lsd ...
109 : !> \return ...
110 : ! **************************************************************************************************
111 7712 : 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 15424 : calc_potential=.FALSE.)
121 7712 : res = (needs%norm_drho .OR. needs%norm_drho_spin)
122 :
123 7712 : END FUNCTION xc_uses_norm_drho
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 222558 : 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 111279 : CPASSERT(ASSOCIATED(rho_r))
164 111279 : CPASSERT(ASSOCIATED(xc_section))
165 111279 : CPASSERT(.NOT. ASSOCIATED(vxc_rho))
166 111279 : CPASSERT(.NOT. ASSOCIATED(vxc_tau))
167 :
168 111279 : virial_xc = 0.0_dp
169 :
170 : CALL section_vals_val_get(xc_section, "FUNCTIONAL_ROUTINE", &
171 111279 : i_val=f_routine)
172 111279 : 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 111279 : 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 111279 : xc_section=xc_section, pw_pool=pw_pool)
187 : CASE default
188 : END SELECT
189 :
190 111279 : 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 272950 : 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 136475 : CALL timeset(routineN, handle)
859 :
860 136475 : CPASSERT(ASSOCIATED(pw_pool))
861 :
862 136475 : nspins = SIZE(rho_r)
863 136475 : lsd = (nspins /= 1)
864 :
865 136475 : xc_fun_sections => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
866 :
867 : ! Create deriv_set object
868 136475 : CALL xc_dset_create(deriv_set, pw_pool)
869 :
870 : ! Create objects for density related stuff
871 : CALL xc_rho_set_create(rho_set, &
872 : rho_r(1)%pw_grid%bounds_local, &
873 : rho_cutoff=section_get_rval(xc_section, "density_cutoff"), &
874 : drho_cutoff=section_get_rval(xc_section, "gradient_cutoff"), &
875 136475 : tau_cutoff=section_get_rval(xc_section, "tau_cutoff"))
876 :
877 : ! Calculate density stuff, for example the gradient of rho, according to the functional needs
878 : CALL xc_rho_set_update(rho_set, rho_r, rho_g, tau, &
879 : xc_functionals_get_needs(xc_fun_sections, lsd, calc_potential), &
880 : section_get_ival(xc_section, "XC_GRID%XC_DERIV"), &
881 : section_get_ival(xc_section, "XC_GRID%XC_SMOOTH_RHO"), &
882 136475 : pw_pool)
883 :
884 : ! Calculate values of the functional on the grid
885 : CALL xc_functionals_eval(xc_fun_sections, &
886 : lsd=lsd, &
887 : rho_set=rho_set, &
888 : deriv_set=deriv_set, &
889 136475 : deriv_order=deriv_order)
890 :
891 136475 : CALL divide_by_norm_drho(deriv_set, rho_set, lsd)
892 :
893 136475 : CALL timestop(handle)
894 :
895 136475 : END SUBROUTINE xc_rho_set_and_dset_create
896 :
897 : ! **************************************************************************************************
898 : !> \brief smooths the cutoff on rho with a function smoothderiv_rho that is 0
899 : !> for rho<rho_cutoff and 1 for rho>rho_cutoff*rho_smooth_cutoff_range:
900 : !> E= integral e_0*smoothderiv_rho => dE/d...= de/d... * smooth,
901 : !> dE/drho = de/drho * smooth + e_0 * dsmooth/drho
902 : !> \param pot the potential to smooth
903 : !> \param rho , rhoa,rhob: the value of the density (used to apply the cutoff)
904 : !> \param rhoa ...
905 : !> \param rhob ...
906 : !> \param rho_cutoff the value at whch the cutoff function must go to 0
907 : !> \param rho_smooth_cutoff_range range of the smoothing
908 : !> \param e_0 value of e_0, if given it is assumed that pot is the derivative
909 : !> wrt. to rho, and needs the dsmooth*e_0 contribution
910 : !> \param e_0_scale_factor ...
911 : !> \author Fawzi Mohamed
912 : ! **************************************************************************************************
913 259537 : SUBROUTINE smooth_cutoff(pot, rho, rhoa, rhob, rho_cutoff, &
914 : rho_smooth_cutoff_range, e_0, e_0_scale_factor)
915 : REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN), &
916 : POINTER :: pot, rho, rhoa, rhob
917 : REAL(kind=dp), INTENT(in) :: rho_cutoff, rho_smooth_cutoff_range
918 : REAL(kind=dp), DIMENSION(:, :, :), OPTIONAL, &
919 : POINTER :: e_0
920 : REAL(kind=dp), INTENT(in), OPTIONAL :: e_0_scale_factor
921 :
922 : INTEGER :: i, j, k
923 : INTEGER, DIMENSION(2, 3) :: bo
924 : REAL(kind=dp) :: my_e_0_scale_factor, my_rho, my_rho_n, my_rho_n2, rho_smooth_cutoff, &
925 : rho_smooth_cutoff_2, rho_smooth_cutoff_range_2
926 :
927 259537 : CPASSERT(ASSOCIATED(pot))
928 1038148 : bo(1, :) = LBOUND(pot)
929 1038148 : bo(2, :) = UBOUND(pot)
930 259537 : my_e_0_scale_factor = 1.0_dp
931 259537 : IF (PRESENT(e_0_scale_factor)) my_e_0_scale_factor = e_0_scale_factor
932 259537 : rho_smooth_cutoff = rho_cutoff*rho_smooth_cutoff_range
933 259537 : rho_smooth_cutoff_2 = (rho_cutoff + rho_smooth_cutoff)/2
934 259537 : rho_smooth_cutoff_range_2 = rho_smooth_cutoff_2 - rho_cutoff
935 :
936 259537 : IF (rho_smooth_cutoff_range > 0.0_dp) THEN
937 2 : IF (PRESENT(e_0)) THEN
938 0 : CPASSERT(ASSOCIATED(e_0))
939 0 : IF (ASSOCIATED(rho)) THEN
940 : !$OMP PARALLEL DO DEFAULT(NONE) &
941 : !$OMP SHARED(bo,e_0,pot,rho,rho_cutoff,rho_smooth_cutoff,rho_smooth_cutoff_2, &
942 : !$OMP rho_smooth_cutoff_range_2,my_e_0_scale_factor) &
943 : !$OMP PRIVATE(k,j,i,my_rho,my_rho_n,my_rho_n2) &
944 0 : !$OMP COLLAPSE(3)
945 : DO k = bo(1, 3), bo(2, 3)
946 : DO j = bo(1, 2), bo(2, 2)
947 : DO i = bo(1, 1), bo(2, 1)
948 : my_rho = rho(i, j, k)
949 : IF (my_rho < rho_smooth_cutoff) THEN
950 : IF (my_rho < rho_cutoff) THEN
951 : pot(i, j, k) = 0.0_dp
952 : ELSEIF (my_rho < rho_smooth_cutoff_2) THEN
953 : my_rho_n = (my_rho - rho_cutoff)/rho_smooth_cutoff_range_2
954 : my_rho_n2 = my_rho_n*my_rho_n
955 : pot(i, j, k) = pot(i, j, k)* &
956 : my_rho_n2*(my_rho_n - 0.5_dp*my_rho_n2) + &
957 : my_e_0_scale_factor*e_0(i, j, k)* &
958 : my_rho_n2*(3.0_dp - 2.0_dp*my_rho_n) &
959 : /rho_smooth_cutoff_range_2
960 : ELSE
961 : my_rho_n = 2.0_dp - (my_rho - rho_cutoff)/rho_smooth_cutoff_range_2
962 : my_rho_n2 = my_rho_n*my_rho_n
963 : pot(i, j, k) = pot(i, j, k)* &
964 : (1.0_dp - my_rho_n2*(my_rho_n - 0.5_dp*my_rho_n2)) &
965 : + my_e_0_scale_factor*e_0(i, j, k)* &
966 : my_rho_n2*(3.0_dp - 2.0_dp*my_rho_n) &
967 : /rho_smooth_cutoff_range_2
968 : END IF
969 : END IF
970 : END DO
971 : END DO
972 : END DO
973 : !$OMP END PARALLEL DO
974 : ELSE
975 : !$OMP PARALLEL DO DEFAULT(NONE) &
976 : !$OMP SHARED(bo,pot,e_0,rhoa,rhob,rho_cutoff,rho_smooth_cutoff,rho_smooth_cutoff_2, &
977 : !$OMP rho_smooth_cutoff_range_2,my_e_0_scale_factor) &
978 : !$OMP PRIVATE(k,j,i,my_rho,my_rho_n,my_rho_n2) &
979 0 : !$OMP COLLAPSE(3)
980 : DO k = bo(1, 3), bo(2, 3)
981 : DO j = bo(1, 2), bo(2, 2)
982 : DO i = bo(1, 1), bo(2, 1)
983 : my_rho = rhoa(i, j, k) + rhob(i, j, k)
984 : IF (my_rho < rho_smooth_cutoff) THEN
985 : IF (my_rho < rho_cutoff) THEN
986 : pot(i, j, k) = 0.0_dp
987 : ELSEIF (my_rho < rho_smooth_cutoff_2) THEN
988 : my_rho_n = (my_rho - rho_cutoff)/rho_smooth_cutoff_range_2
989 : my_rho_n2 = my_rho_n*my_rho_n
990 : pot(i, j, k) = pot(i, j, k)* &
991 : my_rho_n2*(my_rho_n - 0.5_dp*my_rho_n2) + &
992 : my_e_0_scale_factor*e_0(i, j, k)* &
993 : my_rho_n2*(3.0_dp - 2.0_dp*my_rho_n) &
994 : /rho_smooth_cutoff_range_2
995 : ELSE
996 : my_rho_n = 2.0_dp - (my_rho - rho_cutoff)/rho_smooth_cutoff_range_2
997 : my_rho_n2 = my_rho_n*my_rho_n
998 : pot(i, j, k) = pot(i, j, k)* &
999 : (1.0_dp - my_rho_n2*(my_rho_n - 0.5_dp*my_rho_n2)) &
1000 : + my_e_0_scale_factor*e_0(i, j, k)* &
1001 : my_rho_n2*(3.0_dp - 2.0_dp*my_rho_n) &
1002 : /rho_smooth_cutoff_range_2
1003 : END IF
1004 : END IF
1005 : END DO
1006 : END DO
1007 : END DO
1008 : !$OMP END PARALLEL DO
1009 : END IF
1010 : ELSE
1011 2 : IF (ASSOCIATED(rho)) THEN
1012 : !$OMP PARALLEL DO DEFAULT(NONE) &
1013 : !$OMP SHARED(bo,pot,rho_cutoff,rho_smooth_cutoff,rho_smooth_cutoff_2, &
1014 : !$OMP rho_smooth_cutoff_range_2,rho) &
1015 : !$OMP PRIVATE(k,j,i,my_rho,my_rho_n,my_rho_n2) &
1016 2 : !$OMP COLLAPSE(3)
1017 : DO k = bo(1, 3), bo(2, 3)
1018 : DO j = bo(1, 2), bo(2, 2)
1019 : DO i = bo(1, 1), bo(2, 1)
1020 : my_rho = rho(i, j, k)
1021 : IF (my_rho < rho_smooth_cutoff) THEN
1022 : IF (my_rho < rho_cutoff) THEN
1023 : pot(i, j, k) = 0.0_dp
1024 : ELSEIF (my_rho < rho_smooth_cutoff_2) THEN
1025 : my_rho_n = (my_rho - rho_cutoff)/rho_smooth_cutoff_range_2
1026 : my_rho_n2 = my_rho_n*my_rho_n
1027 : pot(i, j, k) = pot(i, j, k)* &
1028 : my_rho_n2*(my_rho_n - 0.5_dp*my_rho_n2)
1029 : ELSE
1030 : my_rho_n = 2.0_dp - (my_rho - rho_cutoff)/rho_smooth_cutoff_range_2
1031 : my_rho_n2 = my_rho_n*my_rho_n
1032 : pot(i, j, k) = pot(i, j, k)* &
1033 : (1.0_dp - my_rho_n2*(my_rho_n - 0.5_dp*my_rho_n2))
1034 : END IF
1035 : END IF
1036 : END DO
1037 : END DO
1038 : END DO
1039 : !$OMP END PARALLEL DO
1040 : ELSE
1041 : !$OMP PARALLEL DO DEFAULT(NONE) &
1042 : !$OMP SHARED(bo,pot,rho_cutoff,rho_smooth_cutoff,rho_smooth_cutoff_2, &
1043 : !$OMP rho_smooth_cutoff_range_2,rhoa,rhob) &
1044 : !$OMP PRIVATE(k,j,i,my_rho,my_rho_n,my_rho_n2) &
1045 0 : !$OMP COLLAPSE(3)
1046 : DO k = bo(1, 3), bo(2, 3)
1047 : DO j = bo(1, 2), bo(2, 2)
1048 : DO i = bo(1, 1), bo(2, 1)
1049 : my_rho = rhoa(i, j, k) + rhob(i, j, k)
1050 : IF (my_rho < rho_smooth_cutoff) THEN
1051 : IF (my_rho < rho_cutoff) THEN
1052 : pot(i, j, k) = 0.0_dp
1053 : ELSEIF (my_rho < rho_smooth_cutoff_2) THEN
1054 : my_rho_n = (my_rho - rho_cutoff)/rho_smooth_cutoff_range_2
1055 : my_rho_n2 = my_rho_n*my_rho_n
1056 : pot(i, j, k) = pot(i, j, k)* &
1057 : my_rho_n2*(my_rho_n - 0.5_dp*my_rho_n2)
1058 : ELSE
1059 : my_rho_n = 2.0_dp - (my_rho - rho_cutoff)/rho_smooth_cutoff_range_2
1060 : my_rho_n2 = my_rho_n*my_rho_n
1061 : pot(i, j, k) = pot(i, j, k)* &
1062 : (1.0_dp - my_rho_n2*(my_rho_n - 0.5_dp*my_rho_n2))
1063 : END IF
1064 : END IF
1065 : END DO
1066 : END DO
1067 : END DO
1068 : !$OMP END PARALLEL DO
1069 : END IF
1070 : END IF
1071 : END IF
1072 259537 : END SUBROUTINE smooth_cutoff
1073 :
1074 64 : SUBROUTINE calc_xc_density(pot, rho, rho_cutoff)
1075 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: pot
1076 : TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(INOUT) :: rho
1077 : REAL(kind=dp), INTENT(in) :: rho_cutoff
1078 :
1079 : INTEGER :: i, j, k, nspins
1080 : INTEGER, DIMENSION(2, 3) :: bo
1081 : REAL(kind=dp) :: eps1, eps2, my_rho, my_pot
1082 :
1083 256 : bo(1, :) = LBOUND(pot%array)
1084 256 : bo(2, :) = UBOUND(pot%array)
1085 64 : nspins = SIZE(rho)
1086 :
1087 64 : eps1 = rho_cutoff*1.E-4_dp
1088 64 : eps2 = rho_cutoff
1089 :
1090 3160 : DO k = bo(1, 3), bo(2, 3)
1091 161272 : DO j = bo(1, 2), bo(2, 2)
1092 4529376 : DO i = bo(1, 1), bo(2, 1)
1093 4368168 : my_pot = pot%array(i, j, k)
1094 4368168 : IF (nspins == 2) THEN
1095 339714 : my_rho = rho(1)%array(i, j, k) + rho(2)%array(i, j, k)
1096 : ELSE
1097 4028454 : my_rho = rho(1)%array(i, j, k)
1098 : END IF
1099 4526280 : IF (my_rho > eps1) THEN
1100 4139578 : pot%array(i, j, k) = my_pot/my_rho
1101 228590 : ELSE IF (my_rho < eps2) THEN
1102 228590 : pot%array(i, j, k) = 0.0_dp
1103 : ELSE
1104 0 : pot%array(i, j, k) = MIN(my_pot/my_rho, my_rho**(1._dp/3._dp))
1105 : END IF
1106 : END DO
1107 : END DO
1108 : END DO
1109 :
1110 64 : END SUBROUTINE calc_xc_density
1111 :
1112 : ! **************************************************************************************************
1113 : !> \brief Exchange and Correlation functional calculations
1114 : !> \param vxc_rho will contain the v_xc part that depend on rho
1115 : !> (if one of the chosen xc functionals has it it is allocated and you
1116 : !> are responsible for it)
1117 : !> \param vxc_tau will contain the kinetic tau part of v_xc
1118 : !> (if one of the chosen xc functionals has it it is allocated and you
1119 : !> are responsible for it)
1120 : !> \param exc the xc energy
1121 : !> \param rho_r the value of the density in the real space
1122 : !> \param rho_g value of the density in the g space (needs to be associated
1123 : !> only for gradient corrections)
1124 : !> \param tau value of the kinetic density tau on the grid (can be null,
1125 : !> used only with meta functionals)
1126 : !> \param xc_section which functional to calculate, and how to do it
1127 : !> \param pw_pool the pool for the grids
1128 : !> \param compute_virial ...
1129 : !> \param virial_xc ...
1130 : !> \param exc_r the value of the xc functional in the real space
1131 : !> \par History
1132 : !> JGH (13-Jun-2002): adaptation to new functionals
1133 : !> Fawzi (11.2002): drho_g(1:3)->drho_g
1134 : !> Fawzi (1.2003). lsd version
1135 : !> Fawzi (11.2003): version using the new xc interface
1136 : !> Fawzi (03.2004): fft free for smoothed density and derivs, gga lsd
1137 : !> Fawzi (04.2004): metafunctionals
1138 : !> mguidon (12.2008) : laplace functionals
1139 : !> \author fawzi; based LDA version of JGH, based on earlier version of apsi
1140 : !> \note
1141 : !> Beware: some really dirty pointer handling!
1142 : !> energy should be kept consistent with xc_exc_calc
1143 : ! **************************************************************************************************
1144 113247 : SUBROUTINE xc_vxc_pw_create(vxc_rho, vxc_tau, exc, rho_r, rho_g, tau, xc_section, &
1145 : pw_pool, compute_virial, virial_xc, exc_r)
1146 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: vxc_rho, vxc_tau
1147 : REAL(KIND=dp), INTENT(out) :: exc
1148 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, tau
1149 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
1150 : TYPE(section_vals_type), POINTER :: xc_section
1151 : TYPE(pw_pool_type), POINTER :: pw_pool
1152 : LOGICAL :: compute_virial
1153 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT) :: virial_xc
1154 : TYPE(pw_r3d_rs_type), INTENT(INOUT), OPTIONAL :: exc_r
1155 :
1156 : CHARACTER(len=*), PARAMETER :: routineN = 'xc_vxc_pw_create'
1157 : INTEGER, DIMENSION(2), PARAMETER :: norm_drho_spin_name = [deriv_norm_drhoa, deriv_norm_drhob]
1158 :
1159 : INTEGER :: handle, idir, ispin, jdir, &
1160 : npoints, nspins, &
1161 : xc_deriv_method_id, xc_rho_smooth_id, deriv_id
1162 : INTEGER, DIMENSION(2, 3) :: bo
1163 : LOGICAL :: dealloc_pw_to_deriv, has_laplace, &
1164 : has_tau, lsd, use_virial, has_gradient, &
1165 : has_derivs, has_rho, dealloc_pw_to_deriv_rho
1166 : REAL(KIND=dp) :: density_smooth_cut_range, drho_cutoff, &
1167 : rho_cutoff
1168 113247 : REAL(kind=dp), DIMENSION(:, :, :), POINTER :: deriv_data, norm_drho, norm_drho_spin, &
1169 226494 : rho, rhoa, rhob
1170 : TYPE(cp_sll_xc_deriv_type), POINTER :: pos
1171 : TYPE(pw_grid_type), POINTER :: pw_grid
1172 792729 : TYPE(pw_r3d_rs_type), DIMENSION(3) :: pw_to_deriv, pw_to_deriv_rho
1173 : TYPE(pw_c1d_gs_type) :: tmp_g, vxc_g
1174 : TYPE(pw_r3d_rs_type) :: v_drho_r, virial_pw
1175 : TYPE(xc_derivative_set_type) :: deriv_set
1176 : TYPE(xc_derivative_type), POINTER :: deriv_att
1177 : TYPE(xc_rho_set_type) :: rho_set
1178 :
1179 113247 : CALL timeset(routineN, handle)
1180 113247 : NULLIFY (norm_drho_spin, norm_drho, pos)
1181 :
1182 113247 : pw_grid => rho_r(1)%pw_grid
1183 :
1184 113247 : CPASSERT(ASSOCIATED(xc_section))
1185 113247 : CPASSERT(ASSOCIATED(pw_pool))
1186 113247 : CPASSERT(.NOT. ASSOCIATED(vxc_rho))
1187 113247 : CPASSERT(.NOT. ASSOCIATED(vxc_tau))
1188 113247 : nspins = SIZE(rho_r)
1189 113247 : lsd = (nspins /= 1)
1190 113247 : IF (lsd) THEN
1191 22791 : CPASSERT(nspins == 2)
1192 : END IF
1193 :
1194 113247 : use_virial = compute_virial
1195 113247 : virial_xc = 0.0_dp
1196 :
1197 1132470 : bo = rho_r(1)%pw_grid%bounds_local
1198 113247 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
1199 :
1200 : ! calculate the potential derivatives
1201 : CALL xc_rho_set_and_dset_create(rho_set=rho_set, deriv_set=deriv_set, &
1202 : deriv_order=1, rho_r=rho_r, rho_g=rho_g, tau=tau, &
1203 : xc_section=xc_section, &
1204 : pw_pool=pw_pool, &
1205 113247 : calc_potential=.TRUE.)
1206 :
1207 : CALL section_vals_val_get(xc_section, "XC_GRID%XC_DERIV", &
1208 113247 : i_val=xc_deriv_method_id)
1209 : CALL section_vals_val_get(xc_section, "XC_GRID%XC_SMOOTH_RHO", &
1210 113247 : i_val=xc_rho_smooth_id)
1211 : CALL section_vals_val_get(xc_section, "DENSITY_SMOOTH_CUTOFF_RANGE", &
1212 113247 : r_val=density_smooth_cut_range)
1213 :
1214 : CALL xc_rho_set_get(rho_set, rho_cutoff=rho_cutoff, &
1215 113247 : drho_cutoff=drho_cutoff)
1216 :
1217 113247 : CALL check_for_derivatives(deriv_set, lsd, has_rho, has_gradient, has_tau, has_laplace)
1218 : ! check for unknown derivatives
1219 113247 : has_derivs = has_rho .OR. has_gradient .OR. has_tau .OR. has_laplace
1220 :
1221 475779 : ALLOCATE (vxc_rho(nspins))
1222 :
1223 : CALL xc_rho_set_get(rho_set, rho=rho, rhoa=rhoa, rhob=rhob, &
1224 113247 : can_return_null=.TRUE.)
1225 :
1226 : ! recover the vxc arrays
1227 113247 : IF (lsd) THEN
1228 22791 : CALL xc_dset_recover_pw(deriv_set, [deriv_rhoa], vxc_rho(1), pw_grid, pw_pool)
1229 22791 : CALL xc_dset_recover_pw(deriv_set, [deriv_rhob], vxc_rho(2), pw_grid, pw_pool)
1230 :
1231 : ELSE
1232 90456 : CALL xc_dset_recover_pw(deriv_set, [deriv_rho], vxc_rho(1), pw_grid, pw_pool)
1233 : END IF
1234 :
1235 113247 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_norm_drho])
1236 113247 : IF (ASSOCIATED(deriv_att)) THEN
1237 61363 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
1238 :
1239 : CALL xc_rho_set_get(rho_set, norm_drho=norm_drho, &
1240 : rho_cutoff=rho_cutoff, &
1241 : drho_cutoff=drho_cutoff, &
1242 61363 : can_return_null=.TRUE.)
1243 61363 : CALL xc_rho_set_recover_pw(rho_set, pw_grid, pw_pool, dealloc_pw_to_deriv_rho, drho=pw_to_deriv_rho)
1244 :
1245 61363 : CPASSERT(ASSOCIATED(deriv_data))
1246 61363 : IF (use_virial) THEN
1247 1568 : CALL pw_pool%create_pw(virial_pw)
1248 1568 : CALL pw_zero(virial_pw)
1249 6272 : DO idir = 1, 3
1250 4704 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(virial_pw,pw_to_deriv_rho,deriv_data,idir)
1251 : virial_pw%array(:, :, :) = pw_to_deriv_rho(idir)%array(:, :, :)*deriv_data(:, :, :)
1252 : !$OMP END PARALLEL WORKSHARE
1253 15680 : DO jdir = 1, idir
1254 : virial_xc(idir, jdir) = -pw_grid%dvol* &
1255 : accurate_dot_product(virial_pw%array(:, :, :), &
1256 9408 : pw_to_deriv_rho(jdir)%array(:, :, :))
1257 14112 : virial_xc(jdir, idir) = virial_xc(idir, jdir)
1258 : END DO
1259 : END DO
1260 1568 : CALL pw_pool%give_back_pw(virial_pw)
1261 : END IF ! use_virial
1262 245452 : DO idir = 1, 3
1263 184089 : CPASSERT(ASSOCIATED(pw_to_deriv_rho(idir)%array))
1264 245452 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(deriv_data,pw_to_deriv_rho,idir)
1265 : pw_to_deriv_rho(idir)%array(:, :, :) = pw_to_deriv_rho(idir)%array(:, :, :)*deriv_data(:, :, :)
1266 : !$OMP END PARALLEL WORKSHARE
1267 : END DO
1268 :
1269 : ! Deallocate pw to save memory
1270 61363 : CALL pw_pool%give_back_cr3d(deriv_att%deriv_data)
1271 :
1272 : END IF
1273 :
1274 113247 : IF ((has_gradient .AND. xc_requires_tmp_g(xc_deriv_method_id)) .OR. pw_grid%spherical) THEN
1275 61133 : CALL pw_pool%create_pw(vxc_g)
1276 61133 : IF (.NOT. pw_grid%spherical) THEN
1277 61133 : CALL pw_pool%create_pw(tmp_g)
1278 : END IF
1279 : END IF
1280 :
1281 249285 : DO ispin = 1, nspins
1282 :
1283 136038 : IF (lsd) THEN
1284 45582 : IF (ispin == 1) THEN
1285 : CALL xc_rho_set_get(rho_set, norm_drhoa=norm_drho_spin, &
1286 22791 : can_return_null=.TRUE.)
1287 22791 : IF (ASSOCIATED(norm_drho_spin)) CALL xc_rho_set_recover_pw( &
1288 14598 : rho_set, pw_grid, pw_pool, dealloc_pw_to_deriv, drhoa=pw_to_deriv)
1289 : ELSE
1290 : CALL xc_rho_set_get(rho_set, norm_drhob=norm_drho_spin, &
1291 22791 : can_return_null=.TRUE.)
1292 22791 : IF (ASSOCIATED(norm_drho_spin)) CALL xc_rho_set_recover_pw( &
1293 14598 : rho_set, pw_grid, pw_pool, dealloc_pw_to_deriv, drhob=pw_to_deriv)
1294 : END IF
1295 :
1296 91164 : deriv_att => xc_dset_get_derivative(deriv_set, [norm_drho_spin_name(ispin)])
1297 45582 : IF (ASSOCIATED(deriv_att)) THEN
1298 : CPASSERT(lsd)
1299 29196 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
1300 :
1301 29196 : IF (use_virial) THEN
1302 104 : CALL pw_pool%create_pw(virial_pw)
1303 104 : CALL pw_zero(virial_pw)
1304 416 : DO idir = 1, 3
1305 312 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(deriv_data,pw_to_deriv,virial_pw,idir)
1306 : virial_pw%array(:, :, :) = pw_to_deriv(idir)%array(:, :, :)*deriv_data(:, :, :)
1307 : !$OMP END PARALLEL WORKSHARE
1308 1040 : DO jdir = 1, idir
1309 : virial_xc(idir, jdir) = virial_xc(idir, jdir) - pw_grid%dvol* &
1310 : accurate_dot_product(virial_pw%array(:, :, :), &
1311 624 : pw_to_deriv(jdir)%array(:, :, :))
1312 936 : virial_xc(jdir, idir) = virial_xc(idir, jdir)
1313 : END DO
1314 : END DO
1315 104 : CALL pw_pool%give_back_pw(virial_pw)
1316 : END IF ! use_virial
1317 :
1318 116784 : DO idir = 1, 3
1319 116784 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(deriv_data,idir,pw_to_deriv)
1320 : pw_to_deriv(idir)%array(:, :, :) = deriv_data(:, :, :)*pw_to_deriv(idir)%array(:, :, :)
1321 : !$OMP END PARALLEL WORKSHARE
1322 : END DO
1323 : END IF ! deriv_att
1324 :
1325 : END IF ! LSD
1326 :
1327 136038 : IF (ASSOCIATED(pw_to_deriv_rho(1)%array)) THEN
1328 75203 : IF (.NOT. ASSOCIATED(pw_to_deriv(1)%array)) THEN
1329 47523 : pw_to_deriv = pw_to_deriv_rho
1330 47523 : dealloc_pw_to_deriv = ((.NOT. lsd) .OR. (ispin == 2))
1331 47523 : dealloc_pw_to_deriv = dealloc_pw_to_deriv .AND. dealloc_pw_to_deriv_rho
1332 : ELSE
1333 : ! This branch is called in case of open-shell systems
1334 : ! Add the contributions from norm_drho and norm_drho_spin
1335 110720 : DO idir = 1, 3
1336 83040 : CALL pw_axpy(pw_to_deriv_rho(idir), pw_to_deriv(idir))
1337 110720 : IF (ispin == 2) THEN
1338 41520 : IF (dealloc_pw_to_deriv_rho) THEN
1339 41520 : CALL pw_pool%give_back_pw(pw_to_deriv_rho(idir))
1340 : END IF
1341 : END IF
1342 : END DO
1343 : END IF
1344 : END IF
1345 :
1346 136038 : IF (ASSOCIATED(pw_to_deriv(1)%array)) THEN
1347 306876 : DO idir = 1, 3
1348 306876 : CALL pw_scale(pw_to_deriv(idir), -1.0_dp)
1349 : END DO
1350 :
1351 76719 : CALL xc_pw_divergence(xc_deriv_method_id, pw_to_deriv, tmp_g, vxc_g, vxc_rho(ispin))
1352 :
1353 76719 : IF (dealloc_pw_to_deriv) THEN
1354 306876 : DO idir = 1, 3
1355 306876 : CALL pw_pool%give_back_pw(pw_to_deriv(idir))
1356 : END DO
1357 : END IF
1358 : END IF
1359 :
1360 : ! Add laplace part to vxc_rho
1361 136038 : IF (has_laplace) THEN
1362 1092 : IF (lsd) THEN
1363 472 : IF (ispin == 1) THEN
1364 : deriv_id = deriv_laplace_rhoa
1365 : ELSE
1366 236 : deriv_id = deriv_laplace_rhob
1367 : END IF
1368 : ELSE
1369 : deriv_id = deriv_laplace_rho
1370 : END IF
1371 :
1372 2184 : CALL xc_dset_recover_pw(deriv_set, [deriv_id], pw_to_deriv(1), pw_grid)
1373 :
1374 1092 : IF (use_virial) CALL virial_laplace(rho_r(ispin), pw_pool, virial_xc, pw_to_deriv(1)%array)
1375 :
1376 1092 : CALL xc_pw_laplace(pw_to_deriv(1), pw_pool, xc_deriv_method_id)
1377 :
1378 1092 : CALL pw_axpy(pw_to_deriv(1), vxc_rho(ispin))
1379 :
1380 1092 : CALL pw_pool%give_back_pw(pw_to_deriv(1))
1381 : END IF
1382 :
1383 136038 : IF (pw_grid%spherical) THEN
1384 : ! filter vxc
1385 0 : CALL pw_transfer(vxc_rho(ispin), vxc_g)
1386 0 : CALL pw_transfer(vxc_g, vxc_rho(ispin))
1387 : END IF
1388 : CALL smooth_cutoff(pot=vxc_rho(ispin)%array, rho=rho, rhoa=rhoa, rhob=rhob, &
1389 : rho_cutoff=rho_cutoff*density_smooth_cut_range, &
1390 136038 : rho_smooth_cutoff_range=density_smooth_cut_range)
1391 :
1392 136038 : v_drho_r = vxc_rho(ispin)
1393 136038 : CALL pw_pool%create_pw(vxc_rho(ispin))
1394 136038 : CALL xc_pw_smooth(v_drho_r, vxc_rho(ispin), xc_rho_smooth_id)
1395 249285 : CALL pw_pool%give_back_pw(v_drho_r)
1396 : END DO
1397 :
1398 113247 : CALL pw_pool%give_back_pw(vxc_g)
1399 113247 : CALL pw_pool%give_back_pw(tmp_g)
1400 :
1401 : ! 0-deriv -> value of exc
1402 : ! this has to be kept consistent with xc_exc_calc
1403 113247 : IF (has_derivs) THEN
1404 112929 : CALL xc_dset_recover_pw(deriv_set, [INTEGER::], v_drho_r, pw_grid)
1405 :
1406 : CALL smooth_cutoff(pot=v_drho_r%array, rho=rho, rhoa=rhoa, rhob=rhob, &
1407 : rho_cutoff=rho_cutoff, &
1408 112929 : rho_smooth_cutoff_range=density_smooth_cut_range)
1409 :
1410 112929 : exc = pw_integrate_function(v_drho_r)
1411 : !
1412 : ! return the xc functional value at the grid points
1413 : !
1414 112929 : IF (PRESENT(exc_r)) THEN
1415 96 : exc_r = v_drho_r
1416 : ELSE
1417 112833 : CALL v_drho_r%release()
1418 : END IF
1419 : ELSE
1420 318 : exc = 0.0_dp
1421 : END IF
1422 :
1423 113247 : CALL xc_rho_set_release(rho_set, pw_pool=pw_pool)
1424 :
1425 : ! tau part
1426 113247 : IF (has_tau) THEN
1427 8374 : ALLOCATE (vxc_tau(nspins))
1428 2580 : IF (lsd) THEN
1429 634 : CALL xc_dset_recover_pw(deriv_set, [deriv_tau_a], vxc_tau(1), pw_grid)
1430 634 : CALL xc_dset_recover_pw(deriv_set, [deriv_tau_b], vxc_tau(2), pw_grid)
1431 : ELSE
1432 1946 : CALL xc_dset_recover_pw(deriv_set, [deriv_tau], vxc_tau(1), pw_grid)
1433 : END IF
1434 5794 : DO ispin = 1, nspins
1435 5794 : CPASSERT(ASSOCIATED(vxc_tau(ispin)%array))
1436 : END DO
1437 : END IF
1438 113247 : CALL xc_dset_release(deriv_set)
1439 :
1440 113247 : CALL timestop(handle)
1441 :
1442 2378187 : END SUBROUTINE xc_vxc_pw_create
1443 :
1444 : ! **************************************************************************************************
1445 : !> \brief calculates just the exchange and correlation energy
1446 : !> (no vxc)
1447 : !> \param rho_r realspace density on the grid
1448 : !> \param rho_g g-space density on the grid
1449 : !> \param tau kinetic energy density on the grid
1450 : !> \param xc_section XC parameters
1451 : !> \param pw_pool pool of plain-wave grids
1452 : !> \return the XC energy
1453 : !> \par History
1454 : !> 11.2003 created [fawzi]
1455 : !> \author fawzi
1456 : !> \note
1457 : !> has to be kept consistent with xc_vxc_pw_create
1458 : ! **************************************************************************************************
1459 21136 : FUNCTION xc_exc_calc(rho_r, rho_g, tau, xc_section, pw_pool) &
1460 : RESULT(exc)
1461 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, tau
1462 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
1463 : TYPE(section_vals_type), POINTER :: xc_section
1464 : TYPE(pw_pool_type), POINTER :: pw_pool
1465 : REAL(kind=dp) :: exc
1466 :
1467 : CHARACTER(len=*), PARAMETER :: routineN = 'xc_exc_calc'
1468 :
1469 : INTEGER :: handle
1470 : REAL(dp) :: density_smooth_cut_range, rho_cutoff
1471 10568 : REAL(dp), DIMENSION(:, :, :), POINTER :: e_0
1472 : TYPE(xc_derivative_set_type) :: deriv_set
1473 : TYPE(xc_derivative_type), POINTER :: deriv
1474 : TYPE(xc_rho_set_type) :: rho_set
1475 :
1476 10568 : CALL timeset(routineN, handle)
1477 10568 : NULLIFY (deriv, e_0)
1478 10568 : exc = 0.0_dp
1479 :
1480 : ! this has to be consistent with what is done in xc_vxc_pw_create
1481 : CALL xc_rho_set_and_dset_create(rho_set=rho_set, &
1482 : deriv_set=deriv_set, deriv_order=0, &
1483 : rho_r=rho_r, rho_g=rho_g, tau=tau, xc_section=xc_section, &
1484 : pw_pool=pw_pool, &
1485 10568 : calc_potential=.FALSE.)
1486 10568 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::])
1487 :
1488 10568 : IF (ASSOCIATED(deriv)) THEN
1489 10568 : CALL xc_derivative_get(deriv, deriv_data=e_0)
1490 :
1491 : CALL section_vals_val_get(xc_section, "DENSITY_CUTOFF", &
1492 10568 : r_val=rho_cutoff)
1493 : CALL section_vals_val_get(xc_section, "DENSITY_SMOOTH_CUTOFF_RANGE", &
1494 10568 : r_val=density_smooth_cut_range)
1495 : CALL smooth_cutoff(pot=e_0, rho=rho_set%rho, &
1496 : rhoa=rho_set%rhoa, rhob=rho_set%rhob, &
1497 : rho_cutoff=rho_cutoff, &
1498 10568 : rho_smooth_cutoff_range=density_smooth_cut_range)
1499 :
1500 10568 : exc = accurate_sum(e_0)*rho_r(1)%pw_grid%dvol
1501 10568 : IF (rho_r(1)%pw_grid%para%mode == PW_MODE_DISTRIBUTED) THEN
1502 10430 : CALL rho_r(1)%pw_grid%para%group%sum(exc)
1503 : END IF
1504 :
1505 10568 : CALL xc_rho_set_release(rho_set, pw_pool=pw_pool)
1506 10568 : CALL xc_dset_release(deriv_set)
1507 : END IF
1508 10568 : CALL timestop(handle)
1509 200792 : END FUNCTION xc_exc_calc
1510 :
1511 : ! **************************************************************************************************
1512 : !> \brief Caller routine to calculate the second order potential in the direction of rho1_r
1513 : !> \param v_xc XC potential, will be allocated, to be integrated with the KS density
1514 : !> \param v_xc_tau ...
1515 : !> \param deriv_set XC derivatives from xc_prep_2nd_deriv
1516 : !> \param rho_set XC rho set from KS rho from xc_prep_2nd_deriv
1517 : !> \param rho1_r first-order density in r space
1518 : !> \param rho1_g first-order density in g space
1519 : !> \param tau1_r ...
1520 : !> \param pw_pool pw pool to create new grids
1521 : !> \param xc_section XC section to calculate the derivatives from
1522 : !> \param gapw whether to carry out GAPW (not possible with numerical derivatives)
1523 : !> \param vxg GAPW potential
1524 : !> \param do_excitations ...
1525 : !> \param do_triplet ...
1526 : !> \param compute_virial ...
1527 : !> \param virial_xc virial terms will be collected here
1528 : ! **************************************************************************************************
1529 14420 : SUBROUTINE xc_calc_2nd_deriv(v_xc, v_xc_tau, deriv_set, rho_set, rho1_r, rho1_g, tau1_r, &
1530 : pw_pool, xc_section, gapw, vxg, &
1531 : do_excitations, do_sf, do_triplet, &
1532 : compute_virial, virial_xc)
1533 :
1534 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_xc, v_xc_tau
1535 : TYPE(xc_derivative_set_type) :: deriv_set
1536 : TYPE(xc_rho_set_type) :: rho_set
1537 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho1_r, tau1_r
1538 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho1_g
1539 : TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool
1540 : TYPE(section_vals_type), INTENT(IN), POINTER :: xc_section
1541 : LOGICAL, INTENT(IN) :: gapw
1542 : REAL(KIND=dp), DIMENSION(:, :, :, :), OPTIONAL, &
1543 : POINTER :: vxg
1544 : LOGICAL, INTENT(IN), OPTIONAL :: do_excitations, do_sf, &
1545 : do_triplet, compute_virial
1546 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT), &
1547 : OPTIONAL :: virial_xc
1548 :
1549 : CHARACTER(len=*), PARAMETER :: routineN = 'xc_calc_2nd_deriv'
1550 :
1551 : INTEGER :: handle, ispin, nspins
1552 : INTEGER, DIMENSION(2, 3) :: bo
1553 : LOGICAL :: lsd, my_compute_virial, &
1554 : my_do_excitations, my_do_sf, &
1555 : my_do_triplet
1556 : REAL(KIND=dp) :: fac
1557 : TYPE(section_vals_type), POINTER :: xc_fun_section
1558 : TYPE(xc_rho_cflags_type) :: needs
1559 : TYPE(xc_rho_set_type) :: rho1_set
1560 :
1561 14420 : CALL timeset(routineN, handle)
1562 :
1563 14420 : my_compute_virial = .FALSE.
1564 14420 : IF (PRESENT(compute_virial)) my_compute_virial = compute_virial
1565 :
1566 14420 : my_do_sf = .FALSE.
1567 14420 : IF (PRESENT(do_sf)) my_do_sf = do_sf
1568 :
1569 14420 : my_do_excitations = .FALSE.
1570 14420 : IF (PRESENT(do_excitations)) my_do_excitations = do_excitations
1571 :
1572 14420 : my_do_triplet = .FALSE.
1573 14420 : IF (PRESENT(do_triplet)) my_do_triplet = do_triplet
1574 :
1575 14420 : nspins = SIZE(rho1_r)
1576 14420 : lsd = (nspins == 2)
1577 14420 : IF (nspins == 1 .AND. my_do_excitations .AND. my_do_triplet) THEN
1578 0 : nspins = 2
1579 0 : lsd = .TRUE.
1580 14420 : ELSE IF (my_do_sf) THEN
1581 104 : nspins = 1
1582 104 : lsd = .TRUE.
1583 : END IF
1584 :
1585 14420 : NULLIFY (v_xc, v_xc_tau)
1586 59674 : ALLOCATE (v_xc(nspins))
1587 30834 : DO ispin = 1, nspins
1588 16414 : CALL pw_pool%create_pw(v_xc(ispin))
1589 30834 : CALL pw_zero(v_xc(ispin))
1590 : END DO
1591 :
1592 14420 : xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
1593 14420 : needs = xc_functionals_get_needs(xc_fun_section, lsd, .TRUE.)
1594 :
1595 14420 : IF (needs%tau .OR. needs%tau_spin) THEN
1596 536 : IF (.NOT. ASSOCIATED(tau1_r)) &
1597 0 : CPABORT("Tau-dependent functionals requires allocated kinetic energy density grid")
1598 1736 : ALLOCATE (v_xc_tau(nspins))
1599 1200 : DO ispin = 1, nspins
1600 664 : CALL pw_pool%create_pw(v_xc_tau(ispin))
1601 15084 : CALL pw_zero(v_xc_tau(ispin))
1602 : END DO
1603 : END IF
1604 :
1605 14420 : IF (section_get_lval(xc_section, "2ND_DERIV_ANALYTICAL")) THEN
1606 : !------!
1607 : ! rho1 !
1608 : !------!
1609 141200 : bo = rho1_r(1)%pw_grid%bounds_local
1610 : ! create the place where to store the argument for the functionals
1611 : CALL xc_rho_set_create(rho1_set, bo, &
1612 : rho_cutoff=section_get_rval(xc_section, "DENSITY_CUTOFF"), &
1613 : drho_cutoff=section_get_rval(xc_section, "GRADIENT_CUTOFF"), &
1614 14120 : tau_cutoff=section_get_rval(xc_section, "TAU_CUTOFF"))
1615 :
1616 : ! calculate the arguments needed by the functionals
1617 : CALL xc_rho_set_update(rho1_set, rho1_r, rho1_g, tau1_r, needs, &
1618 : section_get_ival(xc_section, "XC_GRID%XC_DERIV"), &
1619 : section_get_ival(xc_section, "XC_GRID%XC_SMOOTH_RHO"), &
1620 14120 : pw_pool, spinflip=my_do_sf)
1621 :
1622 14120 : fac = 0._dp
1623 14120 : IF (nspins == 1 .AND. my_do_excitations) THEN
1624 1124 : IF (my_do_triplet) fac = -1.0_dp
1625 : END IF
1626 :
1627 : CALL xc_calc_2nd_deriv_analytical(v_xc, v_xc_tau, deriv_set, rho_set, &
1628 : rho1_set, pw_pool, xc_section, gapw, vxg=vxg, spinflip=my_do_sf, &
1629 14120 : tddfpt_fac=fac, compute_virial=compute_virial, virial_xc=virial_xc)
1630 :
1631 14120 : CALL xc_rho_set_release(rho1_set)
1632 :
1633 : ELSE
1634 300 : IF (gapw) CPABORT("Numerical 2nd derivatives not implemented with GAPW")
1635 :
1636 : CALL xc_calc_2nd_deriv_numerical(v_xc, v_xc_tau, rho_set, rho1_r, rho1_g, tau1_r, &
1637 : pw_pool, xc_section, &
1638 : my_do_excitations .AND. my_do_triplet, &
1639 300 : compute_virial, virial_xc, deriv_set)
1640 : END IF
1641 :
1642 14420 : CALL timestop(handle)
1643 :
1644 317240 : END SUBROUTINE xc_calc_2nd_deriv
1645 :
1646 : ! **************************************************************************************************
1647 : !> \brief calculates 2nd derivative numerically
1648 : !> \param v_xc potential to be calculated (has to be allocated already)
1649 : !> \param v_tau tau-part of the potential to be calculated (has to be allocated already)
1650 : !> \param rho_set KS density from xc_prep_2nd_deriv
1651 : !> \param rho1_r first-order density in r-space
1652 : !> \param rho1_g first-order density in g-space
1653 : !> \param tau1_r first-order kinetic-energy density in r-space
1654 : !> \param pw_pool pw pool for new grids
1655 : !> \param xc_section XC section to calculate the derivatives from
1656 : !> \param do_triplet ...
1657 : !> \param calc_virial whether to calculate virial terms
1658 : !> \param virial_xc collects stress tensor components (no metaGGAs!)
1659 : !> \param deriv_set deriv set from xc_prep_2nd_deriv (only for virials)
1660 : ! **************************************************************************************************
1661 340 : SUBROUTINE xc_calc_2nd_deriv_numerical(v_xc, v_tau, rho_set, rho1_r, rho1_g, tau1_r, &
1662 : pw_pool, xc_section, &
1663 : do_triplet, calc_virial, virial_xc, deriv_set)
1664 :
1665 : TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(IN), POINTER :: v_xc, v_tau
1666 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
1667 : TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(IN), POINTER :: rho1_r, tau1_r
1668 : TYPE(pw_c1d_gs_type), DIMENSION(:), INTENT(IN), POINTER :: rho1_g
1669 : TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool
1670 : TYPE(section_vals_type), INTENT(IN), POINTER :: xc_section
1671 : LOGICAL, INTENT(IN) :: do_triplet
1672 : LOGICAL, INTENT(IN), OPTIONAL :: calc_virial
1673 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT), &
1674 : OPTIONAL :: virial_xc
1675 : TYPE(xc_derivative_set_type), OPTIONAL :: deriv_set
1676 :
1677 : CHARACTER(len=*), PARAMETER :: routineN = 'xc_calc_2nd_deriv_numerical'
1678 : REAL(KIND=dp), DIMENSION(-4:4, 4), PARAMETER :: &
1679 : 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, &
1680 : 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, &
1681 : 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, &
1682 : 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])
1683 :
1684 : INTEGER :: handle, idir, ispin, nspins, istep, nsteps
1685 : INTEGER, DIMENSION(2, 3) :: bo
1686 : LOGICAL :: gradient_f, lsd, my_calc_virial, tau_f, laplace_f, rho_f
1687 : REAL(KIND=dp) :: exc, gradient_cut, h, weight, step, rho_cutoff
1688 340 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: dr1dr, dra1dra, drb1drb
1689 : REAL(KIND=dp), DIMENSION(3, 3) :: virial_dummy
1690 340 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: norm_drho, norm_drho2, norm_drho2a, &
1691 340 : norm_drho2b, norm_drhoa, norm_drhob, &
1692 1020 : rho, rho1, rho1a, rho1b, rhoa, rhob, &
1693 680 : tau_a, tau_b, tau, tau1, tau1a, tau1b, laplace, laplace1, &
1694 340 : laplacea, laplaceb, laplace1a, laplace1b, &
1695 680 : laplace2, laplace2a, laplace2b, deriv_data
1696 8160 : TYPE(cp_3d_r_cp_type), DIMENSION(3) :: drho, drho1, drho1a, drho1b, drhoa, drhob
1697 : TYPE(pw_r3d_rs_type) :: v_drho, v_drhoa, v_drhob
1698 340 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: vxc_rho, vxc_tau
1699 340 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
1700 340 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, tau_r
1701 : TYPE(pw_r3d_rs_type) :: virial_pw, v_laplace, v_laplacea, v_laplaceb
1702 : TYPE(section_vals_type), POINTER :: xc_fun_section
1703 : TYPE(xc_derivative_set_type) :: deriv_set1
1704 : TYPE(xc_rho_cflags_type) :: needs
1705 : TYPE(xc_rho_set_type) :: rho1_set, rho2_set
1706 :
1707 340 : CALL timeset(routineN, handle)
1708 :
1709 340 : my_calc_virial = .FALSE.
1710 340 : IF (PRESENT(calc_virial) .AND. PRESENT(virial_xc)) my_calc_virial = calc_virial
1711 :
1712 340 : nspins = SIZE(v_xc)
1713 :
1714 340 : NULLIFY (tau, tau_r, tau_a, tau_b)
1715 :
1716 340 : h = section_get_rval(xc_section, "STEP_SIZE")
1717 340 : nsteps = section_get_ival(xc_section, "NSTEPS")
1718 340 : IF (nsteps < LBOUND(weights, 2) .OR. nspins > UBOUND(weights, 2)) THEN
1719 0 : CPABORT("The number of steps must be a value from 1 to 4.")
1720 : END IF
1721 :
1722 340 : IF (nspins == 2) THEN
1723 148 : NULLIFY (vxc_rho, rho_g, vxc_tau)
1724 444 : ALLOCATE (rho_r(2))
1725 444 : DO ispin = 1, nspins
1726 444 : CALL pw_pool%create_pw(rho_r(ispin))
1727 : END DO
1728 148 : IF (ASSOCIATED(tau1_r) .AND. ASSOCIATED(v_tau)) THEN
1729 162 : ALLOCATE (tau_r(2))
1730 162 : DO ispin = 1, nspins
1731 162 : CALL pw_pool%create_pw(tau_r(ispin))
1732 : END DO
1733 : END IF
1734 148 : CALL xc_rho_set_get(rho_set, can_return_null=.TRUE., rhoa=rhoa, rhob=rhob, tau_a=tau_a, tau_b=tau_b)
1735 1088 : DO istep = -nsteps, nsteps
1736 940 : IF (istep == 0) CYCLE
1737 792 : weight = weights(istep, nsteps)/h
1738 792 : step = REAL(istep, dp)*h
1739 : CALL calc_resp_potential_numer_ab(rho_r, rho_g, rho1_r, rhoa, rhob, vxc_rho, &
1740 792 : tau_r, tau1_r, tau_a, tau_b, vxc_tau, xc_section, pw_pool, step)
1741 2376 : DO ispin = 1, nspins
1742 1584 : CALL pw_axpy(vxc_rho(ispin), v_xc(ispin), weight)
1743 2376 : IF (ASSOCIATED(vxc_tau) .AND. ASSOCIATED(v_tau)) THEN
1744 456 : CALL pw_axpy(vxc_tau(ispin), v_tau(ispin), weight)
1745 : END IF
1746 : END DO
1747 2376 : DO ispin = 1, nspins
1748 2376 : CALL vxc_rho(ispin)%release()
1749 : END DO
1750 792 : DEALLOCATE (vxc_rho)
1751 940 : IF (ASSOCIATED(vxc_tau)) THEN
1752 684 : DO ispin = 1, nspins
1753 684 : CALL vxc_tau(ispin)%release()
1754 : END DO
1755 228 : DEALLOCATE (vxc_tau)
1756 : END IF
1757 : END DO
1758 192 : ELSE IF (nspins == 1 .AND. do_triplet) THEN
1759 20 : NULLIFY (vxc_rho, vxc_tau, rho_g)
1760 60 : ALLOCATE (rho_r(2))
1761 60 : DO ispin = 1, 2
1762 60 : CALL pw_pool%create_pw(rho_r(ispin))
1763 : END DO
1764 20 : IF (ASSOCIATED(tau1_r) .AND. ASSOCIATED(v_tau)) THEN
1765 0 : ALLOCATE (tau_r(2))
1766 0 : DO ispin = 1, nspins
1767 0 : CALL pw_pool%create_pw(tau_r(ispin))
1768 : END DO
1769 : END IF
1770 20 : CALL xc_rho_set_get(rho_set, can_return_null=.TRUE., rhoa=rhoa, rhob=rhob, tau_a=tau_a, tau_b=tau_b)
1771 160 : DO istep = -nsteps, nsteps
1772 140 : IF (istep == 0) CYCLE
1773 120 : weight = weights(istep, nsteps)/h
1774 120 : step = REAL(istep, dp)*h
1775 : ! K(alpha,alpha)
1776 120 : !$OMP PARALLEL DEFAULT(NONE) SHARED(rho_r,rhoa,rhob,step,rho1_r,tau_r,tau_a,tau_b,tau1_r)
1777 : !$OMP WORKSHARE
1778 : rho_r(1)%array(:, :, :) = rhoa(:, :, :) + step*rho1_r(1)%array(:, :, :)
1779 : !$OMP END WORKSHARE NOWAIT
1780 : !$OMP WORKSHARE
1781 : rho_r(2)%array(:, :, :) = rhob(:, :, :)
1782 : !$OMP END WORKSHARE NOWAIT
1783 : IF (ASSOCIATED(tau1_r)) THEN
1784 : !$OMP WORKSHARE
1785 : tau_r(1)%array(:, :, :) = tau_a(:, :, :) + step*tau1_r(1)%array(:, :, :)
1786 : !$OMP END WORKSHARE NOWAIT
1787 : !$OMP WORKSHARE
1788 : tau_r(2)%array(:, :, :) = tau_b(:, :, :)
1789 : !$OMP END WORKSHARE NOWAIT
1790 : END IF
1791 : !$OMP END PARALLEL
1792 : CALL xc_vxc_pw_create(vxc_rho, vxc_tau, exc, rho_r, rho_g, tau_r, xc_section, &
1793 120 : pw_pool, .FALSE., virial_dummy)
1794 120 : CALL pw_axpy(vxc_rho(1), v_xc(1), weight)
1795 120 : IF (ASSOCIATED(vxc_tau) .AND. ASSOCIATED(v_tau)) THEN
1796 0 : CALL pw_axpy(vxc_tau(1), v_tau(1), weight)
1797 : END IF
1798 360 : DO ispin = 1, 2
1799 360 : CALL vxc_rho(ispin)%release()
1800 : END DO
1801 120 : DEALLOCATE (vxc_rho)
1802 120 : IF (ASSOCIATED(vxc_tau)) THEN
1803 0 : DO ispin = 1, 2
1804 0 : CALL vxc_tau(ispin)%release()
1805 : END DO
1806 0 : DEALLOCATE (vxc_tau)
1807 : END IF
1808 120 : !$OMP PARALLEL DEFAULT(NONE) SHARED(rho_r,rhoa,rhob,step,rho1_r,tau_r,tau_a,tau_b,tau1_r)
1809 : !$OMP WORKSHARE
1810 : ! K(alpha,beta)
1811 : rho_r(1)%array(:, :, :) = rhoa(:, :, :)
1812 : !$OMP END WORKSHARE NOWAIT
1813 : !$OMP WORKSHARE
1814 : rho_r(2)%array(:, :, :) = rhob(:, :, :) + step*rho1_r(1)%array(:, :, :)
1815 : !$OMP END WORKSHARE NOWAIT
1816 : IF (ASSOCIATED(tau1_r)) THEN
1817 : !$OMP WORKSHARE
1818 : tau_r(1)%array(:, :, :) = tau_a(:, :, :)
1819 : !$OMP END WORKSHARE NOWAIT
1820 : !$OMP WORKSHARE
1821 : tau_r(2)%array(:, :, :) = tau_b(:, :, :) + step*tau1_r(1)%array(:, :, :)
1822 : !$OMP END WORKSHARE NOWAIT
1823 : END IF
1824 : !$OMP END PARALLEL
1825 : CALL xc_vxc_pw_create(vxc_rho, vxc_tau, exc, rho_r, rho_g, tau_r, xc_section, &
1826 120 : pw_pool, .FALSE., virial_dummy)
1827 120 : CALL pw_axpy(vxc_rho(1), v_xc(1), weight)
1828 120 : IF (ASSOCIATED(vxc_tau) .AND. ASSOCIATED(v_tau)) THEN
1829 0 : CALL pw_axpy(vxc_tau(1), v_tau(1), weight)
1830 : END IF
1831 360 : DO ispin = 1, 2
1832 360 : CALL vxc_rho(ispin)%release()
1833 : END DO
1834 120 : DEALLOCATE (vxc_rho)
1835 140 : IF (ASSOCIATED(vxc_tau)) THEN
1836 0 : DO ispin = 1, 2
1837 0 : CALL vxc_tau(ispin)%release()
1838 : END DO
1839 0 : DEALLOCATE (vxc_tau)
1840 : END IF
1841 : END DO
1842 : ELSE
1843 172 : NULLIFY (vxc_rho, rho_r, rho_g, vxc_tau, tau_r, tau)
1844 344 : ALLOCATE (rho_r(1))
1845 172 : CALL pw_pool%create_pw(rho_r(1))
1846 172 : IF (ASSOCIATED(tau1_r) .AND. ASSOCIATED(v_tau)) THEN
1847 92 : ALLOCATE (tau_r(1))
1848 46 : CALL pw_pool%create_pw(tau_r(1))
1849 : END IF
1850 172 : CALL xc_rho_set_get(rho_set, can_return_null=.TRUE., rho=rho, tau=tau)
1851 1280 : DO istep = -nsteps, nsteps
1852 1108 : IF (istep == 0) CYCLE
1853 936 : weight = weights(istep, nsteps)/h
1854 936 : step = REAL(istep, dp)*h
1855 936 : !$OMP PARALLEL DEFAULT(NONE) SHARED(rho_r,rho,step,rho1_r,tau1_r,tau,tau_r)
1856 : !$OMP WORKSHARE
1857 : rho_r(1)%array(:, :, :) = rho(:, :, :) + step*rho1_r(1)%array(:, :, :)
1858 : !$OMP END WORKSHARE NOWAIT
1859 : IF (ASSOCIATED(tau1_r) .AND. ASSOCIATED(tau) .AND. ASSOCIATED(tau1_r)) THEN
1860 : !$OMP WORKSHARE
1861 : tau_r(1)%array(:, :, :) = tau(:, :, :) + step*tau1_r(1)%array(:, :, :)
1862 : !$OMP END WORKSHARE NOWAIT
1863 : END IF
1864 : !$OMP END PARALLEL
1865 : CALL xc_vxc_pw_create(vxc_rho, vxc_tau, exc, rho_r, rho_g, tau_r, xc_section, &
1866 936 : pw_pool, .FALSE., virial_dummy)
1867 936 : CALL pw_axpy(vxc_rho(1), v_xc(1), weight)
1868 936 : IF (ASSOCIATED(vxc_tau) .AND. ASSOCIATED(v_tau)) THEN
1869 276 : CALL pw_axpy(vxc_tau(1), v_tau(1), weight)
1870 : END IF
1871 936 : CALL vxc_rho(1)%release()
1872 936 : DEALLOCATE (vxc_rho)
1873 1108 : IF (ASSOCIATED(vxc_tau)) THEN
1874 276 : CALL vxc_tau(1)%release()
1875 276 : DEALLOCATE (vxc_tau)
1876 : END IF
1877 : END DO
1878 : END IF
1879 :
1880 340 : IF (my_calc_virial) THEN
1881 36 : lsd = (nspins == 2)
1882 36 : IF (nspins == 1 .AND. do_triplet) THEN
1883 0 : lsd = .TRUE.
1884 : END IF
1885 :
1886 36 : CALL check_for_derivatives(deriv_set, (nspins == 2), rho_f, gradient_f, tau_f, laplace_f)
1887 :
1888 : ! Calculate the virial terms
1889 : ! Those arising from the first derivatives are treated like in xc_calc_2nd_deriv_analytical
1890 : ! Those arising from the second derivatives are calculated numerically
1891 : ! We assume that all metaGGA functionals require the gradient
1892 36 : IF (gradient_f) THEN
1893 360 : bo = rho_set%local_bounds
1894 :
1895 : ! Create the work grid for the virial terms
1896 36 : CALL allocate_pw(virial_pw, pw_pool, bo)
1897 :
1898 36 : gradient_cut = section_get_rval(xc_section, "GRADIENT_CUTOFF")
1899 :
1900 : ! create the container to store the argument of the functionals
1901 : CALL xc_rho_set_create(rho1_set, bo, &
1902 : rho_cutoff=section_get_rval(xc_section, "DENSITY_CUTOFF"), &
1903 : drho_cutoff=gradient_cut, &
1904 36 : tau_cutoff=section_get_rval(xc_section, "TAU_CUTOFF"))
1905 :
1906 36 : xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
1907 36 : needs = xc_functionals_get_needs(xc_fun_section, lsd, .TRUE.)
1908 :
1909 : ! calculate the arguments needed by the functionals
1910 : CALL xc_rho_set_update(rho1_set, rho1_r, rho1_g, tau1_r, needs, &
1911 : section_get_ival(xc_section, "XC_GRID%XC_DERIV"), &
1912 : section_get_ival(xc_section, "XC_GRID%XC_SMOOTH_RHO"), &
1913 36 : pw_pool)
1914 :
1915 36 : IF (lsd) THEN
1916 : CALL xc_rho_set_get(rho_set, drhoa=drhoa, drhob=drhob, norm_drho=norm_drho, &
1917 : norm_drhoa=norm_drhoa, norm_drhob=norm_drhob, tau_a=tau_a, tau_b=tau_b, &
1918 10 : laplace_rhoa=laplacea, laplace_rhob=laplaceb, can_return_null=.TRUE.)
1919 : CALL xc_rho_set_get(rho1_set, rhoa=rho1a, rhob=rho1b, drhoa=drho1a, drhob=drho1b, laplace_rhoa=laplace1a, &
1920 10 : laplace_rhob=laplace1b, can_return_null=.TRUE.)
1921 :
1922 10 : CALL calc_drho_from_ab(drho, drhoa, drhob)
1923 10 : CALL calc_drho_from_ab(drho1, drho1a, drho1b)
1924 : ELSE
1925 26 : CALL xc_rho_set_get(rho_set, drho=drho, norm_drho=norm_drho, tau=tau, laplace_rho=laplace, can_return_null=.TRUE.)
1926 26 : CALL xc_rho_set_get(rho1_set, rho=rho1, drho=drho1, laplace_rho=laplace1, can_return_null=.TRUE.)
1927 : END IF
1928 :
1929 36 : CALL prepare_dr1dr(dr1dr, drho, drho1)
1930 :
1931 36 : IF (lsd) THEN
1932 10 : CALL prepare_dr1dr(dra1dra, drhoa, drho1a)
1933 10 : CALL prepare_dr1dr(drb1drb, drhob, drho1b)
1934 :
1935 10 : CALL allocate_pw(v_drho, pw_pool, bo)
1936 10 : CALL allocate_pw(v_drhoa, pw_pool, bo)
1937 10 : CALL allocate_pw(v_drhob, pw_pool, bo)
1938 :
1939 10 : IF (ASSOCIATED(norm_drhoa)) CALL apply_drho(deriv_set, [deriv_norm_drhoa], virial_pw, drhoa, drho1a, virial_xc, &
1940 10 : norm_drhoa, gradient_cut, dra1dra, v_drhoa%array)
1941 10 : IF (ASSOCIATED(norm_drhob)) CALL apply_drho(deriv_set, [deriv_norm_drhob], virial_pw, drhob, drho1b, virial_xc, &
1942 10 : norm_drhob, gradient_cut, drb1drb, v_drhob%array)
1943 10 : IF (ASSOCIATED(norm_drho)) CALL apply_drho(deriv_set, [deriv_norm_drho], virial_pw, drho, drho1, virial_xc, &
1944 6 : norm_drho, gradient_cut, dr1dr, v_drho%array)
1945 10 : IF (laplace_f) THEN
1946 2 : CALL xc_derivative_get(xc_dset_get_derivative(deriv_set, [deriv_laplace_rhoa]), deriv_data=deriv_data)
1947 2 : CPASSERT(ASSOCIATED(deriv_data))
1948 15026 : virial_pw%array(:, :, :) = -rho1a(:, :, :)
1949 2 : CALL virial_laplace(virial_pw, pw_pool, virial_xc, deriv_data)
1950 :
1951 2 : CALL allocate_pw(v_laplacea, pw_pool, bo)
1952 :
1953 2 : CALL xc_derivative_get(xc_dset_get_derivative(deriv_set, [deriv_laplace_rhob]), deriv_data=deriv_data)
1954 2 : CPASSERT(ASSOCIATED(deriv_data))
1955 15026 : virial_pw%array(:, :, :) = -rho1b(:, :, :)
1956 2 : CALL virial_laplace(virial_pw, pw_pool, virial_xc, deriv_data)
1957 :
1958 2 : CALL allocate_pw(v_laplaceb, pw_pool, bo)
1959 : END IF
1960 :
1961 : ELSE
1962 :
1963 : ! Create the work grid for the potential of the gradient part
1964 26 : CALL allocate_pw(v_drho, pw_pool, bo)
1965 :
1966 : CALL apply_drho(deriv_set, [deriv_norm_drho], virial_pw, drho, drho1, virial_xc, &
1967 26 : norm_drho, gradient_cut, dr1dr, v_drho%array)
1968 26 : IF (laplace_f) THEN
1969 2 : CALL xc_derivative_get(xc_dset_get_derivative(deriv_set, [deriv_laplace_rho]), deriv_data=deriv_data)
1970 2 : CPASSERT(ASSOCIATED(deriv_data))
1971 28862 : virial_pw%array(:, :, :) = -rho1(:, :, :)
1972 2 : CALL virial_laplace(virial_pw, pw_pool, virial_xc, deriv_data)
1973 :
1974 2 : CALL allocate_pw(v_laplace, pw_pool, bo)
1975 : END IF
1976 :
1977 : END IF
1978 :
1979 36 : IF (lsd) THEN
1980 150260 : rho_r(1)%array = rhoa
1981 150260 : rho_r(2)%array = rhob
1982 : ELSE
1983 701552 : rho_r(1)%array = rho
1984 : END IF
1985 36 : IF (ASSOCIATED(tau1_r)) THEN
1986 8 : IF (lsd) THEN
1987 60104 : tau_r(1)%array = tau_a
1988 60104 : tau_r(2)%array = tau_b
1989 : ELSE
1990 115448 : tau_r(1)%array = tau
1991 : END IF
1992 : END IF
1993 :
1994 : ! Create deriv sets with same densities but different gradients
1995 36 : CALL xc_dset_create(deriv_set1, pw_pool)
1996 :
1997 36 : rho_cutoff = section_get_rval(xc_section, "DENSITY_CUTOFF")
1998 :
1999 : ! create the place where to store the argument for the functionals
2000 : CALL xc_rho_set_create(rho2_set, bo, &
2001 : rho_cutoff=rho_cutoff, &
2002 : drho_cutoff=section_get_rval(xc_section, "GRADIENT_CUTOFF"), &
2003 36 : tau_cutoff=section_get_rval(xc_section, "TAU_CUTOFF"))
2004 :
2005 : ! calculate the arguments needed by the functionals
2006 : CALL xc_rho_set_update(rho2_set, rho_r, rho_g, tau_r, needs, &
2007 : section_get_ival(xc_section, "XC_GRID%XC_DERIV"), &
2008 : section_get_ival(xc_section, "XC_GRID%XC_SMOOTH_RHO"), &
2009 36 : pw_pool)
2010 :
2011 36 : IF (lsd) THEN
2012 : CALL xc_rho_set_get(rho1_set, rhoa=rho1a, rhob=rho1b, tau_a=tau1a, tau_b=tau1b, &
2013 10 : laplace_rhoa=laplace1a, laplace_rhob=laplace1b, can_return_null=.TRUE.)
2014 : CALL xc_rho_set_get(rho2_set, norm_drhoa=norm_drho2a, norm_drhob=norm_drho2b, &
2015 10 : norm_drho=norm_drho2, laplace_rhoa=laplace2a, laplace_rhob=laplace2b, can_return_null=.TRUE.)
2016 :
2017 64 : DO istep = -nsteps, nsteps
2018 54 : IF (istep == 0) CYCLE
2019 44 : weight = weights(istep, nsteps)/h
2020 44 : step = REAL(istep, dp)*h
2021 44 : IF (ASSOCIATED(norm_drhoa)) THEN
2022 44 : CALL get_derivs_rho(norm_drho2a, norm_drhoa, step, xc_fun_section, lsd, rho2_set, deriv_set1)
2023 : CALL update_deriv_rho(deriv_set1, [deriv_rhoa], bo, &
2024 44 : norm_drhoa, gradient_cut, weight, rho1a, v_drhoa%array)
2025 : CALL update_deriv_rho(deriv_set1, [deriv_rhob], bo, &
2026 44 : norm_drhoa, gradient_cut, weight, rho1b, v_drhoa%array)
2027 : CALL update_deriv_rho(deriv_set1, [deriv_norm_drhoa], bo, &
2028 44 : norm_drhoa, gradient_cut, weight, dra1dra, v_drhoa%array)
2029 : CALL update_deriv_drho_ab(deriv_set1, [deriv_norm_drhob], bo, &
2030 44 : norm_drhoa, gradient_cut, weight, dra1dra, drb1drb, v_drhoa%array, v_drhob%array)
2031 : CALL update_deriv_drho_ab(deriv_set1, [deriv_norm_drho], bo, &
2032 44 : norm_drhoa, gradient_cut, weight, dra1dra, dr1dr, v_drhoa%array, v_drho%array)
2033 44 : IF (tau_f) THEN
2034 : CALL update_deriv_rho(deriv_set1, [deriv_tau_a], bo, &
2035 8 : norm_drhoa, gradient_cut, weight, tau1a, v_drhoa%array)
2036 : CALL update_deriv_rho(deriv_set1, [deriv_tau_b], bo, &
2037 8 : norm_drhoa, gradient_cut, weight, tau1b, v_drhoa%array)
2038 : END IF
2039 44 : IF (laplace_f) THEN
2040 : CALL update_deriv_rho(deriv_set1, [deriv_laplace_rhoa], bo, &
2041 4 : norm_drhoa, gradient_cut, weight, laplace1a, v_drhoa%array)
2042 : CALL update_deriv_rho(deriv_set1, [deriv_laplace_rhob], bo, &
2043 4 : norm_drhoa, gradient_cut, weight, laplace1b, v_drhoa%array)
2044 : END IF
2045 : END IF
2046 :
2047 44 : IF (ASSOCIATED(norm_drhob)) THEN
2048 44 : CALL get_derivs_rho(norm_drho2b, norm_drhob, step, xc_fun_section, lsd, rho2_set, deriv_set1)
2049 : CALL update_deriv_rho(deriv_set1, [deriv_rhoa], bo, &
2050 44 : norm_drhob, gradient_cut, weight, rho1a, v_drhob%array)
2051 : CALL update_deriv_rho(deriv_set1, [deriv_rhob], bo, &
2052 44 : norm_drhob, gradient_cut, weight, rho1b, v_drhob%array)
2053 : CALL update_deriv_rho(deriv_set1, [deriv_norm_drhob], bo, &
2054 44 : norm_drhob, gradient_cut, weight, drb1drb, v_drhob%array)
2055 : CALL update_deriv_drho_ab(deriv_set1, [deriv_norm_drhoa], bo, &
2056 44 : norm_drhob, gradient_cut, weight, drb1drb, dra1dra, v_drhob%array, v_drhoa%array)
2057 : CALL update_deriv_drho_ab(deriv_set1, [deriv_norm_drho], bo, &
2058 44 : norm_drhob, gradient_cut, weight, drb1drb, dr1dr, v_drhob%array, v_drho%array)
2059 44 : IF (tau_f) THEN
2060 : CALL update_deriv_rho(deriv_set1, [deriv_tau_a], bo, &
2061 8 : norm_drhob, gradient_cut, weight, tau1a, v_drhob%array)
2062 : CALL update_deriv_rho(deriv_set1, [deriv_tau_b], bo, &
2063 8 : norm_drhob, gradient_cut, weight, tau1b, v_drhob%array)
2064 : END IF
2065 44 : IF (laplace_f) THEN
2066 : CALL update_deriv_rho(deriv_set1, [deriv_laplace_rhoa], bo, &
2067 4 : norm_drhob, gradient_cut, weight, laplace1a, v_drhob%array)
2068 : CALL update_deriv_rho(deriv_set1, [deriv_laplace_rhob], bo, &
2069 4 : norm_drhob, gradient_cut, weight, laplace1b, v_drhob%array)
2070 : END IF
2071 : END IF
2072 :
2073 44 : IF (ASSOCIATED(norm_drho)) THEN
2074 20 : CALL get_derivs_rho(norm_drho2, norm_drho, step, xc_fun_section, lsd, rho2_set, deriv_set1)
2075 : CALL update_deriv_rho(deriv_set1, [deriv_rhoa], bo, &
2076 20 : norm_drho, gradient_cut, weight, rho1a, v_drho%array)
2077 : CALL update_deriv_rho(deriv_set1, [deriv_rhob], bo, &
2078 20 : norm_drho, gradient_cut, weight, rho1b, v_drho%array)
2079 : CALL update_deriv_rho(deriv_set1, [deriv_norm_drho], bo, &
2080 20 : norm_drho, gradient_cut, weight, dr1dr, v_drho%array)
2081 : CALL update_deriv_drho_ab(deriv_set1, [deriv_norm_drhoa], bo, &
2082 20 : norm_drho, gradient_cut, weight, dr1dr, dra1dra, v_drho%array, v_drhoa%array)
2083 : CALL update_deriv_drho_ab(deriv_set1, [deriv_norm_drhob], bo, &
2084 20 : norm_drho, gradient_cut, weight, dr1dr, drb1drb, v_drho%array, v_drhob%array)
2085 20 : IF (tau_f) THEN
2086 : CALL update_deriv_rho(deriv_set1, [deriv_tau_a], bo, &
2087 8 : norm_drho, gradient_cut, weight, tau1a, v_drho%array)
2088 : CALL update_deriv_rho(deriv_set1, [deriv_tau_b], bo, &
2089 8 : norm_drho, gradient_cut, weight, tau1b, v_drho%array)
2090 : END IF
2091 20 : IF (laplace_f) THEN
2092 : CALL update_deriv_rho(deriv_set1, [deriv_laplace_rhoa], bo, &
2093 4 : norm_drho, gradient_cut, weight, laplace1a, v_drho%array)
2094 : CALL update_deriv_rho(deriv_set1, [deriv_laplace_rhob], bo, &
2095 4 : norm_drho, gradient_cut, weight, laplace1b, v_drho%array)
2096 : END IF
2097 : END IF
2098 :
2099 54 : IF (laplace_f) THEN
2100 :
2101 4 : CALL get_derivs_rho(laplace2a, laplacea, step, xc_fun_section, lsd, rho2_set, deriv_set1)
2102 :
2103 : ! Obtain the numerical 2nd derivatives w.r.t. to drho and collect the potential
2104 : CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [deriv_rhoa], bo, &
2105 4 : weight, rho1a, v_laplacea%array)
2106 : CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [deriv_rhob], bo, &
2107 4 : weight, rho1b, v_laplacea%array)
2108 4 : IF (ASSOCIATED(norm_drho)) THEN
2109 : CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [deriv_norm_drho], bo, &
2110 4 : weight, dr1dr, v_laplacea%array)
2111 : END IF
2112 4 : IF (ASSOCIATED(norm_drhoa)) THEN
2113 : CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [deriv_norm_drhoa], bo, &
2114 4 : weight, dra1dra, v_laplacea%array)
2115 : END IF
2116 4 : IF (ASSOCIATED(norm_drhob)) THEN
2117 : CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [deriv_norm_drhob], bo, &
2118 4 : weight, drb1drb, v_laplacea%array)
2119 : END IF
2120 :
2121 4 : IF (ASSOCIATED(tau1a)) THEN
2122 : CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [deriv_tau_a], bo, &
2123 4 : weight, tau1a, v_laplacea%array)
2124 : END IF
2125 4 : IF (ASSOCIATED(tau1b)) THEN
2126 : CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [deriv_tau_b], bo, &
2127 4 : weight, tau1b, v_laplacea%array)
2128 : END IF
2129 :
2130 : CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [deriv_laplace_rhoa], bo, &
2131 4 : weight, laplace1a, v_laplacea%array)
2132 :
2133 : CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [deriv_laplace_rhob], bo, &
2134 4 : weight, laplace1b, v_laplacea%array)
2135 :
2136 : ! The same for the beta spin
2137 4 : CALL get_derivs_rho(laplace2b, laplaceb, step, xc_fun_section, lsd, rho2_set, deriv_set1)
2138 :
2139 : ! Obtain the numerical 2nd derivatives w.r.t. to drho and collect the potential
2140 : CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [deriv_rhoa], bo, &
2141 4 : weight, rho1a, v_laplaceb%array)
2142 : CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [deriv_rhob], bo, &
2143 4 : weight, rho1b, v_laplaceb%array)
2144 4 : IF (ASSOCIATED(norm_drho)) THEN
2145 : CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [deriv_norm_drho], bo, &
2146 4 : weight, dr1dr, v_laplaceb%array)
2147 : END IF
2148 4 : IF (ASSOCIATED(norm_drhoa)) THEN
2149 : CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [deriv_norm_drhoa], bo, &
2150 4 : weight, dra1dra, v_laplaceb%array)
2151 : END IF
2152 4 : IF (ASSOCIATED(norm_drhob)) THEN
2153 : CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [deriv_norm_drhob], bo, &
2154 4 : weight, drb1drb, v_laplaceb%array)
2155 : END IF
2156 :
2157 4 : IF (tau_f) THEN
2158 : CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [deriv_tau_a], bo, &
2159 4 : weight, tau1a, v_laplaceb%array)
2160 : CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [deriv_tau_b], bo, &
2161 4 : weight, tau1b, v_laplaceb%array)
2162 : END IF
2163 :
2164 : CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [deriv_laplace_rhoa], bo, &
2165 4 : weight, laplace1a, v_laplaceb%array)
2166 :
2167 : CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [deriv_laplace_rhob], bo, &
2168 4 : weight, laplace1b, v_laplaceb%array)
2169 : END IF
2170 : END DO
2171 :
2172 10 : CALL virial_drho_drho(virial_pw, drhoa, v_drhoa, virial_xc)
2173 10 : CALL virial_drho_drho(virial_pw, drhob, v_drhob, virial_xc)
2174 10 : CALL virial_drho_drho(virial_pw, drho, v_drho, virial_xc)
2175 :
2176 10 : CALL deallocate_pw(v_drho, pw_pool)
2177 10 : CALL deallocate_pw(v_drhoa, pw_pool)
2178 10 : CALL deallocate_pw(v_drhob, pw_pool)
2179 :
2180 10 : IF (laplace_f) THEN
2181 15026 : virial_pw%array(:, :, :) = -rhoa(:, :, :)
2182 2 : CALL virial_laplace(virial_pw, pw_pool, virial_xc, v_laplacea%array)
2183 2 : CALL deallocate_pw(v_laplacea, pw_pool)
2184 :
2185 15026 : virial_pw%array(:, :, :) = -rhob(:, :, :)
2186 2 : CALL virial_laplace(virial_pw, pw_pool, virial_xc, v_laplaceb%array)
2187 2 : CALL deallocate_pw(v_laplaceb, pw_pool)
2188 : END IF
2189 :
2190 10 : CALL deallocate_pw(virial_pw, pw_pool)
2191 :
2192 40 : DO idir = 1, 3
2193 30 : DEALLOCATE (drho(idir)%array)
2194 40 : DEALLOCATE (drho1(idir)%array)
2195 : END DO
2196 10 : DEALLOCATE (dra1dra, drb1drb)
2197 :
2198 : ELSE
2199 26 : CALL xc_rho_set_get(rho1_set, rho=rho1, tau=tau1, laplace_rho=laplace1, can_return_null=.TRUE.)
2200 26 : CALL xc_rho_set_get(rho2_set, norm_drho=norm_drho2, laplace_rho=laplace2, can_return_null=.TRUE.)
2201 :
2202 200 : DO istep = -nsteps, nsteps
2203 174 : IF (istep == 0) CYCLE
2204 148 : weight = weights(istep, nsteps)/h
2205 148 : step = REAL(istep, dp)*h
2206 148 : CALL get_derivs_rho(norm_drho2, norm_drho, step, xc_fun_section, lsd, rho2_set, deriv_set1)
2207 :
2208 : ! Obtain the numerical 2nd derivatives w.r.t. to drho and collect the potential
2209 : CALL update_deriv_rho(deriv_set1, [deriv_rho], bo, &
2210 148 : norm_drho, gradient_cut, weight, rho1, v_drho%array)
2211 : CALL update_deriv_rho(deriv_set1, [deriv_norm_drho], bo, &
2212 148 : norm_drho, gradient_cut, weight, dr1dr, v_drho%array)
2213 :
2214 148 : IF (tau_f) THEN
2215 : CALL update_deriv_rho(deriv_set1, [deriv_tau], bo, &
2216 24 : norm_drho, gradient_cut, weight, tau1, v_drho%array)
2217 : END IF
2218 174 : IF (laplace_f) THEN
2219 : CALL update_deriv_rho(deriv_set1, [deriv_laplace_rho], bo, &
2220 12 : norm_drho, gradient_cut, weight, laplace1, v_drho%array)
2221 :
2222 12 : CALL get_derivs_rho(laplace2, laplace, step, xc_fun_section, lsd, rho2_set, deriv_set1)
2223 :
2224 : ! Obtain the numerical 2nd derivatives w.r.t. to drho and collect the potential
2225 : CALL update_deriv(deriv_set1, laplace, rho_cutoff, [deriv_rho], bo, &
2226 12 : weight, rho1, v_laplace%array)
2227 : CALL update_deriv(deriv_set1, laplace, rho_cutoff, [deriv_norm_drho], bo, &
2228 12 : weight, dr1dr, v_laplace%array)
2229 :
2230 12 : IF (tau_f) THEN
2231 : CALL update_deriv(deriv_set1, laplace, rho_cutoff, [deriv_tau], bo, &
2232 12 : weight, tau1, v_laplace%array)
2233 : END IF
2234 :
2235 : CALL update_deriv(deriv_set1, laplace, rho_cutoff, [deriv_laplace_rho], bo, &
2236 12 : weight, laplace1, v_laplace%array)
2237 : END IF
2238 : END DO
2239 :
2240 : ! Calculate the virial contribution from the potential
2241 26 : CALL virial_drho_drho(virial_pw, drho, v_drho, virial_xc)
2242 :
2243 26 : CALL deallocate_pw(v_drho, pw_pool)
2244 :
2245 26 : IF (laplace_f) THEN
2246 28862 : virial_pw%array(:, :, :) = -rho(:, :, :)
2247 2 : CALL virial_laplace(virial_pw, pw_pool, virial_xc, v_laplace%array)
2248 2 : CALL deallocate_pw(v_laplace, pw_pool)
2249 : END IF
2250 :
2251 26 : CALL deallocate_pw(virial_pw, pw_pool)
2252 : END IF
2253 :
2254 : END IF
2255 :
2256 36 : CALL xc_dset_release(deriv_set1)
2257 :
2258 36 : DEALLOCATE (dr1dr)
2259 :
2260 36 : CALL xc_rho_set_release(rho1_set)
2261 36 : CALL xc_rho_set_release(rho2_set)
2262 : END IF
2263 :
2264 848 : DO ispin = 1, SIZE(rho_r)
2265 848 : CALL pw_pool%give_back_pw(rho_r(ispin))
2266 : END DO
2267 340 : DEALLOCATE (rho_r)
2268 :
2269 340 : IF (ASSOCIATED(tau_r)) THEN
2270 254 : DO ispin = 1, SIZE(tau_r)
2271 254 : CALL pw_pool%give_back_pw(tau_r(ispin))
2272 : END DO
2273 100 : DEALLOCATE (tau_r)
2274 : END IF
2275 :
2276 340 : CALL timestop(handle)
2277 :
2278 15300 : END SUBROUTINE xc_calc_2nd_deriv_numerical
2279 :
2280 : ! **************************************************************************************************
2281 : !> \brief ...
2282 : !> \param rho_r ...
2283 : !> \param rho_g ...
2284 : !> \param rho1_r ...
2285 : !> \param rhoa ...
2286 : !> \param rhob ...
2287 : !> \param vxc_rho ...
2288 : !> \param tau_r ...
2289 : !> \param tau1_r ...
2290 : !> \param tau_a ...
2291 : !> \param tau_b ...
2292 : !> \param vxc_tau ...
2293 : !> \param xc_section ...
2294 : !> \param pw_pool ...
2295 : !> \param step ...
2296 : ! **************************************************************************************************
2297 792 : SUBROUTINE calc_resp_potential_numer_ab(rho_r, rho_g, rho1_r, rhoa, rhob, vxc_rho, &
2298 : tau_r, tau1_r, tau_a, tau_b, vxc_tau, &
2299 : xc_section, pw_pool, step)
2300 :
2301 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER, INTENT(IN) :: vxc_rho, vxc_tau
2302 : TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(IN) :: rho1_r
2303 : TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(IN), POINTER :: tau1_r
2304 : TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool
2305 : TYPE(section_vals_type), INTENT(IN), POINTER :: xc_section
2306 : REAL(KIND=dp), INTENT(IN) :: step
2307 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER, INTENT(IN) :: rhoa, rhob, tau_a, tau_b
2308 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER, INTENT(IN) :: rho_r
2309 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
2310 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: tau_r
2311 :
2312 : CHARACTER(len=*), PARAMETER :: routineN = 'calc_resp_potential_numer_ab'
2313 :
2314 : INTEGER :: handle
2315 : REAL(KIND=dp) :: exc
2316 : REAL(KIND=dp), DIMENSION(3, 3) :: virial_dummy
2317 :
2318 792 : CALL timeset(routineN, handle)
2319 :
2320 792 : !$OMP PARALLEL DEFAULT(NONE) SHARED(rho_r,rhoa,rhob,step,rho1_r,tau_r,tau_a,tau_b,tau1_r)
2321 : !$OMP WORKSHARE
2322 : rho_r(1)%array(:, :, :) = rhoa(:, :, :) + step*rho1_r(1)%array(:, :, :)
2323 : !$OMP END WORKSHARE NOWAIT
2324 : !$OMP WORKSHARE
2325 : rho_r(2)%array(:, :, :) = rhob(:, :, :) + step*rho1_r(2)%array(:, :, :)
2326 : !$OMP END WORKSHARE NOWAIT
2327 : IF (ASSOCIATED(tau1_r) .AND. ASSOCIATED(tau_r) .AND. ASSOCIATED(tau_a) .AND. ASSOCIATED(tau_b)) THEN
2328 : !$OMP WORKSHARE
2329 : tau_r(1)%array(:, :, :) = tau_a(:, :, :) + step*tau1_r(1)%array(:, :, :)
2330 : !$OMP END WORKSHARE NOWAIT
2331 : !$OMP WORKSHARE
2332 : tau_r(2)%array(:, :, :) = tau_b(:, :, :) + step*tau1_r(2)%array(:, :, :)
2333 : !$OMP END WORKSHARE NOWAIT
2334 : END IF
2335 : !$OMP END PARALLEL
2336 : CALL xc_vxc_pw_create(vxc_rho, vxc_tau, exc, rho_r, rho_g, tau_r, xc_section, &
2337 792 : pw_pool, .FALSE., virial_dummy)
2338 :
2339 792 : CALL timestop(handle)
2340 :
2341 792 : END SUBROUTINE calc_resp_potential_numer_ab
2342 :
2343 : ! **************************************************************************************************
2344 : !> \brief calculates stress tensor and potential contributions from the first derivative
2345 : !> \param deriv_set ...
2346 : !> \param description ...
2347 : !> \param virial_pw ...
2348 : !> \param drho ...
2349 : !> \param drho1 ...
2350 : !> \param virial_xc ...
2351 : !> \param norm_drho ...
2352 : !> \param gradient_cut ...
2353 : !> \param dr1dr ...
2354 : !> \param v_drho ...
2355 : ! **************************************************************************************************
2356 52 : SUBROUTINE apply_drho(deriv_set, description, virial_pw, drho, drho1, virial_xc, norm_drho, gradient_cut, dr1dr, v_drho)
2357 :
2358 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
2359 : INTEGER, DIMENSION(:), INTENT(in) :: description
2360 : TYPE(pw_r3d_rs_type), INTENT(IN) :: virial_pw
2361 : TYPE(cp_3d_r_cp_type), DIMENSION(3), INTENT(IN) :: drho, drho1
2362 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT) :: virial_xc
2363 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: norm_drho
2364 : REAL(KIND=dp), INTENT(IN) :: gradient_cut
2365 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: dr1dr
2366 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: v_drho
2367 :
2368 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_drho'
2369 :
2370 : INTEGER :: handle
2371 52 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: deriv_data
2372 : TYPE(xc_derivative_type), POINTER :: deriv_att
2373 :
2374 52 : CALL timeset(routineN, handle)
2375 :
2376 52 : deriv_att => xc_dset_get_derivative(deriv_set, description)
2377 52 : IF (ASSOCIATED(deriv_att)) THEN
2378 52 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
2379 52 : CALL virial_drho_drho1(virial_pw, drho, drho1, deriv_data, virial_xc)
2380 :
2381 52 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(dr1dr,gradient_cut,norm_drho,v_drho,deriv_data)
2382 : v_drho(:, :, :) = v_drho(:, :, :) + &
2383 : deriv_data(:, :, :)*dr1dr(:, :, :)/MAX(gradient_cut, norm_drho(:, :, :))**2
2384 : !$OMP END PARALLEL WORKSHARE
2385 : END IF
2386 :
2387 52 : CALL timestop(handle)
2388 :
2389 52 : END SUBROUTINE apply_drho
2390 :
2391 : ! **************************************************************************************************
2392 : !> \brief adds potential contributions from derivatives of rho or diagonal terms of norm_drho
2393 : !> \param deriv_set1 ...
2394 : !> \param description ...
2395 : !> \param bo ...
2396 : !> \param norm_drho norm_drho of which derivative is calculated
2397 : !> \param gradient_cut ...
2398 : !> \param h ...
2399 : !> \param rho1 function to contract the derivative with (rho1 for rho, dr1dr for norm_drho)
2400 : !> \param v_drho ...
2401 : ! **************************************************************************************************
2402 728 : SUBROUTINE update_deriv_rho(deriv_set1, description, bo, norm_drho, gradient_cut, weight, rho1, v_drho)
2403 :
2404 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set1
2405 : INTEGER, DIMENSION(:), INTENT(in) :: description
2406 : INTEGER, DIMENSION(2, 3), INTENT(IN) :: bo
2407 : REAL(KIND=dp), DIMENSION(bo(1, 1):bo(2, 1), bo(1, &
2408 : 2):bo(2, 2), bo(1, 3):bo(2, 3)), INTENT(IN) :: norm_drho
2409 : REAL(KIND=dp), INTENT(IN) :: gradient_cut, weight
2410 : REAL(KIND=dp), DIMENSION(bo(1, 1):bo(2, 1), bo(1, &
2411 : 2):bo(2, 2), bo(1, 3):bo(2, 3)), INTENT(IN) :: rho1
2412 : REAL(KIND=dp), DIMENSION(bo(1, 1):bo(2, 1), bo(1, &
2413 : 2):bo(2, 2), bo(1, 3):bo(2, 3)), INTENT(INOUT) :: v_drho
2414 :
2415 : CHARACTER(len=*), PARAMETER :: routineN = 'update_deriv_rho'
2416 :
2417 : INTEGER :: handle, i, j, k
2418 : REAL(KIND=dp) :: de
2419 728 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: deriv_data1
2420 : TYPE(xc_derivative_type), POINTER :: deriv_att1
2421 :
2422 728 : CALL timeset(routineN, handle)
2423 :
2424 : ! Obtain the numerical 2nd derivatives w.r.t. to drho and collect the potential
2425 728 : deriv_att1 => xc_dset_get_derivative(deriv_set1, description)
2426 728 : IF (ASSOCIATED(deriv_att1)) THEN
2427 728 : CALL xc_derivative_get(deriv_att1, deriv_data=deriv_data1)
2428 : !$OMP PARALLEL DO DEFAULT(NONE) &
2429 : !$OMP SHARED(bo,deriv_data1,weight,norm_drho,v_drho,rho1,gradient_cut) &
2430 : !$OMP PRIVATE(i,j,k,de) &
2431 728 : !$OMP COLLAPSE(3)
2432 : DO k = bo(1, 3), bo(2, 3)
2433 : DO j = bo(1, 2), bo(2, 2)
2434 : DO i = bo(1, 1), bo(2, 1)
2435 : de = weight*deriv_data1(i, j, k)/MAX(gradient_cut, norm_drho(i, j, k))**2
2436 : v_drho(i, j, k) = v_drho(i, j, k) - de*rho1(i, j, k)
2437 : END DO
2438 : END DO
2439 : END DO
2440 : !$OMP END PARALLEL DO
2441 : END IF
2442 :
2443 728 : CALL timestop(handle)
2444 :
2445 728 : END SUBROUTINE update_deriv_rho
2446 :
2447 : ! **************************************************************************************************
2448 : !> \brief adds potential contributions from derivatives of a component with positive and negative values
2449 : !> \param deriv_set1 ...
2450 : !> \param description ...
2451 : !> \param bo ...
2452 : !> \param h ...
2453 : !> \param rho1 function to contract the derivative with (rho1 for rho, dr1dr for norm_drho)
2454 : !> \param v ...
2455 : ! **************************************************************************************************
2456 120 : SUBROUTINE update_deriv(deriv_set1, rho, rho_cutoff, description, bo, weight, rho1, v)
2457 :
2458 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set1
2459 : INTEGER, DIMENSION(:), INTENT(in) :: description
2460 : INTEGER, DIMENSION(2, 3), INTENT(IN) :: bo
2461 : REAL(KIND=dp), INTENT(IN) :: weight, rho_cutoff
2462 : REAL(KIND=dp), DIMENSION(bo(1, 1):bo(2, 1), bo(1, &
2463 : 2):bo(2, 2), bo(1, 3):bo(2, 3)), INTENT(IN) :: rho, rho1
2464 : REAL(KIND=dp), DIMENSION(bo(1, 1):bo(2, 1), bo(1, &
2465 : 2):bo(2, 2), bo(1, 3):bo(2, 3)), INTENT(INOUT) :: v
2466 :
2467 : CHARACTER(len=*), PARAMETER :: routineN = 'update_deriv'
2468 :
2469 : INTEGER :: handle, i, j, k
2470 : REAL(KIND=dp) :: de
2471 120 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: deriv_data1
2472 : TYPE(xc_derivative_type), POINTER :: deriv_att1
2473 :
2474 120 : CALL timeset(routineN, handle)
2475 :
2476 : ! Obtain the numerical 2nd derivatives w.r.t. to drho and collect the potential
2477 120 : deriv_att1 => xc_dset_get_derivative(deriv_set1, description)
2478 120 : IF (ASSOCIATED(deriv_att1)) THEN
2479 120 : CALL xc_derivative_get(deriv_att1, deriv_data=deriv_data1)
2480 : !$OMP PARALLEL DO DEFAULT(NONE) &
2481 : !$OMP SHARED(bo,deriv_data1,weight,v,rho1,rho, rho_cutoff) &
2482 : !$OMP PRIVATE(i,j,k,de) &
2483 120 : !$OMP COLLAPSE(3)
2484 : DO k = bo(1, 3), bo(2, 3)
2485 : DO j = bo(1, 2), bo(2, 2)
2486 : DO i = bo(1, 1), bo(2, 1)
2487 : ! We have to consider that the given density (mostly the Laplacian) may have positive and negative values
2488 : de = weight*deriv_data1(i, j, k)/SIGN(MAX(ABS(rho(i, j, k)), rho_cutoff), rho(i, j, k))
2489 : v(i, j, k) = v(i, j, k) + de*rho1(i, j, k)
2490 : END DO
2491 : END DO
2492 : END DO
2493 : !$OMP END PARALLEL DO
2494 : END IF
2495 :
2496 120 : CALL timestop(handle)
2497 :
2498 120 : END SUBROUTINE update_deriv
2499 :
2500 : ! **************************************************************************************************
2501 : !> \brief adds mixed derivatives of norm_drho
2502 : !> \param deriv_set1 ...
2503 : !> \param description ...
2504 : !> \param bo ...
2505 : !> \param norm_drhoa norm_drho of which derivatives is calculated
2506 : !> \param gradient_cut ...
2507 : !> \param h ...
2508 : !> \param dra1dra dr1dr corresponding to norm_drho
2509 : !> \param drb1drb ...
2510 : !> \param v_drhoa potential corresponding to norm_drho
2511 : !> \param v_drhob ...
2512 : ! **************************************************************************************************
2513 216 : SUBROUTINE update_deriv_drho_ab(deriv_set1, description, bo, &
2514 216 : norm_drhoa, gradient_cut, weight, dra1dra, drb1drb, v_drhoa, v_drhob)
2515 :
2516 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set1
2517 : INTEGER, DIMENSION(:), INTENT(in) :: description
2518 : INTEGER, DIMENSION(2, 3), INTENT(IN) :: bo
2519 : REAL(KIND=dp), DIMENSION(bo(1, 1):bo(2, 1), bo(1, &
2520 : 2):bo(2, 2), bo(1, 3):bo(2, 3)), INTENT(IN) :: norm_drhoa
2521 : REAL(KIND=dp), INTENT(IN) :: gradient_cut, weight
2522 : REAL(KIND=dp), DIMENSION(bo(1, 1):bo(2, 1), bo(1, &
2523 : 2):bo(2, 2), bo(1, 3):bo(2, 3)), INTENT(IN) :: dra1dra, drb1drb
2524 : REAL(KIND=dp), DIMENSION(bo(1, 1):bo(2, 1), bo(1, &
2525 : 2):bo(2, 2), bo(1, 3):bo(2, 3)), INTENT(INOUT) :: v_drhoa, v_drhob
2526 :
2527 : CHARACTER(len=*), PARAMETER :: routineN = 'update_deriv_drho_ab'
2528 :
2529 : INTEGER :: handle, i, j, k
2530 : REAL(KIND=dp) :: de
2531 216 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: deriv_data1
2532 : TYPE(xc_derivative_type), POINTER :: deriv_att1
2533 :
2534 216 : CALL timeset(routineN, handle)
2535 :
2536 216 : deriv_att1 => xc_dset_get_derivative(deriv_set1, description)
2537 216 : IF (ASSOCIATED(deriv_att1)) THEN
2538 168 : CALL xc_derivative_get(deriv_att1, deriv_data=deriv_data1)
2539 : !$OMP PARALLEL DO DEFAULT(NONE) &
2540 : !$OMP PRIVATE(k,j,i,de) &
2541 : !$OMP SHARED(bo,drb1drb,dra1dra,deriv_data1,weight,gradient_cut,norm_drhoa,v_drhoa,v_drhob) &
2542 168 : !$OMP COLLAPSE(3)
2543 : DO k = bo(1, 3), bo(2, 3)
2544 : DO j = bo(1, 2), bo(2, 2)
2545 : DO i = bo(1, 1), bo(2, 1)
2546 : ! We introduce a factor of two because we will average between both numerical derivatives
2547 : de = 0.5_dp*weight*deriv_data1(i, j, k)/MAX(gradient_cut, norm_drhoa(i, j, k))**2
2548 : v_drhoa(i, j, k) = v_drhoa(i, j, k) - de*drb1drb(i, j, k)
2549 : v_drhob(i, j, k) = v_drhob(i, j, k) - de*dra1dra(i, j, k)
2550 : END DO
2551 : END DO
2552 : END DO
2553 : !$OMP END PARALLEL DO
2554 : END IF
2555 :
2556 216 : CALL timestop(handle)
2557 :
2558 216 : END SUBROUTINE update_deriv_drho_ab
2559 :
2560 : ! **************************************************************************************************
2561 : !> \brief calculate derivative sets for helper points
2562 : !> \param norm_drho2 norm_drho of new points
2563 : !> \param norm_drho norm_drho of KS density
2564 : !> \param h ...
2565 : !> \param xc_fun_section ...
2566 : !> \param lsd ...
2567 : !> \param rho2_set rho_set for new points
2568 : !> \param deriv_set1 will contain derivatives of the perturbed density
2569 : ! **************************************************************************************************
2570 276 : SUBROUTINE get_derivs_rho(norm_drho2, norm_drho, step, xc_fun_section, lsd, rho2_set, deriv_set1)
2571 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: norm_drho2
2572 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: norm_drho
2573 : REAL(KIND=dp), INTENT(IN) :: step
2574 : TYPE(section_vals_type), INTENT(IN), POINTER :: xc_fun_section
2575 : LOGICAL, INTENT(IN) :: lsd
2576 : TYPE(xc_rho_set_type), INTENT(INOUT) :: rho2_set
2577 : TYPE(xc_derivative_set_type) :: deriv_set1
2578 :
2579 : CHARACTER(len=*), PARAMETER :: routineN = 'get_derivs_rho'
2580 :
2581 : INTEGER :: handle
2582 :
2583 276 : CALL timeset(routineN, handle)
2584 :
2585 : ! Copy the densities, do one step into the direction of drho
2586 276 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(norm_drho,norm_drho2,step)
2587 : norm_drho2 = norm_drho*(1.0_dp + step)
2588 : !$OMP END PARALLEL WORKSHARE
2589 :
2590 276 : CALL xc_dset_zero_all(deriv_set1)
2591 :
2592 : ! Calculate the derivatives of the functional
2593 : CALL xc_functionals_eval(xc_fun_section, &
2594 : lsd=lsd, &
2595 : rho_set=rho2_set, &
2596 : deriv_set=deriv_set1, &
2597 276 : deriv_order=1)
2598 :
2599 : ! Return to the original values
2600 276 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(norm_drho,norm_drho2)
2601 : norm_drho2 = norm_drho
2602 : !$OMP END PARALLEL WORKSHARE
2603 :
2604 276 : CALL divide_by_norm_drho(deriv_set1, rho2_set, lsd)
2605 :
2606 276 : CALL timestop(handle)
2607 :
2608 276 : END SUBROUTINE get_derivs_rho
2609 :
2610 : ! **************************************************************************************************
2611 : !> \brief Calculates the second derivative of E_xc at rho in the direction
2612 : !> rho1 (if you see the second derivative as bilinear form)
2613 : !> partial_rho|_(rho=rho) partial_rho|_(rho=rho) E_xc drho(rho1)drho
2614 : !> The other direction is still undetermined, thus it returns
2615 : !> a potential (partial integration is performed to reduce it to
2616 : !> function of rho, removing the dependence from its partial derivs)
2617 : !> Has to be called after the setup by xc_prep_2nd_deriv.
2618 : !> \param v_xc exchange-correlation potential
2619 : !> \param v_xc_tau ...
2620 : !> \param deriv_set derivatives of the exchange-correlation potential
2621 : !> \param rho_set object containing the density at which the derivatives were calculated
2622 : !> \param rho1_set object containing the density with which to fold
2623 : !> \param pw_pool the pool for the grids
2624 : !> \param xc_section XC parameters
2625 : !> \param gapw Gaussian and augmented plane waves calculation
2626 : !> \param vxg ...
2627 : !> \param tddfpt_fac factor that multiplies the crossterms (tddfpt triplets
2628 : !> on a closed shell system it should be -1, defaults to 1)
2629 : !> \param compute_virial ...
2630 : !> \param virial_xc ...
2631 : !> \note
2632 : !> The old version of this routine was smarter: it handled split_desc(1)
2633 : !> and split_desc(2) separately, thus the code automatically handled all
2634 : !> possible cross terms (you only had to check if it was diagonal to avoid
2635 : !> double counting). I think that is the way to go if you want to add more
2636 : !> terms (tau,rho in LSD,...). The problem with the old code was that it
2637 : !> because of the old functional structure it sometime guessed wrongly
2638 : !> which derivative was where. There were probably still bugs with gradient
2639 : !> corrected functionals (never tested), and it didn't contain first
2640 : !> derivatives with respect to drho (that contribute also to the second
2641 : !> derivative wrt. rho).
2642 : !> The code was a little complex because it really tried to handle any
2643 : !> functional derivative in the most efficient way with the given contents of
2644 : !> rho_set.
2645 : !> Anyway I strongly encourage whoever wants to modify this code to give a
2646 : !> look to the old version. [fawzi]
2647 : ! **************************************************************************************************
2648 33196 : SUBROUTINE xc_calc_2nd_deriv_analytical(v_xc, v_xc_tau, deriv_set, rho_set, rho1_set, &
2649 : pw_pool, xc_section, gapw, vxg, tddfpt_fac, &
2650 : compute_virial, virial_xc, spinflip)
2651 :
2652 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_xc, v_xc_tau
2653 : TYPE(xc_derivative_set_type) :: deriv_set
2654 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set, rho1_set
2655 : TYPE(pw_pool_type), POINTER :: pw_pool
2656 : TYPE(section_vals_type), POINTER :: xc_section
2657 : LOGICAL, INTENT(IN), OPTIONAL :: gapw
2658 : REAL(kind=dp), DIMENSION(:, :, :, :), OPTIONAL, &
2659 : POINTER :: vxg
2660 : REAL(kind=dp), INTENT(in), OPTIONAL :: tddfpt_fac
2661 : LOGICAL, INTENT(IN), OPTIONAL :: compute_virial
2662 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT), &
2663 : OPTIONAL :: virial_xc
2664 : LOGICAL, INTENT(in), OPTIONAL :: spinflip
2665 :
2666 : CHARACTER(len=*), PARAMETER :: routineN = 'xc_calc_2nd_deriv_analytical'
2667 :
2668 : INTEGER :: handle, i, ia, idir, ir, ispin, j, jdir, &
2669 : k, nspins, xc_deriv_method_id
2670 : INTEGER, DIMENSION(2, 3) :: bo
2671 : LOGICAL :: gradient_f, lsd, my_compute_virial, alda0, &
2672 : my_gapw, tau_f, laplace_f, rho_f, do_spinflip
2673 : REAL(KIND=dp) :: fac, gradient_cut, tmp, factor2, s, S_THRESH
2674 33196 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: dr1dr, dra1dra, drb1drb
2675 33196 : REAL(kind=dp), DIMENSION(:, :, :), POINTER :: deriv_data, deriv_data2, &
2676 33196 : e_drhoa, e_drhob, e_drho, norm_drho, norm_drhoa, &
2677 33196 : norm_drhob, rho1, rho1a, rho1b, &
2678 33196 : tau1, tau1a, tau1b, laplace1, laplace1a, laplace1b, &
2679 33196 : rho, rhoa, rhob
2680 630724 : TYPE(cp_3d_r_cp_type), DIMENSION(3) :: drho, drho1, drho1a, drho1b, drhoa, drhob
2681 33196 : TYPE(pw_r3d_rs_type), DIMENSION(:), ALLOCATABLE :: v_drhoa, v_drhob, v_drho, v_laplace
2682 33196 : TYPE(pw_r3d_rs_type), DIMENSION(:, :), ALLOCATABLE :: v_drho_r
2683 : TYPE(pw_r3d_rs_type) :: virial_pw
2684 : TYPE(pw_c1d_gs_type) :: tmp_g, vxc_g
2685 : TYPE(xc_derivative_type), POINTER :: deriv_att
2686 :
2687 33196 : CALL timeset(routineN, handle)
2688 :
2689 33196 : NULLIFY (e_drhoa, e_drhob, e_drho)
2690 :
2691 33196 : my_gapw = .FALSE.
2692 33196 : IF (PRESENT(gapw)) my_gapw = gapw
2693 :
2694 33196 : my_compute_virial = .FALSE.
2695 33196 : IF (PRESENT(compute_virial)) my_compute_virial = compute_virial
2696 :
2697 33196 : CPASSERT(ASSOCIATED(v_xc))
2698 33196 : CPASSERT(ASSOCIATED(xc_section))
2699 33196 : IF (my_gapw) THEN
2700 11166 : CPASSERT(PRESENT(vxg))
2701 : END IF
2702 33196 : IF (my_compute_virial) THEN
2703 358 : CPASSERT(PRESENT(virial_xc))
2704 : END IF
2705 :
2706 : CALL section_vals_val_get(xc_section, "XC_GRID%XC_DERIV", &
2707 33196 : i_val=xc_deriv_method_id)
2708 33196 : CALL xc_rho_set_get(rho_set, drho_cutoff=gradient_cut)
2709 33196 : nspins = SIZE(v_xc)
2710 33196 : lsd = ASSOCIATED(rho_set%rhoa)
2711 33196 : fac = 0.0_dp
2712 33196 : factor2 = 1.0_dp
2713 33196 : IF (PRESENT(tddfpt_fac)) fac = tddfpt_fac
2714 33196 : IF (PRESENT(tddfpt_fac)) factor2 = tddfpt_fac
2715 33196 : do_spinflip = .FALSE.
2716 33196 : IF (PRESENT(spinflip)) do_spinflip = spinflip
2717 :
2718 331960 : bo = rho_set%local_bounds
2719 :
2720 33196 : CALL check_for_derivatives(deriv_set, lsd, rho_f, gradient_f, tau_f, laplace_f)
2721 :
2722 33196 : alda0 = .false.
2723 33196 : IF (gradient_f) THEN
2724 : S_THRESH = 1.0E-04
2725 : ELSE
2726 12278 : S_THRESH = 1.0E-10
2727 : END IF
2728 :
2729 33196 : IF (tau_f) THEN
2730 656 : CPASSERT(ASSOCIATED(v_xc_tau))
2731 : END IF
2732 :
2733 33196 : IF (gradient_f) THEN
2734 219060 : ALLOCATE (v_drho_r(3, nspins), v_drho(nspins))
2735 43812 : DO ispin = 1, nspins
2736 91576 : DO idir = 1, 3
2737 91576 : CALL allocate_pw(v_drho_r(idir, ispin), pw_pool, bo)
2738 : END DO
2739 43812 : CALL allocate_pw(v_drho(ispin), pw_pool, bo)
2740 : END DO
2741 :
2742 20918 : IF (xc_requires_tmp_g(xc_deriv_method_id) .AND. .NOT. my_gapw) THEN
2743 12300 : IF (ASSOCIATED(pw_pool)) THEN
2744 12300 : CALL pw_pool%create_pw(tmp_g)
2745 12300 : CALL pw_pool%create_pw(vxc_g)
2746 : ELSE
2747 : ! remember to refix for gapw
2748 0 : CPABORT("XC_DERIV method is not implemented in GAPW")
2749 : END IF
2750 : END IF
2751 : END IF
2752 :
2753 70418 : DO ispin = 1, nspins
2754 1065080064 : v_xc(ispin)%array = 0.0_dp
2755 : END DO
2756 :
2757 33196 : IF (tau_f) THEN
2758 1386 : DO ispin = 1, nspins
2759 14897758 : v_xc_tau(ispin)%array = 0.0_dp
2760 : END DO
2761 : END IF
2762 :
2763 33196 : IF (laplace_f .AND. my_gapw) &
2764 0 : CPABORT("Laplace-dependent functional not implemented with GAPW!")
2765 :
2766 33196 : IF (my_compute_virial .AND. (gradient_f .OR. laplace_f)) CALL allocate_pw(virial_pw, pw_pool, bo)
2767 :
2768 33196 : IF (lsd) THEN
2769 :
2770 : !-------------------!
2771 : ! UNrestricted case !
2772 : !-------------------!
2773 :
2774 5168 : IF (do_spinflip) THEN
2775 414 : CALL xc_rho_set_get(rho1_set, rhoa=rho1a)
2776 414 : CALL xc_rho_set_get(rho_set, rhoa=rhoa, rhob=rhob)
2777 : ELSE
2778 4754 : CALL xc_rho_set_get(rho1_set, rhoa=rho1a, rhob=rho1b)
2779 : END IF
2780 :
2781 5168 : IF (gradient_f) THEN
2782 : CALL xc_rho_set_get(rho_set, drhoa=drhoa, drhob=drhob, &
2783 2738 : norm_drho=norm_drho, norm_drhoa=norm_drhoa, norm_drhob=norm_drhob)
2784 2738 : IF (do_spinflip) THEN
2785 262 : CALL xc_rho_set_get(rho1_set, drhoa=drho1a)
2786 262 : CALL calc_drho_from_a(drho1, drho1a)
2787 : ELSE
2788 2476 : CALL xc_rho_set_get(rho1_set, drhoa=drho1a, drhob=drho1b)
2789 2476 : CALL calc_drho_from_ab(drho1, drho1a, drho1b)
2790 : END IF
2791 :
2792 2738 : CALL calc_drho_from_ab(drho, drhoa, drhob)
2793 :
2794 2738 : CALL prepare_dr1dr(dra1dra, drhoa, drho1a)
2795 2738 : IF (do_spinflip) THEN
2796 262 : CALL prepare_dr1dr(drb1drb, drhob, drho1a)
2797 262 : CALL prepare_dr1dr(dr1dr, drho, drho1a)
2798 2476 : ELSE IF (nspins /= 1) THEN
2799 1766 : CALL prepare_dr1dr(drb1drb, drhob, drho1b)
2800 1766 : CALL prepare_dr1dr(dr1dr, drho, drho1)
2801 : ELSE
2802 710 : CALL prepare_dr1dr(drb1drb, drhob, drho1b)
2803 710 : CALL prepare_dr1dr_ab(dr1dr, drhoa, drhob, drho1a, drho1b, fac)
2804 : END IF
2805 :
2806 20380 : ALLOCATE (v_drhoa(nspins), v_drhob(nspins))
2807 7452 : DO ispin = 1, nspins
2808 4714 : CALL allocate_pw(v_drhoa(ispin), pw_pool, bo)
2809 7452 : CALL allocate_pw(v_drhob(ispin), pw_pool, bo)
2810 : END DO
2811 :
2812 : END IF
2813 :
2814 5168 : IF (laplace_f) THEN
2815 38 : CALL xc_rho_set_get(rho1_set, laplace_rhoa=laplace1a, laplace_rhob=laplace1b)
2816 :
2817 190 : ALLOCATE (v_laplace(nspins))
2818 114 : DO ispin = 1, nspins
2819 114 : CALL allocate_pw(v_laplace(ispin), pw_pool, bo)
2820 : END DO
2821 :
2822 38 : IF (my_compute_virial) CALL xc_rho_set_get(rho_set, rhoa=rhoa, rhob=rhob)
2823 : END IF
2824 :
2825 5168 : IF (tau_f) THEN
2826 74 : CALL xc_rho_set_get(rho1_set, tau_a=tau1a, tau_b=tau1b)
2827 : END IF
2828 :
2829 5168 : IF (do_spinflip) THEN
2830 :
2831 : ! vxc contributions
2832 : ! vxc = (vxc^{\alpha}-vxc^{\beta})*rho1/(rhoa-rhob)
2833 : ! Alpha LDA contribution
2834 : ! | d e_xc d e_xc | rho1a
2835 : ! vxca = |-------- - --------|*-------------
2836 : ! | drhoa drhob | |rhoa - rhob|
2837 414 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_rhoa])
2838 414 : IF (ASSOCIATED(deriv_att)) THEN
2839 414 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
2840 414 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_rhob])
2841 414 : IF (ASSOCIATED(deriv_att)) THEN
2842 414 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data2)
2843 : !$OMP PARALLEL DO PRIVATE(k,j,i,s) DEFAULT(NONE)&
2844 414 : !$OMP SHARED(bo,v_xc,deriv_data,deriv_data2,rho1a,rhoa,rhob,S_THRESH) COLLAPSE(3)
2845 : DO k = bo(1, 3), bo(2, 3)
2846 : DO j = bo(1, 2), bo(2, 2)
2847 : DO i = bo(1, 1), bo(2, 1)
2848 : s = MAX(ABS(rhoa(i, j, k) - rhob(i, j, k)), S_THRESH)
2849 : v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
2850 : (deriv_data(i, j, k) - deriv_data2(i, j, k))*rho1a(i, j, k)/s
2851 : END DO
2852 : END DO
2853 : END DO
2854 : !$OMP END PARALLEL DO
2855 : END IF
2856 : END IF
2857 : ! GGA contributions to the spin-flip xcKernel
2858 : ! GGA contribution
2859 : ! | d e_xc d e_xc | 1
2860 : ! vxca += |----------* dra1dra - ----------*drb1drb|*-------------
2861 : ! | d|drhoa| d|drhob| | |rhoa - rhob|
2862 : IF (.NOT. alda0) THEN
2863 414 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa])
2864 414 : IF (ASSOCIATED(deriv_att)) THEN
2865 262 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
2866 262 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob])
2867 262 : IF (ASSOCIATED(deriv_att)) THEN
2868 262 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data2)
2869 : !$OMP PARALLEL DO PRIVATE(k,j,i,s) DEFAULT(NONE)&
2870 262 : !$OMP SHARED(bo,deriv_data,deriv_data2,dra1dra,drb1drb,v_xc,rhoa,rhob,S_THRESH) COLLAPSE(3)
2871 : DO k = bo(1, 3), bo(2, 3)
2872 : DO j = bo(1, 2), bo(2, 2)
2873 : DO i = bo(1, 1), bo(2, 1)
2874 : s = MAX(ABS(rhoa(i, j, k) - rhob(i, j, k)), S_THRESH)
2875 : v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
2876 : (deriv_data(i, j, k)*dra1dra(i, j, k) - &
2877 : deriv_data2(i, j, k)*drb1drb(i, j, k))/s
2878 : END DO
2879 : END DO
2880 : END DO
2881 : !$OMP END PARALLEL DO
2882 : END IF
2883 : END IF
2884 : END IF
2885 :
2886 4754 : ELSE IF (nspins /= 1) THEN
2887 :
2888 : ! Compute \sum_{\tau}fxc^{\sigma\tau}*\rho^{\tau}(1) over the grid points
2889 33764 : $:add_2nd_derivative_terms(arguments_openshell)
2890 :
2891 : ELSE
2892 :
2893 : ! Compute (fxc^{\alpha\alpha}+-fxc^{\beta\beta})*\rho(1) over the grid points
2894 1038 : $:add_2nd_derivative_terms(arguments_triplet_outer, arguments_triplet_inner)
2895 :
2896 : END IF
2897 :
2898 5168 : IF (gradient_f) THEN
2899 2738 : IF (.NOT. do_spinflip) THEN
2900 :
2901 2476 : IF (my_compute_virial) THEN
2902 10 : CALL virial_drho_drho(virial_pw, drhoa, v_drhoa(1), virial_xc)
2903 10 : CALL virial_drho_drho(virial_pw, drhob, v_drhob(2), virial_xc)
2904 40 : DO idir = 1, 3
2905 30 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(drho,idir,v_drho,virial_pw)
2906 : virial_pw%array(:, :, :) = drho(idir)%array(:, :, :)*(v_drho(1)%array(:, :, :) + v_drho(2)%array(:, :, :))
2907 : !$OMP END PARALLEL WORKSHARE
2908 100 : DO jdir = 1, idir
2909 : tmp = -0.5_dp*virial_pw%pw_grid%dvol*accurate_dot_product(virial_pw%array(:, :, :), &
2910 60 : drho(jdir)%array(:, :, :))
2911 60 : virial_xc(jdir, idir) = virial_xc(jdir, idir) + tmp
2912 90 : virial_xc(idir, jdir) = virial_xc(jdir, idir)
2913 : END DO
2914 : END DO
2915 : END IF ! my_compute_virial
2916 :
2917 2476 : IF (my_gapw) THEN
2918 : !$OMP PARALLEL DO DEFAULT(NONE) &
2919 : !$OMP PRIVATE(ia,idir,ispin,ir) &
2920 : !$OMP SHARED(bo,nspins,vxg,drhoa,drhob,v_drhoa,v_drhob,v_drho, &
2921 764 : !$OMP e_drhoa,e_drhob,e_drho,drho1a,drho1b,fac,drho,drho1) COLLAPSE(3)
2922 : DO ir = bo(1, 2), bo(2, 2)
2923 : DO ia = bo(1, 1), bo(2, 1)
2924 : DO idir = 1, 3
2925 : DO ispin = 1, nspins
2926 : vxg(idir, ia, ir, ispin) = &
2927 : -(v_drhoa(ispin)%array(ia, ir, 1)*drhoa(idir)%array(ia, ir, 1) + &
2928 : v_drhob(ispin)%array(ia, ir, 1)*drhob(idir)%array(ia, ir, 1) + &
2929 : v_drho(ispin)%array(ia, ir, 1)*drho(idir)%array(ia, ir, 1))
2930 : END DO
2931 : IF (ASSOCIATED(e_drhoa)) THEN
2932 : vxg(idir, ia, ir, 1) = vxg(idir, ia, ir, 1) + &
2933 : e_drhoa(ia, ir, 1)*drho1a(idir)%array(ia, ir, 1)
2934 : END IF
2935 : IF (nspins /= 1 .AND. ASSOCIATED(e_drhob)) THEN
2936 : vxg(idir, ia, ir, 2) = vxg(idir, ia, ir, 2) + &
2937 : e_drhob(ia, ir, 1)*drho1b(idir)%array(ia, ir, 1)
2938 : END IF
2939 : IF (ASSOCIATED(e_drho)) THEN
2940 : IF (nspins /= 1) THEN
2941 : vxg(idir, ia, ir, 1) = vxg(idir, ia, ir, 1) + &
2942 : e_drho(ia, ir, 1)*drho1(idir)%array(ia, ir, 1)
2943 : vxg(idir, ia, ir, 2) = vxg(idir, ia, ir, 2) + &
2944 : e_drho(ia, ir, 1)*drho1(idir)%array(ia, ir, 1)
2945 : ELSE
2946 : vxg(idir, ia, ir, 1) = vxg(idir, ia, ir, 1) + &
2947 : e_drho(ia, ir, 1)*(drho1a(idir)%array(ia, ir, 1) + &
2948 : fac*drho1b(idir)%array(ia, ir, 1))
2949 : END IF
2950 : END IF
2951 : END DO
2952 : END DO
2953 : END DO
2954 : !$OMP END PARALLEL DO
2955 : ELSE
2956 :
2957 : ! partial integration
2958 6848 : DO idir = 1, 3
2959 :
2960 14466 : DO ispin = 1, nspins
2961 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) &
2962 14466 : !$OMP SHARED(v_drho_r,v_drhoa,v_drhob,v_drho,drhoa,drhob,drho,ispin,idir)
2963 : v_drho_r(idir, ispin)%array(:, :, :) = &
2964 : v_drhoa(ispin)%array(:, :, :)*drhoa(idir)%array(:, :, :) + &
2965 : v_drhob(ispin)%array(:, :, :)*drhob(idir)%array(:, :, :) + &
2966 : v_drho(ispin)%array(:, :, :)*drho(idir)%array(:, :, :)
2967 : !$OMP END PARALLEL WORKSHARE
2968 : END DO
2969 5136 : IF (ASSOCIATED(e_drhoa)) THEN
2970 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) &
2971 5136 : !$OMP SHARED(v_drho_r,e_drhoa,drho1a,idir)
2972 : v_drho_r(idir, 1)%array(:, :, :) = v_drho_r(idir, 1)%array(:, :, :) - &
2973 : e_drhoa(:, :, :)*drho1a(idir)%array(:, :, :)
2974 : !$OMP END PARALLEL WORKSHARE
2975 : END IF
2976 5136 : IF (nspins /= 1 .AND. ASSOCIATED(e_drhob)) THEN
2977 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE)&
2978 4194 : !$OMP SHARED(v_drho_r,e_drhob,drho1b,idir)
2979 : v_drho_r(idir, 2)%array(:, :, :) = v_drho_r(idir, 2)%array(:, :, :) - &
2980 : e_drhob(:, :, :)*drho1b(idir)%array(:, :, :)
2981 : !$OMP END PARALLEL WORKSHARE
2982 : END IF
2983 6848 : IF (ASSOCIATED(e_drho)) THEN
2984 : !$OMP PARALLEL DO PRIVATE(k,j,i) DEFAULT(NONE)&
2985 5136 : !$OMP SHARED(bo,v_drho_r,e_drho,drho1a,drho1b,drho1,fac,idir,nspins) COLLAPSE(3)
2986 : DO k = bo(1, 3), bo(2, 3)
2987 : DO j = bo(1, 2), bo(2, 2)
2988 : DO i = bo(1, 1), bo(2, 1)
2989 : IF (nspins /= 1) THEN
2990 : v_drho_r(idir, 1)%array(i, j, k) = v_drho_r(idir, 1)%array(i, j, k) - &
2991 : e_drho(i, j, k)*drho1(idir)%array(i, j, k)
2992 : v_drho_r(idir, 2)%array(i, j, k) = v_drho_r(idir, 2)%array(i, j, k) - &
2993 : e_drho(i, j, k)*drho1(idir)%array(i, j, k)
2994 : ELSE
2995 : v_drho_r(idir, 1)%array(i, j, k) = v_drho_r(idir, 1)%array(i, j, k) - &
2996 : e_drho(i, j, k)*(drho1a(idir)%array(i, j, k) + &
2997 : fac*drho1b(idir)%array(i, j, k))
2998 : END IF
2999 : END DO
3000 : END DO
3001 : END DO
3002 : !$OMP END PARALLEL DO
3003 : END IF
3004 : END DO
3005 :
3006 : ! partial integration
3007 4822 : DO ispin = 1, nspins
3008 4822 : CALL xc_pw_divergence(xc_deriv_method_id, v_drho_r(:, ispin), tmp_g, vxc_g, v_xc(ispin))
3009 : END DO ! ispin
3010 :
3011 : END IF
3012 :
3013 : END IF ! .NOT.do_spinflip
3014 :
3015 10952 : DO idir = 1, 3
3016 8214 : DEALLOCATE (drho(idir)%array)
3017 10952 : DEALLOCATE (drho1(idir)%array)
3018 : END DO
3019 :
3020 7452 : DO ispin = 1, nspins
3021 4714 : CALL deallocate_pw(v_drhoa(ispin), pw_pool)
3022 7452 : CALL deallocate_pw(v_drhob(ispin), pw_pool)
3023 : END DO
3024 :
3025 2738 : DEALLOCATE (v_drhoa, v_drhob)
3026 :
3027 : END IF ! gradient_f
3028 :
3029 5168 : IF (laplace_f .AND. my_compute_virial) THEN
3030 15026 : virial_pw%array(:, :, :) = -rhoa(:, :, :)
3031 2 : CALL virial_laplace(virial_pw, pw_pool, virial_xc, v_laplace(1)%array)
3032 15026 : virial_pw%array(:, :, :) = -rhob(:, :, :)
3033 2 : CALL virial_laplace(virial_pw, pw_pool, virial_xc, v_laplace(2)%array)
3034 : END IF
3035 :
3036 : ELSE
3037 :
3038 : !-----------------!
3039 : ! restricted case !
3040 : !-----------------!
3041 :
3042 28028 : CALL xc_rho_set_get(rho1_set, rho=rho1)
3043 :
3044 28028 : IF (gradient_f) THEN
3045 18180 : CALL xc_rho_set_get(rho_set, drho=drho, norm_drho=norm_drho)
3046 18180 : CALL xc_rho_set_get(rho1_set, drho=drho1)
3047 18180 : CALL prepare_dr1dr(dr1dr, drho, drho1)
3048 : END IF
3049 :
3050 28028 : IF (laplace_f) THEN
3051 136 : CALL xc_rho_set_get(rho1_set, laplace_rho=laplace1)
3052 :
3053 544 : ALLOCATE (v_laplace(nspins))
3054 272 : DO ispin = 1, nspins
3055 272 : CALL allocate_pw(v_laplace(ispin), pw_pool, bo)
3056 : END DO
3057 :
3058 136 : IF (my_compute_virial) CALL xc_rho_set_get(rho_set, rho=rho)
3059 : END IF
3060 :
3061 28028 : IF (tau_f) THEN
3062 582 : CALL xc_rho_set_get(rho1_set, tau=tau1)
3063 : END IF
3064 :
3065 322544 : $:add_2nd_derivative_terms(arguments_closedshell)
3066 :
3067 28028 : IF (gradient_f) THEN
3068 :
3069 18180 : IF (my_compute_virial) THEN
3070 222 : CALL virial_drho_drho(virial_pw, drho, v_drho(1), virial_xc)
3071 : END IF ! my_compute_virial
3072 :
3073 18180 : IF (my_gapw) THEN
3074 :
3075 29672 : DO idir = 1, 3
3076 : !$OMP PARALLEL DO DEFAULT(NONE) &
3077 : !$OMP PRIVATE(ia,ir) &
3078 : !$OMP SHARED(bo,vxg,drho,v_drho,e_drho,drho1,idir,factor2) &
3079 29672 : !$OMP COLLAPSE(2)
3080 : DO ia = bo(1, 1), bo(2, 1)
3081 : DO ir = bo(1, 2), bo(2, 2)
3082 : vxg(idir, ia, ir, 1) = -drho(idir)%array(ia, ir, 1)*v_drho(1)%array(ia, ir, 1)
3083 : IF (ASSOCIATED(e_drho)) THEN
3084 : vxg(idir, ia, ir, 1) = vxg(idir, ia, ir, 1) + factor2*drho1(idir)%array(ia, ir, 1)*e_drho(ia, ir, 1)
3085 : END IF
3086 : END DO
3087 : END DO
3088 : !$OMP END PARALLEL DO
3089 : END DO
3090 :
3091 : ELSE
3092 : ! partial integration
3093 43048 : DO idir = 1, 3
3094 43048 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(v_drho_r,drho,v_drho,drho1,e_drho,idir)
3095 : v_drho_r(idir, 1)%array(:, :, :) = drho(idir)%array(:, :, :)*v_drho(1)%array(:, :, :) - &
3096 : drho1(idir)%array(:, :, :)*e_drho(:, :, :)
3097 : !$OMP END PARALLEL WORKSHARE
3098 : END DO
3099 :
3100 10762 : CALL xc_pw_divergence(xc_deriv_method_id, v_drho_r(:, 1), tmp_g, vxc_g, v_xc(1))
3101 : END IF
3102 :
3103 : END IF
3104 :
3105 28028 : IF (laplace_f .AND. my_compute_virial) THEN
3106 294530 : virial_pw%array(:, :, :) = -rho(:, :, :)
3107 14 : CALL virial_laplace(virial_pw, pw_pool, virial_xc, v_laplace(1)%array)
3108 : END IF
3109 :
3110 : END IF
3111 :
3112 33196 : IF (laplace_f) THEN
3113 386 : DO ispin = 1, nspins
3114 212 : CALL xc_pw_laplace(v_laplace(ispin), pw_pool, xc_deriv_method_id)
3115 386 : CALL pw_axpy(v_laplace(ispin), v_xc(ispin))
3116 : END DO
3117 : END IF
3118 :
3119 33196 : IF (gradient_f) THEN
3120 :
3121 43812 : DO ispin = 1, nspins
3122 22894 : CALL deallocate_pw(v_drho(ispin), pw_pool)
3123 112494 : DO idir = 1, 3
3124 91576 : CALL deallocate_pw(v_drho_r(idir, ispin), pw_pool)
3125 : END DO
3126 : END DO
3127 20918 : DEALLOCATE (v_drho, v_drho_r)
3128 :
3129 : END IF
3130 :
3131 33196 : IF (laplace_f) THEN
3132 386 : DO ispin = 1, nspins
3133 386 : CALL deallocate_pw(v_laplace(ispin), pw_pool)
3134 : END DO
3135 174 : DEALLOCATE (v_laplace)
3136 : END IF
3137 :
3138 33196 : IF (ASSOCIATED(tmp_g%pw_grid) .AND. ASSOCIATED(pw_pool)) THEN
3139 12300 : CALL pw_pool%give_back_pw(tmp_g)
3140 : END IF
3141 :
3142 33196 : IF (ASSOCIATED(vxc_g%pw_grid) .AND. ASSOCIATED(pw_pool)) THEN
3143 12300 : CALL pw_pool%give_back_pw(vxc_g)
3144 : END IF
3145 :
3146 33196 : IF (my_compute_virial .AND. (gradient_f .OR. laplace_f)) THEN
3147 232 : CALL deallocate_pw(virial_pw, pw_pool)
3148 : END IF
3149 :
3150 33196 : CALL timestop(handle)
3151 :
3152 99588 : END SUBROUTINE xc_calc_2nd_deriv_analytical
3153 :
3154 : ! **************************************************************************************************
3155 : !> \brief Calculates the third functional derivative of the exchange-correlation functional, E_xc.
3156 : !> Any GGA functional can be written as:
3157 : !>
3158 : !> E_xc[\rho] = \int e_xc(\rho,\nabla\rho)dr
3159 : !>
3160 : !> This routine gives you back the contraction of the derivatives of e_xc with respect to the
3161 : !> alpha or beta density or with respect to the norm of their gradients contracted with rho1.
3162 : !> For example, the alpha component would be (d stands for total derivative):
3163 : !>
3164 : !> d^3 e_xc
3165 : !> v_xc(1) = \sum_{s,s'}^{a,b} ---------------------\rhos1\rho1s'
3166 : !> d\rhoa d\rhos d\rhos'
3167 : !>
3168 : !> \param v_xc Third derivative of the exchange-correlation functional
3169 : !> \param v_xc_tau ...
3170 : !> \param deriv_set derivatives of the exchange-correlation potential, e_xc
3171 : !> \param rho_set object containing the density at which the derivatives were calculated, \rho
3172 : !> \param rho1_set object containing the density with which to fold, \rho1s
3173 : !> \param pw_pool the pool for the grids
3174 : !> \param xc_section XC parameters
3175 : !> \par History
3176 : !> * 07.2024 Created [LHS]
3177 : ! **************************************************************************************************
3178 0 : SUBROUTINE xc_calc_3rd_deriv_analytical(v_xc, v_xc_tau, deriv_set, rho_set, rho1_set, &
3179 : pw_pool, xc_section, spinflip)
3180 :
3181 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_xc, v_xc_tau
3182 : TYPE(xc_derivative_set_type) :: deriv_set
3183 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set, rho1_set
3184 : TYPE(pw_pool_type), POINTER :: pw_pool
3185 : TYPE(section_vals_type), POINTER :: xc_section
3186 : LOGICAL, INTENT(in), OPTIONAL :: spinflip
3187 :
3188 : CHARACTER(len=*), PARAMETER :: routineN = 'xc_calc_3rd_deriv_analytical'
3189 :
3190 : INTEGER :: handle, i, idir, ispin, j, &
3191 : k, nspins, xc_deriv_method_id
3192 : INTEGER, DIMENSION(2, 3) :: bo
3193 : LOGICAL :: lsd, do_spinflip, alda0, &
3194 : rho_f, gradient_f, tau_f, laplace_f
3195 : REAL(KIND=dp) :: s, S_THRESH, S_THRESH2, gradient_cut
3196 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: dr1dr, dra1dra, drb1drb
3197 0 : REAL(kind=dp), DIMENSION(:, :, :), POINTER :: deriv_data, deriv_data2, e_drhoa, e_drhob, &
3198 0 : e_drho, norm_drho, norm_drhoa, &
3199 0 : norm_drhob, rho1a, rho1b, &
3200 0 : rhoa, rhob
3201 0 : TYPE(cp_3d_r_cp_type), DIMENSION(3) :: drho, drho1, drho1a, drho1b, drhoa, drhob
3202 0 : TYPE(pw_r3d_rs_type), DIMENSION(:), ALLOCATABLE :: v_drhoa, v_drhob, v_drho
3203 0 : TYPE(pw_r3d_rs_type), DIMENSION(:, :), ALLOCATABLE :: v_drho_r
3204 : TYPE(pw_c1d_gs_type) :: tmp_g, vxc_g
3205 : TYPE(xc_derivative_type), POINTER :: deriv_att
3206 :
3207 0 : CALL timeset(routineN, handle)
3208 :
3209 0 : NULLIFY (e_drhoa, e_drhob, e_drho)
3210 :
3211 0 : CPASSERT(ASSOCIATED(v_xc))
3212 0 : CPASSERT(ASSOCIATED(xc_section))
3213 :
3214 : ! Initialize parameters
3215 : CALL section_vals_val_get(xc_section, "XC_GRID%XC_DERIV", &
3216 0 : i_val=xc_deriv_method_id)
3217 : !
3218 0 : nspins = SIZE(v_xc)
3219 0 : lsd = ASSOCIATED(rho_set%rhoa)
3220 : !
3221 0 : do_spinflip = .FALSE.
3222 0 : IF (PRESENT(spinflip)) do_spinflip = spinflip
3223 : !
3224 0 : bo = rho_set%local_bounds
3225 : !
3226 0 : CALL check_for_derivatives(deriv_set, lsd, rho_f, gradient_f, tau_f, laplace_f)
3227 : !
3228 0 : CALL xc_rho_set_get(rho_set, drho_cutoff=gradient_cut)
3229 : !
3230 : !S_THRESH has to be the same as S_THRESH in xc_calc_2nd_deriv_analytical
3231 0 : alda0 = .false.
3232 0 : S_THRESH = 1.0E-04
3233 0 : S_THRESH2 = 1.0E-07
3234 :
3235 : ! Initialize potential
3236 0 : DO ispin = 1, nspins
3237 : !CALL pw_zero(v_xc(ispin))
3238 0 : v_xc(ispin)%array = 0.0_dp
3239 : END DO
3240 :
3241 : ! Create GGA fields
3242 0 : IF (gradient_f) THEN
3243 0 : ALLOCATE (v_drho_r(3, nspins), v_drho(nspins))
3244 0 : DO ispin = 1, nspins
3245 0 : DO idir = 1, 3
3246 0 : CALL allocate_pw(v_drho_r(idir, ispin), pw_pool, bo)
3247 : END DO
3248 0 : CALL allocate_pw(v_drho(ispin), pw_pool, bo)
3249 : END DO
3250 :
3251 0 : IF (xc_requires_tmp_g(xc_deriv_method_id)) THEN
3252 0 : IF (ASSOCIATED(pw_pool)) THEN
3253 0 : CALL pw_pool%create_pw(tmp_g)
3254 0 : CALL pw_pool%create_pw(vxc_g)
3255 : ELSE
3256 : ! remember to refix for gapw
3257 0 : CPABORT("XC_DERIV method is not implemented in GAPW")
3258 : END IF
3259 : END IF
3260 :
3261 : END IF
3262 :
3263 : ! Initialize mGGA potential
3264 0 : IF (tau_f) THEN
3265 0 : CPASSERT(ASSOCIATED(v_xc_tau))
3266 0 : DO ispin = 1, nspins
3267 0 : v_xc_tau(ispin)%array = 0.0_dp
3268 : END DO
3269 : END IF
3270 :
3271 0 : IF (lsd) THEN
3272 :
3273 : !-------------------!
3274 : ! UNrestricted case !
3275 : !-------------------!
3276 :
3277 0 : IF (do_spinflip) THEN
3278 0 : CALL xc_rho_set_get(rho1_set, rhoa=rho1a)
3279 0 : CALL xc_rho_set_get(rho_set, rhoa=rhoa, rhob=rhob)
3280 : ELSE
3281 0 : CALL xc_rho_set_get(rho1_set, rhoa=rho1a, rhob=rho1b)
3282 : END IF
3283 :
3284 0 : IF (gradient_f) THEN
3285 : CALL xc_rho_set_get(rho_set, drhoa=drhoa, drhob=drhob, &
3286 0 : norm_drho=norm_drho, norm_drhoa=norm_drhoa, norm_drhob=norm_drhob)
3287 0 : IF (do_spinflip) THEN
3288 0 : CALL xc_rho_set_get(rho1_set, drhoa=drho1a)
3289 0 : CALL calc_drho_from_a(drho1, drho1a)
3290 : ELSE
3291 0 : CALL xc_rho_set_get(rho1_set, drhoa=drho1a, drhob=drho1b)
3292 0 : CALL calc_drho_from_ab(drho1, drho1a, drho1b)
3293 : END IF
3294 :
3295 0 : CALL calc_drho_from_ab(drho, drhoa, drhob)
3296 :
3297 0 : CALL prepare_dr1dr(dra1dra, drhoa, drho1a)
3298 0 : IF (do_spinflip) THEN
3299 0 : CALL prepare_dr1dr(drb1drb, drhob, drho1a)
3300 0 : CALL prepare_dr1dr(dr1dr, drho, drho1a)
3301 0 : ELSE IF (nspins /= 1) THEN
3302 0 : CALL prepare_dr1dr(drb1drb, drhob, drho1b)
3303 0 : CALL prepare_dr1dr(dr1dr, drho, drho1)
3304 : ELSE
3305 0 : CPABORT("Exchange-correlation's third derivative for closed-shell not yet implemented")
3306 : END IF
3307 :
3308 : ! Create vectors for partial integration term
3309 0 : ALLOCATE (v_drhoa(nspins), v_drhob(nspins))
3310 0 : DO ispin = 1, nspins
3311 0 : CALL allocate_pw(v_drhoa(ispin), pw_pool, bo)
3312 0 : CALL allocate_pw(v_drhob(ispin), pw_pool, bo)
3313 : END DO
3314 :
3315 : END IF
3316 :
3317 0 : IF (laplace_f) THEN
3318 0 : CPABORT("Exchange-correlation's laplace analytic third derivative not implemented")
3319 : END IF
3320 :
3321 0 : IF (tau_f) THEN
3322 0 : CPABORT("Exchange-correlation's mGGA analytic third derivative not implemented")
3323 : END IF
3324 :
3325 0 : IF (nspins /= 1) THEN
3326 :
3327 0 : IF (.NOT. do_spinflip) THEN
3328 : ! Analytic third derivative of the excchange-correlation functional
3329 0 : CPABORT("Exchange-correlation's analytic third derivative not implemented")
3330 :
3331 : ELSE
3332 :
3333 : ! vxc contributions
3334 : ! vxca = (vxc^{\alpha}-vxc^{\beta})*rho1/|rhoa-rhob|^2
3335 : ! vxcb =-(vxc^{\alpha}-vxc^{\beta})*rho1/|rhoa-rhob|^2
3336 : ! Alpha LDA contribution
3337 : ! | d e_xc d e_xc | rho1a
3338 : ! vxca = rho1a*|-------- - --------|*---------------
3339 : ! | drhoa drhob | |rhoa - rhob|^2
3340 0 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_rhoa])
3341 0 : IF (ASSOCIATED(deriv_att)) THEN
3342 0 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
3343 0 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_rhob])
3344 0 : IF (ASSOCIATED(deriv_att)) THEN
3345 0 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data2)
3346 : !$OMP PARALLEL DO PRIVATE(k,j,i,s) DEFAULT(NONE)&
3347 0 : !$OMP SHARED(bo,v_xc,deriv_data,deriv_data2,rho1a,rhoa,rhob,S_THRESH2) COLLAPSE(3)
3348 : DO k = bo(1, 3), bo(2, 3)
3349 : DO j = bo(1, 2), bo(2, 2)
3350 : DO i = bo(1, 1), bo(2, 1)
3351 : s = rhoa(i, j, k) - rhob(i, j, k)
3352 : s = -SIGN(MAX(s**2, S_THRESH2), s)
3353 : v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
3354 : (deriv_data(i, j, k) - deriv_data2(i, j, k))*rho1a(i, j, k)**2/s
3355 : END DO
3356 : END DO
3357 : END DO
3358 : !$OMP END PARALLEL DO
3359 : END IF
3360 : END IF
3361 : ! GGA contributions to the spin-flip xcKernel
3362 : ! Alpha GGA contributions
3363 : ! | d e_xc d e_xc | rho1a
3364 : ! vxca += + |----------*dra1dra - ----------*drb1drb|*---------------
3365 : ! | d|drhoa| d|drhob| | |rhoa - rhob|^2
3366 : IF (.NOT. alda0) THEN
3367 0 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa])
3368 0 : IF (ASSOCIATED(deriv_att)) THEN
3369 0 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
3370 0 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob])
3371 0 : IF (ASSOCIATED(deriv_att)) THEN
3372 0 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data2)
3373 : !$OMP PARALLEL DO PRIVATE(k,j,i,s) DEFAULT(NONE)&
3374 0 : !$OMP SHARED(bo,deriv_data,deriv_data2,dra1dra,drb1drb,v_xc,rho1a,rhoa,rhob,S_THRESH2) COLLAPSE(3)
3375 : DO k = bo(1, 3), bo(2, 3)
3376 : DO j = bo(1, 2), bo(2, 2)
3377 : DO i = bo(1, 1), bo(2, 1)
3378 : s = rhoa(i, j, k) - rhob(i, j, k)
3379 : s = -SIGN(MAX(s**2, S_THRESH2), s)
3380 : v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
3381 : (deriv_data(i, j, k)*dra1dra(i, j, k) - deriv_data2(i, j, k)*drb1drb(i, j, k)) &
3382 : *rho1a(i, j, k)/s
3383 : END DO
3384 : END DO
3385 : END DO
3386 : !$OMP END PARALLEL DO
3387 : END IF
3388 : END IF
3389 : END IF
3390 : ! Beta contribution = - alpha
3391 : !$OMP PARALLEL DO PRIVATE(k,j,i,s) DEFAULT(NONE)&
3392 0 : !$OMP SHARED(bo,v_xc) COLLAPSE(3)
3393 : DO k = bo(1, 3), bo(2, 3)
3394 : DO j = bo(1, 2), bo(2, 2)
3395 : DO i = bo(1, 1), bo(2, 1)
3396 : v_xc(2)%array(i, j, k) = -v_xc(1)%array(i, j, k)
3397 : END DO
3398 : END DO
3399 : END DO
3400 : !$OMP END PARALLEL DO
3401 : ! fxc contributions
3402 : ! vxca = rho1*(fxc^{\alpha\alpha}-fxc^{\alpha\beta})*rho1/|rhoa-rhob|
3403 : ! vxcb = rho1*(fxc^{\beta\alpha}-fxc^{\beta\beta})*rho1/|rhoa-rhob|
3404 : ! Alpha LDA contribution
3405 : ! | d^2 e_xc d^2 e_xc | rho1a
3406 : ! vxca += rho1a*|------------- - -------------|*-------------
3407 : ! | drhoa drhoa drhoa drhob | |rhoa - rhob|
3408 0 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhob])
3409 0 : IF (ASSOCIATED(deriv_att)) THEN
3410 0 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data2)
3411 0 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa])
3412 0 : IF (ASSOCIATED(deriv_att)) THEN
3413 0 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
3414 : !$OMP PARALLEL DO PRIVATE(k,j,i,s) DEFAULT(NONE)&
3415 0 : !$OMP SHARED(bo,deriv_data,deriv_data2,rho1a,v_xc,rhoa,rhob,S_THRESH) COLLAPSE(3)
3416 : DO k = bo(1, 3), bo(2, 3)
3417 : DO j = bo(1, 2), bo(2, 2)
3418 : DO i = bo(1, 1), bo(2, 1)
3419 : s = MAX(ABS(rhoa(i, j, k) - rhob(i, j, k)), S_THRESH)
3420 : v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
3421 : rho1a(i, j, k)**2*(deriv_data(i, j, k) - deriv_data2(i, j, k))/s
3422 : END DO
3423 : END DO
3424 : END DO
3425 : !$OMP END PARALLEL DO
3426 : END IF
3427 : END IF
3428 : ! Beta LDA contribution
3429 : ! | d^2 e_xc d^2 e_xc | rho1a
3430 : ! vxcb += rho1a*|------------- - -------------|*-------------
3431 : ! | drhob drhoa drhob drhob | |rhoa - rhob|
3432 0 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob])
3433 0 : IF (ASSOCIATED(deriv_att)) THEN
3434 0 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data2)
3435 0 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhoa])
3436 0 : IF (ASSOCIATED(deriv_att)) THEN
3437 0 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
3438 : !$OMP PARALLEL DO PRIVATE(k,j,i,s) DEFAULT(NONE)&
3439 0 : !$OMP SHARED(bo,deriv_data,deriv_data2,rho1a,v_xc,rhoa,rhob,S_THRESH) COLLAPSE(3)
3440 : DO k = bo(1, 3), bo(2, 3)
3441 : DO j = bo(1, 2), bo(2, 2)
3442 : DO i = bo(1, 1), bo(2, 1)
3443 : s = MAX(ABS(rhoa(i, j, k) - rhob(i, j, k)), S_THRESH)
3444 : v_xc(2)%array(i, j, k) = v_xc(2)%array(i, j, k) + &
3445 : rho1a(i, j, k)**2*(deriv_data(i, j, k) - deriv_data2(i, j, k))/s
3446 : END DO
3447 : END DO
3448 : END DO
3449 : !$OMP END PARALLEL DO
3450 : END IF
3451 : END IF
3452 : ! Alpha GGA contribution
3453 : IF (.NOT. alda0) THEN
3454 : ! rho1a | d^2 e_xc d^2 e_xc |
3455 : ! vxca += + -------------*|----------------*dra1dra - ----------------*drb1drb|
3456 : ! |rhoa - rhob| | drhoa d|drhoa| drhoa d|drhob| |
3457 0 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_norm_drhob])
3458 0 : IF (ASSOCIATED(deriv_att)) THEN
3459 0 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data2)
3460 0 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_norm_drhoa])
3461 0 : IF (ASSOCIATED(deriv_att)) THEN
3462 0 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
3463 : !$OMP PARALLEL DO PRIVATE(k,j,i,s) DEFAULT(NONE)&
3464 0 : !$OMP SHARED(bo,deriv_data,deriv_data2,dra1dra,drb1drb,rho1a,v_xc,rhoa,rhob,S_THRESH) COLLAPSE(3)
3465 : DO k = bo(1, 3), bo(2, 3)
3466 : DO j = bo(1, 2), bo(2, 2)
3467 : DO i = bo(1, 1), bo(2, 1)
3468 : s = MAX(ABS(rhoa(i, j, k) - rhob(i, j, k)), S_THRESH)
3469 : v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
3470 : (deriv_data(i, j, k)*dra1dra(i, j, k) - deriv_data2(i, j, k)*drb1drb(i, j, k))* &
3471 : rho1a(i, j, k)/s
3472 : END DO
3473 : END DO
3474 : END DO
3475 : !$OMP END PARALLEL DO
3476 : END IF
3477 : END IF
3478 : ! Beta GGA contribution
3479 : ! rho1a | d^2 e_xc d^2 e_xc |
3480 : ! vxcb += + -------------*|----------------*dra1dra - ----------------*drb1drb|
3481 : ! |rhoa - rhob| | drhob d|drhoa| drhob d|drhob| |
3482 0 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_norm_drhob])
3483 0 : IF (ASSOCIATED(deriv_att)) THEN
3484 0 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data2)
3485 0 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_norm_drhoa])
3486 0 : IF (ASSOCIATED(deriv_att)) THEN
3487 0 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
3488 : !$OMP PARALLEL DO PRIVATE(k,j,i,s) DEFAULT(NONE)&
3489 0 : !$OMP SHARED(bo,deriv_data,deriv_data2,dra1dra,drb1drb,rho1a,v_xc,rhoa,rhob,S_THRESH) COLLAPSE(3)
3490 : DO k = bo(1, 3), bo(2, 3)
3491 : DO j = bo(1, 2), bo(2, 2)
3492 : DO i = bo(1, 1), bo(2, 1)
3493 : s = MAX(ABS(rhoa(i, j, k) - rhob(i, j, k)), S_THRESH)
3494 : v_xc(2)%array(i, j, k) = v_xc(2)%array(i, j, k) + &
3495 : (deriv_data(i, j, k)*dra1dra(i, j, k) - deriv_data2(i, j, k)*drb1drb(i, j, k))* &
3496 : rho1a(i, j, k)/s
3497 : END DO
3498 : END DO
3499 : END DO
3500 : !$OMP END PARALLEL DO
3501 : END IF
3502 : END IF
3503 : !
3504 : !
3505 : ! Calculate the vector for the partial integration term of GGA functionals
3506 : ! First contribution alpha
3507 : ! | d^2 e_xc d^2 e_xc |
3508 : ! v_drhoa(1) += -|---------------- - ----------------|*rho1a
3509 : ! | d|drhoa| drhoa d|drhoa| drhob |
3510 0 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_rhob])
3511 0 : IF (ASSOCIATED(deriv_att)) THEN
3512 0 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data2)
3513 0 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_rhoa])
3514 0 : IF (ASSOCIATED(deriv_att)) THEN
3515 0 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
3516 : !$OMP PARALLEL DO PRIVATE(k,j,i) DEFAULT(NONE)&
3517 0 : !$OMP SHARED(bo,deriv_data,deriv_data2,rho1a,v_drhoa) COLLAPSE(3)
3518 : DO k = bo(1, 3), bo(2, 3)
3519 : DO j = bo(1, 2), bo(2, 2)
3520 : DO i = bo(1, 1), bo(2, 1)
3521 : v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
3522 : (deriv_data(i, j, k) - deriv_data2(i, j, k))*rho1a(i, j, k)
3523 : END DO
3524 : END DO
3525 : END DO
3526 : !$OMP END PARALLEL DO
3527 : END IF
3528 : END IF
3529 : ! First contribution beta
3530 : ! | d^2 e_xc d^2 e_xc |
3531 : ! v_drhob(2) += +|---------------- - ----------------|*rho1a
3532 : ! | d|drhob| drhob d|drhob| drhoa |
3533 0 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob, deriv_rhoa])
3534 0 : IF (ASSOCIATED(deriv_att)) THEN
3535 0 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data2)
3536 0 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob, deriv_rhob])
3537 0 : IF (ASSOCIATED(deriv_att)) THEN
3538 0 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
3539 : !$OMP PARALLEL DO PRIVATE(k,j,i) DEFAULT(NONE)&
3540 0 : !$OMP SHARED(bo,deriv_data,deriv_data2,rho1a,v_drhob) COLLAPSE(3)
3541 : DO k = bo(1, 3), bo(2, 3)
3542 : DO j = bo(1, 2), bo(2, 2)
3543 : DO i = bo(1, 1), bo(2, 1)
3544 : v_drhob(2)%array(i, j, k) = v_drhob(2)%array(i, j, k) + &
3545 : (deriv_data(i, j, k) - deriv_data2(i, j, k))*rho1a(i, j, k)
3546 : END DO
3547 : END DO
3548 : END DO
3549 : !$OMP END PARALLEL DO
3550 : END IF
3551 : END IF
3552 : ! First contribution spinless
3553 : ! | d^2 e_xc d^2 e_xc |
3554 : ! v_drho(1) += -|--------------- - ---------------|*rho1a
3555 : ! | d|drho| drhoa d|drho| drhob |
3556 : !
3557 : ! | d^2 e_xc d^2 e_xc |
3558 : ! v_drho(2) += -|--------------- - ---------------|*rho1a
3559 : ! | d|drho| drhoa d|drho| drhob |
3560 0 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_rhoa])
3561 0 : IF (ASSOCIATED(deriv_att)) THEN
3562 0 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
3563 0 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_rhob])
3564 0 : IF (ASSOCIATED(deriv_att)) THEN
3565 0 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data2)
3566 : !$OMP PARALLEL DO PRIVATE(k,j,i) DEFAULT(NONE)&
3567 0 : !$OMP SHARED(bo,deriv_data,deriv_data2,rho1a,v_drho) COLLAPSE(3)
3568 : DO k = bo(1, 3), bo(2, 3)
3569 : DO j = bo(1, 2), bo(2, 2)
3570 : DO i = bo(1, 1), bo(2, 1)
3571 : v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
3572 : (deriv_data(i, j, k) - deriv_data2(i, j, k))*rho1a(i, j, k)
3573 : v_drho(2)%array(i, j, k) = v_drho(2)%array(i, j, k) - &
3574 : (deriv_data(i, j, k) - deriv_data2(i, j, k))*rho1a(i, j, k)
3575 : END DO
3576 : END DO
3577 : END DO
3578 : !$OMP END PARALLEL DO
3579 : END IF
3580 : END IF
3581 : ! Second contribution
3582 0 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_norm_drhob])
3583 0 : IF (ASSOCIATED(deriv_att)) THEN
3584 0 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data2)
3585 : ! d^2 e_xc d^2 e_xc
3586 : ! v_drhoa(1) += - -------------------*dra1dra + ------------------*drb1drb
3587 : ! d|drhoa| d|drhoa| d|drhoa| d|drhob|
3588 0 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_norm_drhoa])
3589 0 : IF (ASSOCIATED(deriv_att)) THEN
3590 0 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
3591 : !$OMP PARALLEL DO PRIVATE(k,j,i) DEFAULT(NONE)&
3592 0 : !$OMP SHARED(bo,deriv_data,deriv_data2,dra1dra,drb1drb,v_drhoa) COLLAPSE(3)
3593 : DO k = bo(1, 3), bo(2, 3)
3594 : DO j = bo(1, 2), bo(2, 2)
3595 : DO i = bo(1, 1), bo(2, 1)
3596 : v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
3597 : deriv_data(i, j, k)*dra1dra(i, j, k) + deriv_data2(i, j, k)*drb1drb(i, j, k)
3598 : END DO
3599 : END DO
3600 : END DO
3601 : !$OMP END PARALLEL DO
3602 : END IF
3603 : ! d^2 e_xc d^2 e_xc
3604 : ! v_drhob(2) += - -------------------*dra1dra + -------------------*drb1drb
3605 : ! d|drhoa| d|drhob| d|drhob| d|drhob|
3606 0 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob, deriv_norm_drhob])
3607 0 : IF (ASSOCIATED(deriv_att)) THEN
3608 0 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
3609 : !$OMP PARALLEL DO PRIVATE(k,j,i) DEFAULT(NONE)&
3610 0 : !$OMP SHARED(bo,deriv_data,deriv_data2,dra1dra,drb1drb,v_drhob) COLLAPSE(3)
3611 : DO k = bo(1, 3), bo(2, 3)
3612 : DO j = bo(1, 2), bo(2, 2)
3613 : DO i = bo(1, 1), bo(2, 1)
3614 : v_drhob(2)%array(i, j, k) = v_drhob(2)%array(i, j, k) - &
3615 : deriv_data(i, j, k)*dra1dra(i, j, k) + deriv_data2(i, j, k)*drb1drb(i, j, k)
3616 : END DO
3617 : END DO
3618 : END DO
3619 : !$OMP END PARALLEL DO
3620 : END IF
3621 : END IF
3622 : ! d^2 e_xc d^2 e_xc
3623 : ! v_drho(1) += - ------------------*dra1dra + ------------------*drb1drb
3624 : ! d|drho| d|drhoa| d|drho| d|drhob|
3625 : !
3626 : ! d^2 e_xc d^2 e_xc
3627 : ! v_drho(2) += - ------------------*dra1dra + ------------------*drb1drb
3628 : ! d|drho| d|drhoa| d|drho| d|drhob|
3629 0 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_norm_drhob])
3630 0 : IF (ASSOCIATED(deriv_att)) THEN
3631 0 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data2)
3632 0 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_norm_drhoa])
3633 0 : IF (ASSOCIATED(deriv_att)) THEN
3634 0 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
3635 : !$OMP PARALLEL DO PRIVATE(k,j,i) DEFAULT(NONE)&
3636 0 : !$OMP SHARED(bo,deriv_data,deriv_data2,dra1dra,drb1drb,v_drho) COLLAPSE(3)
3637 : DO k = bo(1, 3), bo(2, 3)
3638 : DO j = bo(1, 2), bo(2, 2)
3639 : DO i = bo(1, 1), bo(2, 1)
3640 : v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
3641 : deriv_data(i, j, k)*dra1dra(i, j, k) + &
3642 : deriv_data2(i, j, k)*drb1drb(i, j, k)
3643 : v_drho(2)%array(i, j, k) = v_drho(2)%array(i, j, k) - &
3644 : deriv_data(i, j, k)*dra1dra(i, j, k) + &
3645 : deriv_data2(i, j, k)*drb1drb(i, j, k)
3646 : END DO
3647 : END DO
3648 : END DO
3649 : !$OMP END PARALLEL DO
3650 : END IF
3651 : END IF
3652 : !
3653 :
3654 : ! Last GGA contribution
3655 : ! Alpha contribution
3656 : ! d e_xc
3657 : ! v_drhoa(1) += + ----------*dra1dra
3658 : ! d|drhoa|
3659 0 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa])
3660 0 : IF (ASSOCIATED(deriv_att)) THEN
3661 0 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
3662 0 : CALL xc_derivative_get(deriv_att, deriv_data=e_drhoa)
3663 :
3664 0 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(dra1dra,gradient_cut,norm_drhoa,v_drhoa,deriv_data)
3665 : v_drhoa(1)%array(:, :, :) = v_drhoa(1)%array(:, :, :) + &
3666 : deriv_data(:, :, :)*dra1dra(:, :, :)/MAX(gradient_cut, norm_drhoa(:, :, :))**2
3667 : !$OMP END PARALLEL WORKSHARE
3668 : END IF
3669 : ! Beta contribution
3670 : ! d e_xc
3671 : ! v_drhob(2) += - ----------*drb1drb
3672 : ! d|drhob|
3673 0 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob])
3674 0 : IF (ASSOCIATED(deriv_att)) THEN
3675 0 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
3676 0 : CALL xc_derivative_get(deriv_att, deriv_data=e_drhob)
3677 :
3678 0 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(drb1drb,gradient_cut,norm_drhob,v_drhob,deriv_data)
3679 : v_drhob(2)%array(:, :, :) = v_drhob(2)%array(:, :, :) - &
3680 : deriv_data(:, :, :)*drb1drb(:, :, :)/MAX(gradient_cut, norm_drhob(:, :, :))**2
3681 : !$OMP END PARALLEL WORKSHARE
3682 : END IF
3683 : END IF ! If ALDA0
3684 : END IF
3685 :
3686 : ELSE
3687 :
3688 : ! Analytic third derivative for closed-shell
3689 0 : CPABORT("Exchange-correlation's analytic third derivative not implemented")
3690 :
3691 : END IF
3692 :
3693 0 : IF (gradient_f) THEN
3694 : IF (.NOT. alda0) THEN
3695 :
3696 : ! partial integration
3697 0 : DO idir = 1, 3
3698 :
3699 : ! GGA contributions to the spin-flip xc-Kernel
3700 : !
3701 : ! v_drhoa(1)*drhoa(:)*rhoa1 v_drhob(1)*drhob(:)*rhoa1 v_drho(1)*drho(:)*rhoa1
3702 : ! v_drho_r(:,1) = --------------------------- + --------------------------- + -------------------------
3703 : ! |rhoa - rhob| |rhoa - rhob| |rhoa - rhob|
3704 : !
3705 : ! v_drhoa(2)*drhoa(:)*rhoa1 v_drhob(2)*drhob(:)*rhoa1 v_drho(2)*drho(:)*rhoa1
3706 : ! v_drho_r(:,2) = --------------------------- + --------------------------- + -------------------------
3707 : ! |rhoa - rhob| |rhoa - rhob| |rhoa - rhob|
3708 0 : IF (do_spinflip) THEN
3709 : !$OMP PARALLEL DO PRIVATE(k,j,i,s) DEFAULT(NONE)&
3710 0 : !$OMP SHARED(bo,v_drho_r,v_drho,v_drhoa,v_drhob,rhoa,rhob,drho,drhoa,drhob,rho1a,idir,S_THRESH) COLLAPSE(3)
3711 : DO k = bo(1, 3), bo(2, 3)
3712 : DO j = bo(1, 2), bo(2, 2)
3713 : DO i = bo(1, 1), bo(2, 1)
3714 : s = MAX(ABS(rhoa(i, j, k) - rhob(i, j, k)), S_THRESH)
3715 : DO ispin = 1, 2
3716 : v_drho_r(idir, ispin)%array(i, j, k) = v_drho_r(idir, ispin)%array(i, j, k) + &
3717 : v_drhoa(ispin)%array(i, j, k)*drhoa(idir)%array(i, j, k)*rho1a(i, j, k)/s + &
3718 : v_drhob(ispin)%array(i, j, k)*drhob(idir)%array(i, j, k)*rho1a(i, j, k)/s + &
3719 : v_drho(ispin)%array(i, j, k)*drho(idir)%array(i, j, k)*rho1a(i, j, k)/s
3720 : END DO
3721 : END DO
3722 : END DO
3723 : END DO
3724 : !$OMP END PARALLEL DO
3725 : END IF
3726 : ! Last GGA contribution
3727 : ! Alpha contribution
3728 : ! rho1a d e_xc
3729 : ! v_drho_r(:,1) += - -------------*----------*drho1a(:)
3730 : ! |rhoa - rhob| d|drhoa|
3731 0 : IF (ASSOCIATED(e_drhoa)) THEN
3732 : !$OMP PARALLEL DO PRIVATE(k,j,i,s) DEFAULT(NONE)&
3733 0 : !$OMP SHARED(bo,e_drhoa,v_drho_r,drho1a,rho1a,rhoa,rhob,S_THRESH,idir) COLLAPSE(3)
3734 : DO k = bo(1, 3), bo(2, 3)
3735 : DO j = bo(1, 2), bo(2, 2)
3736 : DO i = bo(1, 1), bo(2, 1)
3737 : s = MAX(ABS(rhoa(i, j, k) - rhob(i, j, k)), S_THRESH)
3738 : v_drho_r(idir, 1)%array(i, j, k) = v_drho_r(idir, 1)%array(i, j, k) - &
3739 : e_drhoa(i, j, k)*drho1a(idir)%array(i, j, k)*rho1a(i, j, k)/s
3740 : END DO
3741 : END DO
3742 : END DO
3743 : !$OMP END PARALLEL DO
3744 : END IF
3745 : ! Beta contribution
3746 : ! rho1a d e_xc
3747 : ! v_drho_r(:,2) += + -------------*----------*drho1a(:)
3748 : ! |rhoa - rhob| d|drhob|
3749 0 : IF (ASSOCIATED(e_drhob)) THEN
3750 : !$OMP PARALLEL DO PRIVATE(k,j,i,s) DEFAULT(NONE)&
3751 0 : !$OMP SHARED(bo,e_drhob,v_drho_r,drho1a,rho1a,rhoa,rhob,S_THRESH,idir) COLLAPSE(3)
3752 : DO k = bo(1, 3), bo(2, 3)
3753 : DO j = bo(1, 2), bo(2, 2)
3754 : DO i = bo(1, 1), bo(2, 1)
3755 : s = MAX(ABS(rhoa(i, j, k) - rhob(i, j, k)), S_THRESH)
3756 : v_drho_r(idir, 2)%array(i, j, k) = v_drho_r(idir, 2)%array(i, j, k) + &
3757 : e_drhob(i, j, k)*drho1a(idir)%array(i, j, k)*rho1a(i, j, k)/s
3758 : END DO
3759 : END DO
3760 : END DO
3761 : !$OMP END PARALLEL DO
3762 : END IF
3763 : END DO
3764 :
3765 : ! partial integration: v_xc = v_xc - \nabla \cdot vdrho_r
3766 0 : DO ispin = 1, nspins
3767 0 : CALL xc_pw_divergence(xc_deriv_method_id, v_drho_r(:, ispin), tmp_g, vxc_g, v_xc(ispin))
3768 : END DO ! ispin
3769 : END IF ! ALDA0
3770 :
3771 0 : DO idir = 1, 3
3772 0 : DEALLOCATE (drho(idir)%array)
3773 0 : DEALLOCATE (drho1(idir)%array)
3774 : END DO
3775 :
3776 0 : DO ispin = 1, nspins
3777 0 : CALL deallocate_pw(v_drhoa(ispin), pw_pool)
3778 0 : CALL deallocate_pw(v_drhob(ispin), pw_pool)
3779 : END DO
3780 :
3781 0 : DEALLOCATE (v_drhoa, v_drhob)
3782 :
3783 : END IF ! gradient_f
3784 :
3785 : ELSE
3786 :
3787 : !-----------------!
3788 : ! restricted case !
3789 : !-----------------!
3790 0 : CPABORT("Exchange-correlation's analytic third derivative not implemented")
3791 :
3792 : END IF
3793 :
3794 0 : IF (gradient_f) THEN
3795 :
3796 0 : DO ispin = 1, nspins
3797 0 : CALL deallocate_pw(v_drho(ispin), pw_pool)
3798 0 : DO idir = 1, 3
3799 0 : CALL deallocate_pw(v_drho_r(idir, ispin), pw_pool)
3800 : END DO
3801 : END DO
3802 0 : DEALLOCATE (v_drho, v_drho_r)
3803 :
3804 : END IF
3805 :
3806 0 : IF (ASSOCIATED(tmp_g%pw_grid) .AND. ASSOCIATED(pw_pool)) THEN
3807 0 : CALL pw_pool%give_back_pw(tmp_g)
3808 : END IF
3809 :
3810 0 : IF (ASSOCIATED(vxc_g%pw_grid) .AND. ASSOCIATED(pw_pool)) THEN
3811 0 : CALL pw_pool%give_back_pw(vxc_g)
3812 : END IF
3813 :
3814 0 : CALL timestop(handle)
3815 0 : END SUBROUTINE xc_calc_3rd_deriv_analytical
3816 :
3817 : ! **************************************************************************************************
3818 : !> \brief allocates grids using pw_pool (if associated) or with bounds
3819 : !> \param pw ...
3820 : !> \param pw_pool ...
3821 : !> \param bo ...
3822 : ! **************************************************************************************************
3823 101546 : SUBROUTINE allocate_pw(pw, pw_pool, bo)
3824 : TYPE(pw_r3d_rs_type), INTENT(OUT) :: pw
3825 : TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool
3826 : INTEGER, DIMENSION(2, 3), INTENT(IN) :: bo
3827 :
3828 101546 : IF (ASSOCIATED(pw_pool)) THEN
3829 65082 : CALL pw_pool%create_pw(pw)
3830 65082 : CALL pw_zero(pw)
3831 : ELSE
3832 182320 : ALLOCATE (pw%array(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
3833 93056128 : pw%array = 0.0_dp
3834 : END IF
3835 :
3836 101546 : END SUBROUTINE allocate_pw
3837 :
3838 : ! **************************************************************************************************
3839 : !> \brief deallocates grid allocated with allocate_pw
3840 : !> \param pw ...
3841 : !> \param pw_pool ...
3842 : ! **************************************************************************************************
3843 101546 : SUBROUTINE deallocate_pw(pw, pw_pool)
3844 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: pw
3845 : TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool
3846 :
3847 101546 : IF (ASSOCIATED(pw_pool)) THEN
3848 65082 : CALL pw_pool%give_back_pw(pw)
3849 : ELSE
3850 36464 : CALL pw%release()
3851 : END IF
3852 :
3853 101546 : END SUBROUTINE deallocate_pw
3854 :
3855 : ! **************************************************************************************************
3856 : !> \brief updates virial from first derivative w.r.t. norm_drho
3857 : !> \param virial_pw ...
3858 : !> \param drho ...
3859 : !> \param drho1 ...
3860 : !> \param deriv_data ...
3861 : !> \param virial_xc ...
3862 : ! **************************************************************************************************
3863 304 : SUBROUTINE virial_drho_drho1(virial_pw, drho, drho1, deriv_data, virial_xc)
3864 : TYPE(pw_r3d_rs_type), INTENT(IN) :: virial_pw
3865 : TYPE(cp_3d_r_cp_type), DIMENSION(3), INTENT(IN) :: drho, drho1
3866 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: deriv_data
3867 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT) :: virial_xc
3868 :
3869 : INTEGER :: idir, jdir
3870 : REAL(KIND=dp) :: tmp
3871 :
3872 1216 : DO idir = 1, 3
3873 912 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(drho,idir,virial_pw,deriv_data)
3874 : virial_pw%array(:, :, :) = drho(idir)%array(:, :, :)*deriv_data(:, :, :)
3875 : !$OMP END PARALLEL WORKSHARE
3876 3952 : DO jdir = 1, 3
3877 : tmp = virial_pw%pw_grid%dvol*accurate_dot_product(virial_pw%array(:, :, :), &
3878 2736 : drho1(jdir)%array(:, :, :))
3879 2736 : virial_xc(jdir, idir) = virial_xc(jdir, idir) + tmp
3880 3648 : virial_xc(idir, jdir) = virial_xc(idir, jdir) + tmp
3881 : END DO
3882 : END DO
3883 :
3884 304 : END SUBROUTINE virial_drho_drho1
3885 :
3886 : ! **************************************************************************************************
3887 : !> \brief Adds virial contribution from second order potential parts
3888 : !> \param virial_pw ...
3889 : !> \param drho ...
3890 : !> \param v_drho ...
3891 : !> \param virial_xc ...
3892 : ! **************************************************************************************************
3893 298 : SUBROUTINE virial_drho_drho(virial_pw, drho, v_drho, virial_xc)
3894 : TYPE(pw_r3d_rs_type), INTENT(IN) :: virial_pw
3895 : TYPE(cp_3d_r_cp_type), DIMENSION(3), INTENT(IN) :: drho
3896 : TYPE(pw_r3d_rs_type), INTENT(IN) :: v_drho
3897 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT) :: virial_xc
3898 :
3899 : INTEGER :: idir, jdir
3900 : REAL(KIND=dp) :: tmp
3901 :
3902 1192 : DO idir = 1, 3
3903 894 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(drho,idir,v_drho,virial_pw)
3904 : virial_pw%array(:, :, :) = drho(idir)%array(:, :, :)*v_drho%array(:, :, :)
3905 : !$OMP END PARALLEL WORKSHARE
3906 2980 : DO jdir = 1, idir
3907 : tmp = -virial_pw%pw_grid%dvol*accurate_dot_product(virial_pw%array(:, :, :), &
3908 1788 : drho(jdir)%array(:, :, :))
3909 1788 : virial_xc(jdir, idir) = virial_xc(jdir, idir) + tmp
3910 2682 : virial_xc(idir, jdir) = virial_xc(jdir, idir)
3911 : END DO
3912 : END DO
3913 :
3914 298 : END SUBROUTINE virial_drho_drho
3915 :
3916 : ! **************************************************************************************************
3917 : !> \brief ...
3918 : !> \param rho_r ...
3919 : !> \param pw_pool ...
3920 : !> \param virial_xc ...
3921 : !> \param deriv_data ...
3922 : ! **************************************************************************************************
3923 150 : SUBROUTINE virial_laplace(rho_r, pw_pool, virial_xc, deriv_data)
3924 : TYPE(pw_r3d_rs_type), TARGET :: rho_r
3925 : TYPE(pw_pool_type), POINTER, INTENT(IN) :: pw_pool
3926 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT) :: virial_xc
3927 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: deriv_data
3928 :
3929 : CHARACTER(len=*), PARAMETER :: routineN = 'virial_laplace'
3930 :
3931 : INTEGER :: handle, idir, jdir
3932 : TYPE(pw_r3d_rs_type), POINTER :: virial_pw
3933 : TYPE(pw_c1d_gs_type), POINTER :: tmp_g, rho_g
3934 : INTEGER, DIMENSION(3) :: my_deriv
3935 :
3936 150 : CALL timeset(routineN, handle)
3937 :
3938 150 : NULLIFY (virial_pw, tmp_g, rho_g)
3939 150 : ALLOCATE (virial_pw, tmp_g, rho_g)
3940 150 : CALL pw_pool%create_pw(virial_pw)
3941 150 : CALL pw_pool%create_pw(tmp_g)
3942 150 : CALL pw_pool%create_pw(rho_g)
3943 150 : CALL pw_zero(virial_pw)
3944 150 : CALL pw_transfer(rho_r, rho_g)
3945 600 : DO idir = 1, 3
3946 1500 : DO jdir = idir, 3
3947 900 : CALL pw_copy(rho_g, tmp_g)
3948 :
3949 900 : my_deriv = 0
3950 900 : my_deriv(idir) = 1
3951 900 : my_deriv(jdir) = my_deriv(jdir) + 1
3952 :
3953 900 : CALL pw_derive(tmp_g, my_deriv)
3954 900 : CALL pw_transfer(tmp_g, virial_pw)
3955 : virial_xc(idir, jdir) = virial_xc(idir, jdir) - 2.0_dp*virial_pw%pw_grid%dvol* &
3956 : accurate_dot_product(virial_pw%array(:, :, :), &
3957 900 : deriv_data(:, :, :))
3958 1350 : virial_xc(jdir, idir) = virial_xc(idir, jdir)
3959 : END DO
3960 : END DO
3961 150 : CALL pw_pool%give_back_pw(virial_pw)
3962 150 : CALL pw_pool%give_back_pw(tmp_g)
3963 150 : CALL pw_pool%give_back_pw(rho_g)
3964 150 : DEALLOCATE (virial_pw, tmp_g, rho_g)
3965 :
3966 150 : CALL timestop(handle)
3967 :
3968 150 : END SUBROUTINE virial_laplace
3969 :
3970 : ! **************************************************************************************************
3971 : !> \brief Prepare objects for the calculation of the 2nd derivatives of the density functional.
3972 : !> The calculation must then be performed with xc_calc_2nd_deriv.
3973 : !> \param deriv_set object containing the XC derivatives (out)
3974 : !> \param rho_set object that will contain the density at which the
3975 : !> derivatives were calculated
3976 : !> \param rho_r the place where you evaluate the derivative
3977 : !> \param pw_pool the pool for the grids
3978 : !> \param xc_section which functional should be used and how to calculate it
3979 : !> \param tau_r kinetic energy density in real space
3980 : ! **************************************************************************************************
3981 12660 : SUBROUTINE xc_prep_2nd_deriv(deriv_set, &
3982 : rho_set, rho_r, pw_pool, xc_section, tau_r)
3983 :
3984 : TYPE(xc_derivative_set_type) :: deriv_set
3985 : TYPE(xc_rho_set_type) :: rho_set
3986 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
3987 : TYPE(pw_pool_type), POINTER :: pw_pool
3988 : TYPE(section_vals_type), POINTER :: xc_section
3989 : TYPE(pw_r3d_rs_type), DIMENSION(:), OPTIONAL, POINTER :: tau_r
3990 :
3991 : CHARACTER(len=*), PARAMETER :: routineN = 'xc_prep_2nd_deriv'
3992 :
3993 : INTEGER :: handle, nspins
3994 : LOGICAL :: lsd
3995 12660 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
3996 12660 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: tau
3997 :
3998 12660 : CALL timeset(routineN, handle)
3999 :
4000 12660 : CPASSERT(ASSOCIATED(xc_section))
4001 12660 : CPASSERT(ASSOCIATED(pw_pool))
4002 :
4003 12660 : nspins = SIZE(rho_r)
4004 12660 : lsd = (nspins /= 1)
4005 :
4006 12660 : NULLIFY (rho_g, tau)
4007 12660 : IF (PRESENT(tau_r)) &
4008 11996 : tau => tau_r
4009 :
4010 12660 : IF (section_get_lval(xc_section, "2ND_DERIV_ANALYTICAL")) THEN
4011 : CALL xc_rho_set_and_dset_create(rho_set, deriv_set, 2, &
4012 : rho_r, rho_g, tau, xc_section, pw_pool, &
4013 12576 : calc_potential=.TRUE.)
4014 : ELSE
4015 : CALL xc_rho_set_and_dset_create(rho_set, deriv_set, 1, &
4016 : rho_r, rho_g, tau, xc_section, pw_pool, &
4017 84 : calc_potential=.TRUE.)
4018 : END IF
4019 :
4020 12660 : CALL timestop(handle)
4021 :
4022 12660 : END SUBROUTINE xc_prep_2nd_deriv
4023 :
4024 : ! **************************************************************************************************
4025 : !> \brief Prepare deriv_set for the calculation of the 3rd derivatives of the density functional.
4026 : !> The calculation must then be performed with xc_calc_3rd_deriv.
4027 : !> \param deriv_set object containing the XC derivatives (out)
4028 : !> \param rho_set object that will contain the density at which the
4029 : !> derivatives were calculated
4030 : !> \param rho_r the place where you evaluate the derivative
4031 : !> \param pw_pool the pool for the grids
4032 : !> \param xc_section which functional should be used and how to calculate it
4033 : !> \param tau_r kinetic energy density in real space
4034 : !> \param do_sf Flag to activate the noncollinear kernel for spin flip calculations
4035 : !> \par History
4036 : !> * 07.2024 Created [LHS]
4037 : ! **************************************************************************************************
4038 0 : SUBROUTINE xc_prep_3rd_deriv(deriv_set, rho_set, rho_r, pw_pool, xc_section, tau_r, do_sf)
4039 :
4040 : TYPE(xc_derivative_set_type) :: deriv_set
4041 : TYPE(xc_rho_set_type) :: rho_set
4042 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
4043 : TYPE(pw_pool_type), POINTER :: pw_pool
4044 : TYPE(section_vals_type), POINTER :: xc_section
4045 : TYPE(pw_r3d_rs_type), DIMENSION(:), OPTIONAL, POINTER :: tau_r
4046 : LOGICAL, OPTIONAL :: do_sf
4047 :
4048 : CHARACTER(len=*), PARAMETER :: routineN = 'xc_prep_3rd_deriv'
4049 :
4050 : INTEGER :: handle, nspins
4051 : LOGICAL :: lsd, my_do_sf
4052 0 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
4053 0 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: tau
4054 :
4055 0 : CALL timeset(routineN, handle)
4056 :
4057 0 : CPASSERT(ASSOCIATED(xc_section))
4058 0 : CPASSERT(ASSOCIATED(pw_pool))
4059 :
4060 0 : nspins = SIZE(rho_r)
4061 0 : lsd = (nspins /= 1)
4062 :
4063 0 : NULLIFY (rho_g, tau)
4064 0 : IF (PRESENT(tau_r)) &
4065 0 : tau => tau_r
4066 :
4067 0 : my_do_sf = .FALSE.
4068 0 : IF (PRESENT(do_sf)) my_do_sf = do_sf
4069 :
4070 0 : IF (do_sf) THEN
4071 : CALL xc_rho_set_and_dset_create(rho_set, deriv_set, 2, &
4072 : rho_r, rho_g, tau, xc_section, pw_pool, &
4073 0 : calc_potential=.TRUE.)
4074 : ELSE
4075 : CALL xc_rho_set_and_dset_create(rho_set, deriv_set, 3, &
4076 : rho_r, rho_g, tau, xc_section, pw_pool, &
4077 0 : calc_potential=.TRUE.)
4078 : END IF
4079 :
4080 0 : CALL timestop(handle)
4081 :
4082 0 : END SUBROUTINE xc_prep_3rd_deriv
4083 :
4084 : ! **************************************************************************************************
4085 : !> \brief divides derivatives from deriv_set by norm_drho
4086 : !> \param deriv_set ...
4087 : !> \param rho_set ...
4088 : !> \param lsd ...
4089 : ! **************************************************************************************************
4090 147945 : SUBROUTINE divide_by_norm_drho(deriv_set, rho_set, lsd)
4091 :
4092 : TYPE(xc_derivative_set_type), INTENT(INOUT) :: deriv_set
4093 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
4094 : LOGICAL, INTENT(IN) :: lsd
4095 :
4096 147945 : INTEGER, DIMENSION(:), POINTER :: split_desc
4097 : INTEGER :: idesc
4098 : INTEGER, DIMENSION(2, 3) :: bo
4099 : REAL(KIND=dp) :: drho_cutoff
4100 147945 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: norm_drho, norm_drhoa, norm_drhob
4101 1775340 : TYPE(cp_3d_r_cp_type), DIMENSION(3) :: drho, drhoa, drhob
4102 : TYPE(cp_sll_xc_deriv_type), POINTER :: pos
4103 : TYPE(xc_derivative_type), POINTER :: deriv_att
4104 :
4105 : ! check for unknown derivatives and divide by norm_drho where necessary
4106 :
4107 1479450 : bo = rho_set%local_bounds
4108 : CALL xc_rho_set_get(rho_set, drho_cutoff=drho_cutoff, norm_drho=norm_drho, &
4109 : norm_drhoa=norm_drhoa, norm_drhob=norm_drhob, &
4110 147945 : drho=drho, drhoa=drhoa, drhob=drhob, can_return_null=.TRUE.)
4111 :
4112 147945 : pos => deriv_set%derivs
4113 635029 : DO WHILE (cp_sll_xc_deriv_next(pos, el_att=deriv_att))
4114 487084 : CALL xc_derivative_get(deriv_att, split_desc=split_desc)
4115 1042504 : DO idesc = 1, SIZE(split_desc)
4116 487084 : SELECT CASE (split_desc(idesc))
4117 : CASE (deriv_norm_drho)
4118 123329 : IF (ASSOCIATED(norm_drho)) THEN
4119 123329 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(deriv_att,norm_drho,drho_cutoff)
4120 : deriv_att%deriv_data(:, :, :) = deriv_att%deriv_data(:, :, :)/ &
4121 : MAX(norm_drho(:, :, :), drho_cutoff)
4122 : !$OMP END PARALLEL WORKSHARE
4123 0 : ELSE IF (ASSOCIATED(drho(1)%array)) THEN
4124 0 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(deriv_att,drho,drho_cutoff)
4125 : deriv_att%deriv_data(:, :, :) = deriv_att%deriv_data(:, :, :)/ &
4126 : MAX(SQRT(drho(1)%array(:, :, :)**2 + &
4127 : drho(2)%array(:, :, :)**2 + &
4128 : drho(3)%array(:, :, :)**2), drho_cutoff)
4129 : !$OMP END PARALLEL WORKSHARE
4130 0 : ELSE IF (ASSOCIATED(drhoa(1)%array) .AND. ASSOCIATED(drhob(1)%array)) THEN
4131 0 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(deriv_att,drhoa,drhob,drho_cutoff)
4132 : deriv_att%deriv_data(:, :, :) = deriv_att%deriv_data(:, :, :)/ &
4133 : MAX(SQRT((drhoa(1)%array(:, :, :) + drhob(1)%array(:, :, :))**2 + &
4134 : (drhoa(2)%array(:, :, :) + drhob(2)%array(:, :, :))**2 + &
4135 : (drhoa(3)%array(:, :, :) + drhob(3)%array(:, :, :))**2), drho_cutoff)
4136 : !$OMP END PARALLEL WORKSHARE
4137 : ELSE
4138 0 : CPABORT("Normalization of derivative requires any of norm_drho, drho or drhoa+drhob!")
4139 : END IF
4140 : CASE (deriv_norm_drhoa)
4141 21422 : IF (ASSOCIATED(norm_drhoa)) THEN
4142 21422 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(deriv_att,norm_drhoa,drho_cutoff)
4143 : deriv_att%deriv_data(:, :, :) = deriv_att%deriv_data(:, :, :)/ &
4144 : MAX(norm_drhoa(:, :, :), drho_cutoff)
4145 : !$OMP END PARALLEL WORKSHARE
4146 0 : ELSE IF (ASSOCIATED(drhoa(1)%array)) THEN
4147 0 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(deriv_att,drhoa,drho_cutoff)
4148 : deriv_att%deriv_data(:, :, :) = deriv_att%deriv_data(:, :, :)/ &
4149 : MAX(SQRT(drhoa(1)%array(:, :, :)**2 + &
4150 : drhoa(2)%array(:, :, :)**2 + &
4151 : drhoa(3)%array(:, :, :)**2), drho_cutoff)
4152 : !$OMP END PARALLEL WORKSHARE
4153 : ELSE
4154 0 : CPABORT("Normalization of derivative requires any of norm_drhoa or drhoa!")
4155 : END IF
4156 : CASE (deriv_norm_drhob)
4157 21418 : IF (ASSOCIATED(norm_drhob)) THEN
4158 21418 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(deriv_att,norm_drhob,drho_cutoff)
4159 : deriv_att%deriv_data(:, :, :) = deriv_att%deriv_data(:, :, :)/ &
4160 : MAX(norm_drhob(:, :, :), drho_cutoff)
4161 : !$OMP END PARALLEL WORKSHARE
4162 0 : ELSE IF (ASSOCIATED(drhob(1)%array)) THEN
4163 0 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(deriv_att,drhob,drho_cutoff)
4164 : deriv_att%deriv_data(:, :, :) = deriv_att%deriv_data(:, :, :)/ &
4165 : MAX(SQRT(drhob(1)%array(:, :, :)**2 + &
4166 : drhob(2)%array(:, :, :)**2 + &
4167 : drhob(3)%array(:, :, :)**2), drho_cutoff)
4168 : !$OMP END PARALLEL WORKSHARE
4169 : ELSE
4170 0 : CPABORT("Normalization of derivative requires any of norm_drhob or drhob!")
4171 : END IF
4172 : CASE (deriv_rho, deriv_tau, deriv_laplace_rho)
4173 167518 : IF (lsd) &
4174 0 : CPABORT(TRIM(id_to_desc(split_desc(idesc)))//" not handled in lsd!'")
4175 : CASE (deriv_rhoa, deriv_rhob, deriv_tau_a, deriv_tau_b, deriv_laplace_rhoa, deriv_laplace_rhob)
4176 : CASE default
4177 407475 : CPABORT("Unknown derivative id")
4178 : END SELECT
4179 : END DO
4180 : END DO
4181 :
4182 147945 : END SUBROUTINE divide_by_norm_drho
4183 :
4184 : ! **************************************************************************************************
4185 : !> \brief allocates and calculates drho from given spin densities drhoa, drhob
4186 : !> \param drho ...
4187 : !> \param drhoa ...
4188 : !> \param drhob ...
4189 : ! **************************************************************************************************
4190 20936 : SUBROUTINE calc_drho_from_ab(drho, drhoa, drhob)
4191 : TYPE(cp_3d_r_cp_type), DIMENSION(3), INTENT(OUT) :: drho
4192 : TYPE(cp_3d_r_cp_type), DIMENSION(3), INTENT(IN) :: drhoa, drhob
4193 :
4194 : CHARACTER(len=*), PARAMETER :: routineN = 'calc_drho_from_ab'
4195 :
4196 : INTEGER :: handle, idir
4197 :
4198 5234 : CALL timeset(routineN, handle)
4199 :
4200 20936 : DO idir = 1, 3
4201 15702 : NULLIFY (drho(idir)%array)
4202 : ALLOCATE (drho(idir)%array(LBOUND(drhoa(1)%array, 1):UBOUND(drhoa(1)%array, 1), &
4203 : LBOUND(drhoa(1)%array, 2):UBOUND(drhoa(1)%array, 2), &
4204 188424 : LBOUND(drhoa(1)%array, 3):UBOUND(drhoa(1)%array, 3)))
4205 20936 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(drho,drhoa,drhob,idir)
4206 : drho(idir)%array(:, :, :) = drhoa(idir)%array(:, :, :) + drhob(idir)%array(:, :, :)
4207 : !$OMP END PARALLEL WORKSHARE
4208 : END DO
4209 :
4210 5234 : CALL timestop(handle)
4211 :
4212 5234 : END SUBROUTINE calc_drho_from_ab
4213 :
4214 : ! **************************************************************************************************
4215 : !> \brief allocates and calculates drho from given spin densities drhoa, drhob
4216 : !> \param drho ...
4217 : !> \param drhoa ...
4218 : !> \param drhob ...
4219 : ! **************************************************************************************************
4220 1048 : SUBROUTINE calc_drho_from_a(drho, drhoa)
4221 : TYPE(cp_3d_r_cp_type), DIMENSION(3), INTENT(OUT) :: drho
4222 : TYPE(cp_3d_r_cp_type), DIMENSION(3), INTENT(IN) :: drhoa
4223 :
4224 : CHARACTER(len=*), PARAMETER :: routineN = 'calc_drho_from_a'
4225 :
4226 : INTEGER :: handle, idir
4227 :
4228 262 : CALL timeset(routineN, handle)
4229 :
4230 1048 : DO idir = 1, 3
4231 786 : NULLIFY (drho(idir)%array)
4232 : ALLOCATE (drho(idir)%array(LBOUND(drhoa(1)%array, 1):UBOUND(drhoa(1)%array, 1), &
4233 : LBOUND(drhoa(1)%array, 2):UBOUND(drhoa(1)%array, 2), &
4234 9432 : LBOUND(drhoa(1)%array, 3):UBOUND(drhoa(1)%array, 3)))
4235 1048 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(drho,drhoa,idir)
4236 : drho(idir)%array(:, :, :) = drhoa(idir)%array(:, :, :)
4237 : !$OMP END PARALLEL WORKSHARE
4238 : END DO
4239 :
4240 262 : CALL timestop(handle)
4241 :
4242 262 : END SUBROUTINE calc_drho_from_a
4243 :
4244 : ! **************************************************************************************************
4245 : !> \brief allocates and calculates dot products of two density gradients
4246 : !> \param dr1dr ...
4247 : !> \param drho ...
4248 : !> \param drho1 ...
4249 : ! **************************************************************************************************
4250 25740 : SUBROUTINE prepare_dr1dr(dr1dr, drho, drho1)
4251 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
4252 : INTENT(OUT) :: dr1dr
4253 : TYPE(cp_3d_r_cp_type), DIMENSION(3), INTENT(IN) :: drho, drho1
4254 :
4255 : CHARACTER(len=*), PARAMETER :: routineN = 'prepare_dr1dr'
4256 :
4257 : INTEGER :: handle, idir
4258 :
4259 25740 : CALL timeset(routineN, handle)
4260 :
4261 0 : ALLOCATE (dr1dr(LBOUND(drho(1)%array, 1):UBOUND(drho(1)%array, 1), &
4262 : LBOUND(drho(1)%array, 2):UBOUND(drho(1)%array, 2), &
4263 308880 : LBOUND(drho(1)%array, 3):UBOUND(drho(1)%array, 3)))
4264 :
4265 25740 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(dr1dr,drho,drho1)
4266 : dr1dr(:, :, :) = drho(1)%array(:, :, :)*drho1(1)%array(:, :, :)
4267 : !$OMP END PARALLEL WORKSHARE
4268 77220 : DO idir = 2, 3
4269 77220 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(dr1dr,drho,drho1,idir)
4270 : dr1dr(:, :, :) = dr1dr(:, :, :) + drho(idir)%array(:, :, :)*drho1(idir)%array(:, :, :)
4271 : !$OMP END PARALLEL WORKSHARE
4272 : END DO
4273 :
4274 25740 : CALL timestop(handle)
4275 :
4276 25740 : END SUBROUTINE prepare_dr1dr
4277 :
4278 : ! **************************************************************************************************
4279 : !> \brief allocates and calculates dot product of two densities for triplets
4280 : !> \param dr1dr ...
4281 : !> \param drhoa ...
4282 : !> \param drhob ...
4283 : !> \param drho1a ...
4284 : !> \param drho1b ...
4285 : !> \param fac ...
4286 : ! **************************************************************************************************
4287 710 : SUBROUTINE prepare_dr1dr_ab(dr1dr, drhoa, drhob, drho1a, drho1b, fac)
4288 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
4289 : INTENT(OUT) :: dr1dr
4290 : TYPE(cp_3d_r_cp_type), DIMENSION(3), INTENT(IN) :: drhoa, drhob, drho1a, drho1b
4291 : REAL(KIND=dp), INTENT(IN) :: fac
4292 :
4293 : CHARACTER(len=*), PARAMETER :: routineN = 'prepare_dr1dr_ab'
4294 :
4295 : INTEGER :: handle, idir
4296 :
4297 710 : CALL timeset(routineN, handle)
4298 :
4299 0 : ALLOCATE (dr1dr(LBOUND(drhoa(1)%array, 1):UBOUND(drhoa(1)%array, 1), &
4300 : LBOUND(drhoa(1)%array, 2):UBOUND(drhoa(1)%array, 2), &
4301 8520 : LBOUND(drhoa(1)%array, 3):UBOUND(drhoa(1)%array, 3)))
4302 :
4303 710 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(fac,dr1dr,drho1a,drho1b,drhoa,drhob)
4304 : dr1dr(:, :, :) = drhoa(1)%array(:, :, :)*(drho1a(1)%array(:, :, :) + &
4305 : fac*drho1b(1)%array(:, :, :)) + &
4306 : drhob(1)%array(:, :, :)*(fac*drho1a(1)%array(:, :, :) + &
4307 : drho1b(1)%array(:, :, :))
4308 : !$OMP END PARALLEL WORKSHARE
4309 2130 : DO idir = 2, 3
4310 2130 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(fac,dr1dr,drho1a,drho1b,drhoa,drhob,idir)
4311 : dr1dr(:, :, :) = dr1dr(:, :, :) + &
4312 : drhoa(idir)%array(:, :, :)*(drho1a(idir)%array(:, :, :) + &
4313 : fac*drho1b(idir)%array(:, :, :)) + &
4314 : drhob(idir)%array(:, :, :)*(fac*drho1a(idir)%array(:, :, :) + &
4315 : drho1b(idir)%array(:, :, :))
4316 : !$OMP END PARALLEL WORKSHARE
4317 : END DO
4318 :
4319 710 : CALL timestop(handle)
4320 :
4321 710 : END SUBROUTINE prepare_dr1dr_ab
4322 :
4323 : ! **************************************************************************************************
4324 : !> \brief checks for gradients
4325 : !> \param deriv_set ...
4326 : !> \param lsd ...
4327 : !> \param gradient_f ...
4328 : !> \param tau_f ...
4329 : !> \param laplace_f ...
4330 : ! **************************************************************************************************
4331 146479 : SUBROUTINE check_for_derivatives(deriv_set, lsd, rho_f, gradient_f, tau_f, laplace_f)
4332 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
4333 : LOGICAL, INTENT(IN) :: lsd
4334 : LOGICAL, INTENT(OUT) :: rho_f, gradient_f, tau_f, laplace_f
4335 :
4336 : CHARACTER(len=*), PARAMETER :: routineN = 'check_for_derivatives'
4337 :
4338 : INTEGER :: handle, iorder, order
4339 146479 : INTEGER, DIMENSION(:), POINTER :: split_desc
4340 : TYPE(cp_sll_xc_deriv_type), POINTER :: pos
4341 : TYPE(xc_derivative_type), POINTER :: deriv_att
4342 :
4343 146479 : CALL timeset(routineN, handle)
4344 :
4345 146479 : rho_f = .FALSE.
4346 146479 : gradient_f = .FALSE.
4347 146479 : tau_f = .FALSE.
4348 146479 : laplace_f = .FALSE.
4349 : ! check for unknown derivatives
4350 146479 : pos => deriv_set%derivs
4351 676273 : DO WHILE (cp_sll_xc_deriv_next(pos, el_att=deriv_att))
4352 : CALL xc_derivative_get(deriv_att, order=order, &
4353 529794 : split_desc=split_desc)
4354 676273 : IF (lsd) THEN
4355 340529 : DO iorder = 1, size(split_desc)
4356 166793 : SELECT CASE (split_desc(iorder))
4357 : CASE (deriv_rhoa, deriv_rhob)
4358 89684 : rho_f = .TRUE.
4359 : CASE (deriv_norm_drho, deriv_norm_drhoa, deriv_norm_drhob)
4360 79980 : gradient_f = .TRUE.
4361 : CASE (deriv_tau_a, deriv_tau_b)
4362 2760 : tau_f = .TRUE.
4363 : CASE (deriv_laplace_rhoa, deriv_laplace_rhob)
4364 1312 : laplace_f = .TRUE.
4365 : CASE (deriv_rho, deriv_tau, deriv_laplace_rho)
4366 0 : CPABORT("Derivative not handled in lsd!")
4367 : CASE default
4368 173736 : CPABORT("Unknown derivative id")
4369 : END SELECT
4370 : END DO
4371 : ELSE
4372 674192 : DO iorder = 1, size(split_desc)
4373 363001 : SELECT CASE (split_desc(iorder))
4374 : CASE (deriv_rho)
4375 183770 : rho_f = .TRUE.
4376 : CASE (deriv_tau)
4377 4996 : tau_f = .TRUE.
4378 : CASE (deriv_norm_drho)
4379 120987 : gradient_f = .TRUE.
4380 : CASE (deriv_laplace_rho)
4381 1438 : laplace_f = .TRUE.
4382 : CASE default
4383 311191 : CPABORT("Unknown derivative id")
4384 : END SELECT
4385 : END DO
4386 : END IF
4387 : END DO
4388 :
4389 146479 : CALL timestop(handle)
4390 :
4391 146479 : END SUBROUTINE check_for_derivatives
4392 :
4393 : END MODULE xc
|
