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