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