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