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