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