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 : MODULE xc_atom
10 :
11 : USE cp_linked_list_xc_deriv, ONLY: cp_sll_xc_deriv_next,&
12 : cp_sll_xc_deriv_type
13 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
14 : section_vals_type
15 : USE kinds, ONLY: dp
16 : USE pw_pool_types, ONLY: pw_pool_type
17 : USE pw_types, ONLY: pw_r3d_rs_type
18 : USE xc, ONLY: divide_by_norm_drho,&
19 : xc_calc_2nd_deriv_analytical
20 : USE xc_derivative_desc, ONLY: &
21 : deriv_norm_drho, deriv_norm_drhoa, deriv_norm_drhob, deriv_rho, deriv_rhoa, deriv_rhob, &
22 : deriv_tau, deriv_tau_a, deriv_tau_b
23 : USE xc_derivative_set_types, ONLY: xc_derivative_set_type,&
24 : xc_dset_get_derivative
25 : USE xc_derivative_types, ONLY: xc_derivative_get,&
26 : xc_derivative_type
27 : USE xc_derivatives, ONLY: xc_functionals_eval
28 : USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
29 : USE xc_rho_set_types, ONLY: xc_rho_set_get,&
30 : xc_rho_set_type
31 : #include "../base/base_uses.f90"
32 :
33 : IMPLICIT NONE
34 :
35 : PRIVATE
36 :
37 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_atom'
38 :
39 : PUBLIC :: vxc_of_r_new, vxc_of_r_epr, xc_rho_set_atom_update, xc_2nd_deriv_of_r, fill_rho_set
40 :
41 : CONTAINS
42 :
43 : ! **************************************************************************************************
44 : !> \brief ...
45 : !> \param xc_fun_section ...
46 : !> \param rho_set ...
47 : !> \param deriv_set ...
48 : !> \param deriv_order ...
49 : !> \param needs ...
50 : !> \param w ...
51 : !> \param lsd ...
52 : !> \param na ...
53 : !> \param nr ...
54 : !> \param exc ...
55 : !> \param vxc ...
56 : !> \param vxg ...
57 : !> \param vtau ...
58 : !> \param energy_only ...
59 : !> \param adiabatic_rescale_factor ...
60 : ! **************************************************************************************************
61 77274 : SUBROUTINE vxc_of_r_new(xc_fun_section, rho_set, deriv_set, deriv_order, needs, w, &
62 : lsd, na, nr, exc, vxc, vxg, vtau, &
63 : energy_only, adiabatic_rescale_factor)
64 :
65 : ! This routine updates rho_set by giving to it the rho and drho that are needed.
66 : ! Since for the local densities rho1_h and rho1_s local grids are used it is not possible
67 : ! to call xc_rho_set_update.
68 : ! As input of this routine one gets rho and drho on a one dimensional grid.
69 : ! The grid is the angular grid corresponding to a given point ir_pnt on the radial grid.
70 : ! The derivatives are calculated on this one dimensional grid, the results are stored in
71 : ! exc, vxc(1:na,ir_pnt,ispin), vxg(1:na,ir_pnt,ispin), vxg_cross(1:na,ir_pnt,ispin)
72 : ! Afterwords the arrays containing the derivatives are put to zero so that the routine
73 : ! can safely be called for the next radial point ir_pnt
74 :
75 : TYPE(section_vals_type), POINTER :: xc_fun_section
76 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
77 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
78 : INTEGER, INTENT(in) :: deriv_order
79 : TYPE(xc_rho_cflags_type), INTENT(IN) :: needs
80 : REAL(dp), DIMENSION(:, :), INTENT(IN) :: w
81 : LOGICAL, INTENT(IN) :: lsd
82 : INTEGER, INTENT(in) :: na, nr
83 : REAL(dp) :: exc
84 : REAL(dp), DIMENSION(:, :, :), POINTER :: vxc
85 : REAL(dp), DIMENSION(:, :, :, :), POINTER :: vxg
86 : REAL(dp), DIMENSION(:, :, :), POINTER :: vtau
87 : LOGICAL, INTENT(IN), OPTIONAL :: energy_only
88 : REAL(dp), INTENT(IN), OPTIONAL :: adiabatic_rescale_factor
89 :
90 : CHARACTER(LEN=*), PARAMETER :: routineN = 'vxc_of_r_new'
91 :
92 : INTEGER :: handle, ia, idir, ir
93 : LOGICAL :: gradient_f, my_only_energy
94 : REAL(dp) :: my_adiabatic_rescale_factor
95 77274 : REAL(dp), DIMENSION(:, :, :), POINTER :: deriv_data
96 : REAL(KIND=dp) :: drho_cutoff
97 : TYPE(xc_derivative_type), POINTER :: deriv_att
98 :
99 77274 : CALL timeset(routineN, handle)
100 77274 : my_only_energy = .FALSE.
101 77274 : IF (PRESENT(energy_only)) my_only_energy = energy_only
102 :
103 77274 : IF (PRESENT(adiabatic_rescale_factor)) THEN
104 77274 : my_adiabatic_rescale_factor = adiabatic_rescale_factor
105 : ELSE
106 : my_adiabatic_rescale_factor = 1.0_dp
107 : END IF
108 :
109 : gradient_f = (needs%drho_spin .OR. needs%norm_drho_spin .OR. &
110 77274 : needs%drho .OR. needs%norm_drho)
111 :
112 : ! Calculate the derivatives
113 : CALL xc_functionals_eval(xc_fun_section, &
114 : lsd=lsd, &
115 : rho_set=rho_set, &
116 : deriv_set=deriv_set, &
117 77274 : deriv_order=deriv_order)
118 :
119 77274 : CALL xc_rho_set_get(rho_set, drho_cutoff=drho_cutoff)
120 :
121 77274 : NULLIFY (deriv_data)
122 :
123 : ! EXC energy
124 77274 : deriv_att => xc_dset_get_derivative(deriv_set, [INTEGER::])
125 77274 : exc = 0.0_dp
126 77274 : IF (ASSOCIATED(deriv_att)) THEN
127 77202 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
128 4278382 : DO ir = 1, nr
129 214505862 : DO ia = 1, na
130 214428660 : exc = exc + deriv_data(ia, ir, 1)*w(ia, ir)
131 : END DO
132 : END DO
133 77202 : NULLIFY (deriv_data)
134 : END IF
135 : ! Calculate the potential only if needed
136 77274 : IF (.NOT. my_only_energy) THEN
137 : ! Derivative with respect to the density
138 74020 : IF (lsd) THEN
139 9852 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_rhoa])
140 9852 : IF (ASSOCIATED(deriv_att)) THEN
141 9848 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
142 62775616 : vxc(:, :, 1) = deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
143 9848 : NULLIFY (deriv_data)
144 : END IF
145 9852 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_rhob])
146 9852 : IF (ASSOCIATED(deriv_att)) THEN
147 9848 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
148 62775616 : vxc(:, :, 2) = deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
149 9848 : NULLIFY (deriv_data)
150 : END IF
151 9852 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_rho])
152 9852 : IF (ASSOCIATED(deriv_att)) THEN
153 0 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
154 0 : vxc(:, :, 1) = vxc(:, :, 1) + deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
155 0 : vxc(:, :, 2) = vxc(:, :, 2) + deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
156 0 : NULLIFY (deriv_data)
157 : END IF
158 : ELSE
159 64168 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_rho])
160 64168 : IF (ASSOCIATED(deriv_att)) THEN
161 64100 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
162 348247000 : vxc(:, :, 1) = deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
163 64100 : NULLIFY (deriv_data)
164 : END IF
165 : END IF ! lsd
166 :
167 : ! Derivatives with respect to the gradient
168 74020 : IF (lsd) THEN
169 9852 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa])
170 9852 : IF (ASSOCIATED(deriv_att)) THEN
171 6344 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
172 : !$OMP PARALLEL DO DEFAULT(NONE) COLLAPSE(2) SCHEDULE(STATIC) &
173 : !$OMP PRIVATE(ia, idir, ir) &
174 6344 : !$OMP SHARED(deriv_data, drho_cutoff, my_adiabatic_rescale_factor, na, nr, rho_set, vxg, w)
175 : DO ir = 1, nr
176 : DO ia = 1, na
177 : DO idir = 1, 3
178 : IF (rho_set%norm_drhoa(ia, ir, 1) > drho_cutoff) THEN
179 : vxg(idir, ia, ir, 1) = rho_set%drhoa(idir)%array(ia, ir, 1)* &
180 : deriv_data(ia, ir, 1)*w(ia, ir)/ &
181 : rho_set%norm_drhoa(ia, ir, 1)*my_adiabatic_rescale_factor
182 : ELSE
183 : vxg(idir, ia, ir, 1) = 0.0_dp
184 : END IF
185 : END DO
186 : END DO
187 : END DO
188 : !$OMP END PARALLEL DO
189 6344 : NULLIFY (deriv_data)
190 : END IF
191 9852 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob])
192 9852 : IF (ASSOCIATED(deriv_att)) THEN
193 6344 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
194 : !$OMP PARALLEL DO DEFAULT(NONE) COLLAPSE(2) SCHEDULE(STATIC) &
195 : !$OMP PRIVATE(ia, idir, ir) &
196 6344 : !$OMP SHARED(deriv_data, drho_cutoff, my_adiabatic_rescale_factor, na, nr, rho_set, vxg, w)
197 : DO ir = 1, nr
198 : DO ia = 1, na
199 : DO idir = 1, 3
200 : IF (rho_set%norm_drhob(ia, ir, 1) > drho_cutoff) THEN
201 : vxg(idir, ia, ir, 2) = rho_set%drhob(idir)%array(ia, ir, 1)* &
202 : deriv_data(ia, ir, 1)*w(ia, ir)/ &
203 : rho_set%norm_drhob(ia, ir, 1)*my_adiabatic_rescale_factor
204 : ELSE
205 : vxg(idir, ia, ir, 2) = 0.0_dp
206 : END IF
207 : END DO
208 : END DO
209 : END DO
210 : !$OMP END PARALLEL DO
211 6344 : NULLIFY (deriv_data)
212 : END IF
213 : ! Cross Terms
214 9852 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_norm_drho])
215 9852 : IF (ASSOCIATED(deriv_att)) THEN
216 6056 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
217 : !$OMP PARALLEL DO DEFAULT(NONE) COLLAPSE(2) SCHEDULE(STATIC) &
218 : !$OMP PRIVATE(ia, idir, ir) &
219 6056 : !$OMP SHARED(deriv_data, drho_cutoff, my_adiabatic_rescale_factor, na, nr, rho_set, vxg, w)
220 : DO ir = 1, nr
221 : DO ia = 1, na
222 : DO idir = 1, 3
223 : IF (rho_set%norm_drho(ia, ir, 1) > drho_cutoff) THEN
224 : vxg(idir, ia, ir, 1:2) = &
225 : vxg(idir, ia, ir, 1:2) + ( &
226 : rho_set%drhoa(idir)%array(ia, ir, 1) + &
227 : rho_set%drhob(idir)%array(ia, ir, 1))* &
228 : deriv_data(ia, ir, 1)*w(ia, ir)/rho_set%norm_drho(ia, ir, 1)* &
229 : my_adiabatic_rescale_factor
230 : END IF
231 : END DO
232 : END DO
233 : END DO
234 : !$OMP END PARALLEL DO
235 6056 : NULLIFY (deriv_data)
236 : END IF
237 : ELSE
238 64168 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_norm_drho])
239 64168 : IF (ASSOCIATED(deriv_att)) THEN
240 40964 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
241 : !$OMP PARALLEL DO DEFAULT(NONE) COLLAPSE(2) SCHEDULE(STATIC) &
242 : !$OMP PRIVATE(ia, idir, ir) &
243 40964 : !$OMP SHARED(deriv_data, drho_cutoff, my_adiabatic_rescale_factor, na, nr, rho_set, vxg, w)
244 : DO ir = 1, nr
245 : DO ia = 1, na
246 : IF (rho_set%norm_drho(ia, ir, 1) > drho_cutoff) THEN
247 : DO idir = 1, 3
248 : vxg(idir, ia, ir, 1) = rho_set%drho(idir)%array(ia, ir, 1)* &
249 : deriv_data(ia, ir, 1)*w(ia, ir)/ &
250 : rho_set%norm_drho(ia, ir, 1)*my_adiabatic_rescale_factor
251 : END DO
252 : ELSE
253 : vxg(1:3, ia, ir, 1) = 0.0_dp
254 : END IF
255 : END DO
256 : END DO
257 : !$OMP END PARALLEL DO
258 40964 : NULLIFY (deriv_data)
259 : END IF
260 : END IF ! lsd
261 : ! Derivative with respect to tau
262 74020 : IF (lsd) THEN
263 9852 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_tau_a])
264 9852 : IF (ASSOCIATED(deriv_att)) THEN
265 16 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
266 81632 : vtau(:, :, 1) = deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
267 16 : NULLIFY (deriv_data)
268 : END IF
269 9852 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_tau_b])
270 9852 : IF (ASSOCIATED(deriv_att)) THEN
271 16 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
272 81632 : vtau(:, :, 2) = deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
273 16 : NULLIFY (deriv_data)
274 : END IF
275 9852 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_tau])
276 9852 : IF (ASSOCIATED(deriv_att)) THEN
277 0 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
278 0 : vtau(:, :, 1) = vtau(:, :, 1) + deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
279 0 : vtau(:, :, 2) = vtau(:, :, 2) + deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
280 0 : NULLIFY (deriv_data)
281 : END IF
282 : ELSE
283 64168 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_tau])
284 64168 : IF (ASSOCIATED(deriv_att)) THEN
285 1604 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
286 8823208 : vtau(:, :, 1) = deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
287 1604 : NULLIFY (deriv_data)
288 : END IF
289 : END IF ! lsd
290 : END IF ! only_energy
291 :
292 77274 : CALL timestop(handle)
293 :
294 77274 : END SUBROUTINE vxc_of_r_new
295 :
296 : ! **************************************************************************************************
297 : !> \brief Specific EPR version of vxc_of_r_new
298 : !> \param xc_fun_section ...
299 : !> \param rho_set ...
300 : !> \param deriv_set ...
301 : !> \param needs ...
302 : !> \param w ...
303 : !> \param lsd ...
304 : !> \param na ...
305 : !> \param nr ...
306 : !> \param exc ...
307 : !> \param vxc ...
308 : !> \param vxg ...
309 : !> \param vtau ...
310 : ! **************************************************************************************************
311 30 : SUBROUTINE vxc_of_r_epr(xc_fun_section, rho_set, deriv_set, needs, w, &
312 : lsd, na, nr, exc, vxc, vxg, vtau)
313 :
314 : TYPE(section_vals_type), POINTER :: xc_fun_section
315 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
316 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
317 : TYPE(xc_rho_cflags_type), INTENT(IN) :: needs
318 : REAL(dp), DIMENSION(:, :), INTENT(IN) :: w
319 : LOGICAL, INTENT(IN) :: lsd
320 : INTEGER, INTENT(in) :: na, nr
321 : REAL(dp) :: exc
322 : REAL(dp), DIMENSION(:, :, :), POINTER :: vxc
323 : REAL(dp), DIMENSION(:, :, :, :), POINTER :: vxg
324 : REAL(dp), DIMENSION(:, :, :), POINTER :: vtau
325 :
326 : CHARACTER(LEN=*), PARAMETER :: routineN = 'vxc_of_r_epr'
327 :
328 : INTEGER :: handle, ia, idir, ir, my_deriv_order
329 : LOGICAL :: gradient_f
330 : REAL(dp) :: my_adiabatic_rescale_factor
331 30 : REAL(dp), DIMENSION(:, :, :), POINTER :: deriv_data
332 : REAL(KIND=dp) :: drho_cutoff
333 : TYPE(xc_derivative_type), POINTER :: deriv_att
334 :
335 30 : CALL timeset(routineN, handle)
336 :
337 : MARK_USED(vxc)
338 : MARK_USED(vtau)
339 :
340 30 : my_adiabatic_rescale_factor = 1.0_dp
341 30 : my_deriv_order = 2
342 :
343 : gradient_f = (needs%drho_spin .OR. needs%norm_drho_spin .OR. &
344 30 : needs%drho .OR. needs%norm_drho)
345 :
346 : ! Calculate the derivatives
347 : CALL xc_functionals_eval(xc_fun_section, &
348 : lsd=lsd, &
349 : rho_set=rho_set, &
350 : deriv_set=deriv_set, &
351 30 : deriv_order=my_deriv_order)
352 :
353 30 : CALL xc_rho_set_get(rho_set, drho_cutoff=drho_cutoff)
354 :
355 30 : NULLIFY (deriv_data)
356 :
357 : ! nabla v_xc (using the vxg arrays)
358 : ! there's no point doing this when lsd = false
359 30 : IF (lsd) THEN
360 30 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa])
361 30 : IF (ASSOCIATED(deriv_att)) THEN
362 30 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
363 1530 : DO ir = 1, nr
364 76530 : DO ia = 1, na
365 301500 : DO idir = 1, 3
366 : vxg(idir, ia, ir, 1) = rho_set%drhoa(idir)%array(ia, ir, 1)* &
367 300000 : deriv_data(ia, ir, 1)
368 : END DO !idir
369 : END DO !ia
370 : END DO !ir
371 30 : NULLIFY (deriv_data)
372 : END IF
373 30 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob])
374 30 : IF (ASSOCIATED(deriv_att)) THEN
375 30 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
376 1530 : DO ir = 1, nr
377 76530 : DO ia = 1, na
378 301500 : DO idir = 1, 3
379 : vxg(idir, ia, ir, 2) = rho_set%drhob(idir)%array(ia, ir, 1)* &
380 300000 : deriv_data(ia, ir, 1)
381 : END DO !idir
382 : END DO !ia
383 : END DO !ir
384 30 : NULLIFY (deriv_data)
385 : END IF
386 : END IF
387 : ! EXC energy ! is that needed for epr?
388 30 : deriv_att => xc_dset_get_derivative(deriv_set, [INTEGER::])
389 30 : exc = 0.0_dp
390 30 : IF (ASSOCIATED(deriv_att)) THEN
391 30 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
392 1530 : DO ir = 1, nr
393 76530 : DO ia = 1, na
394 76500 : exc = exc + deriv_data(ia, ir, 1)*w(ia, ir)
395 : END DO
396 : END DO
397 30 : NULLIFY (deriv_data)
398 : END IF
399 :
400 30 : CALL timestop(handle)
401 :
402 30 : END SUBROUTINE vxc_of_r_epr
403 :
404 : ! **************************************************************************************************
405 : !> \brief ...
406 : !> \param rho_set ...
407 : !> \param rho1_set ...
408 : !> \param xc_section ...
409 : !> \param deriv_set ...
410 : !> \param w ...
411 : !> \param vxc ...
412 : !> \param vxg ...
413 : !> \param vtau ...
414 : !> \param do_triplet ...
415 : !> \param do_sf ...
416 : ! **************************************************************************************************
417 20732 : SUBROUTINE xc_2nd_deriv_of_r(rho_set, rho1_set, xc_section, &
418 20732 : deriv_set, w, vxc, vxg, vtau, do_triplet, do_sf)
419 :
420 : ! As input of this routine one gets rho and drho on a one dimensional grid.
421 : ! The grid is the angular grid corresponding to a given point ir on the radial grid.
422 : ! The derivatives are calculated on this one dimensional grid, the results are stored in
423 : ! vxc(1:na,ir,ispin), vxg(1:na,ir,ispin), vxg_cross(1:na,ir,ispin)
424 : ! Afterwords the arrays containing the derivatives are put to zero so that the routine
425 : ! can safely be called for the next radial point ir
426 :
427 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set, rho1_set
428 : TYPE(section_vals_type), POINTER :: xc_section
429 : TYPE(xc_derivative_set_type), INTENT(INOUT) :: deriv_set
430 : REAL(dp), DIMENSION(:, :), INTENT(IN) :: w
431 : REAL(dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: vxc
432 : REAL(dp), DIMENSION(:, :, :, :), POINTER :: vxg
433 : REAL(dp), CONTIGUOUS, DIMENSION(:, :, :), &
434 : OPTIONAL, POINTER :: vtau
435 : LOGICAL, INTENT(IN), OPTIONAL :: do_triplet, do_sf
436 :
437 : CHARACTER(LEN=*), PARAMETER :: routineN = 'xc_2nd_deriv_of_r'
438 :
439 : INTEGER :: handle, ispin, nspins
440 : LOGICAL :: lsd, my_do_sf
441 : REAL(dp) :: drho_cutoff, my_fac_triplet
442 : TYPE(cp_sll_xc_deriv_type), POINTER :: pos
443 : TYPE(pw_pool_type), POINTER :: pw_pool
444 20732 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: vxc_pw, vxc_tau_pw
445 : TYPE(section_vals_type), POINTER :: xc_fun_section
446 : TYPE(xc_derivative_type), POINTER :: deriv_att
447 :
448 20732 : CALL timeset(routineN, handle)
449 :
450 20732 : nspins = SIZE(vxc, 3)
451 20732 : lsd = (nspins == 2)
452 20732 : IF (ASSOCIATED(rho_set%rhoa)) THEN
453 1514 : lsd = .TRUE.
454 : END IF
455 20732 : my_fac_triplet = 1.0_dp
456 20732 : IF (PRESENT(do_triplet)) THEN
457 11184 : IF (do_triplet) my_fac_triplet = -1.0_dp
458 : END IF
459 :
460 20732 : my_do_sf = .FALSE.
461 20732 : IF (PRESENT(do_sf)) my_do_sf = do_sf
462 :
463 20732 : CALL xc_rho_set_get(rho_set, drho_cutoff=drho_cutoff)
464 : xc_fun_section => section_vals_get_subs_vals(xc_section, &
465 20732 : "XC_FUNCTIONAL")
466 :
467 : ! Calculate the derivatives
468 : CALL xc_functionals_eval(xc_fun_section, &
469 : lsd=lsd, &
470 : rho_set=rho_set, &
471 : deriv_set=deriv_set, &
472 20732 : deriv_order=2)
473 :
474 20732 : CALL divide_by_norm_drho(deriv_set, rho_set, lsd)
475 :
476 : ! multiply by w
477 20732 : pos => deriv_set%derivs
478 137552 : DO WHILE (cp_sll_xc_deriv_next(pos, el_att=deriv_att))
479 298028552 : deriv_att%deriv_data(:, :, 1) = w(:, :)*deriv_att%deriv_data(:, :, 1)
480 : END DO
481 :
482 20732 : NULLIFY (pw_pool)
483 83428 : ALLOCATE (vxc_pw(nspins))
484 41964 : DO ispin = 1, nspins
485 41964 : vxc_pw(ispin)%array => vxc(:, :, ispin:ispin)
486 : END DO
487 :
488 20732 : NULLIFY (vxc_tau_pw)
489 20732 : IF (PRESENT(vtau)) THEN
490 20732 : IF (ASSOCIATED(vtau)) THEN
491 0 : ALLOCATE (vxc_tau_pw(nspins))
492 0 : DO ispin = 1, nspins
493 0 : vxc_tau_pw(ispin)%array => vtau(:, :, ispin:ispin)
494 : END DO
495 : END IF
496 : END IF
497 :
498 : CALL xc_calc_2nd_deriv_analytical(vxc_pw, vxc_tau_pw, deriv_set, rho_set, rho1_set, pw_pool, &
499 : xc_section, gapw=.TRUE., vxg=vxg, &
500 20732 : tddfpt_fac=my_fac_triplet, spinflip=do_sf)
501 :
502 20732 : DEALLOCATE (vxc_pw)
503 20732 : IF (ASSOCIATED(vxc_tau_pw)) DEALLOCATE (vxc_tau_pw)
504 :
505 : ! zero the derivative data for the next call
506 20732 : pos => deriv_set%derivs
507 137552 : DO WHILE (cp_sll_xc_deriv_next(pos, el_att=deriv_att))
508 298145372 : deriv_att%deriv_data = 0.0_dp
509 : END DO
510 :
511 20732 : CALL timestop(handle)
512 :
513 41464 : END SUBROUTINE xc_2nd_deriv_of_r
514 :
515 : ! **************************************************************************************************
516 : !> \brief ...
517 : !> \param rho_set ...
518 : !> \param needs ...
519 : !> \param nspins ...
520 : !> \param bo ...
521 : ! **************************************************************************************************
522 175197 : SUBROUTINE xc_rho_set_atom_update(rho_set, needs, nspins, bo)
523 :
524 : ! This routine allocates the storage arrays for rho and drho
525 : ! In calculate_vxc_atom this is called once for each atomic_kind,
526 : ! After the loop over all the atoms of the kind and over all the points
527 : ! of the radial grid for each atom, rho_set is deallocated.
528 : ! Within the same kind, at each new point on the radial grid, the rho_set
529 : ! arrays rho and drho are overwritten.
530 :
531 : TYPE(xc_rho_set_type), INTENT(INOUT) :: rho_set
532 : TYPE(xc_rho_cflags_type), INTENT(IN) :: needs
533 : INTEGER, INTENT(IN) :: nspins
534 : INTEGER, DIMENSION(2, 3), INTENT(IN) :: bo
535 :
536 : INTEGER :: idir
537 :
538 331829 : SELECT CASE (nspins)
539 : CASE (1)
540 : ! What is this for?
541 156632 : IF (needs%rho_1_3) THEN
542 4031 : NULLIFY (rho_set%rho_1_3)
543 20155 : ALLOCATE (rho_set%rho_1_3(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
544 4031 : rho_set%owns%rho_1_3 = .TRUE.
545 4031 : rho_set%has%rho_1_3 = .FALSE.
546 : END IF
547 : ! Allocate the storage space for the density
548 156632 : IF (needs%rho) THEN
549 156632 : NULLIFY (rho_set%rho)
550 783160 : ALLOCATE (rho_set%rho(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
551 156632 : rho_set%owns%rho = .TRUE.
552 156632 : rho_set%has%rho = .FALSE.
553 : END IF
554 : ! Allocate the storage space for the norm of the gradient of the density
555 156632 : IF (needs%norm_drho) THEN
556 113650 : NULLIFY (rho_set%norm_drho)
557 568250 : ALLOCATE (rho_set%norm_drho(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
558 113650 : rho_set%owns%norm_drho = .TRUE.
559 113650 : rho_set%has%norm_drho = .FALSE.
560 : END IF
561 : ! Allocate the storage space for the three components of the gradient of the density
562 156632 : IF (needs%drho) THEN
563 369200 : DO idir = 1, 3
564 276900 : NULLIFY (rho_set%drho(idir)%array)
565 1476800 : ALLOCATE (rho_set%drho(idir)%array(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
566 : END DO
567 92300 : rho_set%owns%drho = .TRUE.
568 92300 : rho_set%has%drho = .FALSE.
569 : END IF
570 : CASE (2)
571 : ! Allocate the storage space for the total density
572 18565 : IF (needs%rho) THEN
573 : ! this should never be the case unless you use LDA functionals with LSD
574 0 : NULLIFY (rho_set%rho)
575 0 : ALLOCATE (rho_set%rho(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
576 0 : rho_set%owns%rho = .TRUE.
577 0 : rho_set%has%rho = .FALSE.
578 : END IF
579 : ! What is this for?
580 18565 : IF (needs%rho_1_3) THEN
581 0 : NULLIFY (rho_set%rho_1_3)
582 0 : ALLOCATE (rho_set%rho_1_3(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
583 0 : rho_set%owns%rho_1_3 = .TRUE.
584 0 : rho_set%has%rho_1_3 = .FALSE.
585 : END IF
586 : ! What is this for?
587 18565 : IF (needs%rho_spin_1_3) THEN
588 2440 : NULLIFY (rho_set%rhoa_1_3, rho_set%rhob_1_3)
589 12200 : ALLOCATE (rho_set%rhoa_1_3(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
590 9760 : ALLOCATE (rho_set%rhob_1_3(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
591 2440 : rho_set%owns%rho_spin_1_3 = .TRUE.
592 2440 : rho_set%has%rho_spin_1_3 = .FALSE.
593 : END IF
594 : ! Allocate the storage space for the spin densities rhoa and rhob
595 18565 : IF (needs%rho_spin) THEN
596 18565 : NULLIFY (rho_set%rhoa, rho_set%rhob)
597 92825 : ALLOCATE (rho_set%rhoa(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
598 74260 : ALLOCATE (rho_set%rhob(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
599 18565 : rho_set%owns%rho_spin = .TRUE.
600 18565 : rho_set%has%rho_spin = .FALSE.
601 : END IF
602 : ! Allocate the storage space for the norm of the gradient of the total density
603 18565 : IF (needs%norm_drho) THEN
604 12259 : NULLIFY (rho_set%norm_drho)
605 61295 : ALLOCATE (rho_set%norm_drho(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
606 12259 : rho_set%owns%norm_drho = .TRUE.
607 12259 : rho_set%has%norm_drho = .FALSE.
608 : END IF
609 : ! Allocate the storage space for the norm of the gradient of rhoa and of rhob separatedly
610 18565 : IF (needs%norm_drho_spin) THEN
611 12547 : NULLIFY (rho_set%norm_drhoa, rho_set%norm_drhob)
612 62735 : ALLOCATE (rho_set%norm_drhoa(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
613 50188 : ALLOCATE (rho_set%norm_drhob(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
614 12547 : rho_set%owns%norm_drho_spin = .TRUE.
615 12547 : rho_set%has%norm_drho_spin = .FALSE.
616 : END IF
617 : ! Allocate the storage space for the components of the gradient for the total rho
618 18565 : IF (needs%drho) THEN
619 0 : DO idir = 1, 3
620 0 : NULLIFY (rho_set%drho(idir)%array)
621 0 : ALLOCATE (rho_set%drho(idir)%array(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
622 : END DO
623 0 : rho_set%owns%drho = .TRUE.
624 0 : rho_set%has%drho = .FALSE.
625 : END IF
626 : ! Allocate the storage space for the components of the gradient for rhoa and rhob
627 193762 : IF (needs%drho_spin) THEN
628 46672 : DO idir = 1, 3
629 35004 : NULLIFY (rho_set%drhoa(idir)%array, rho_set%drhob(idir)%array)
630 175020 : ALLOCATE (rho_set%drhoa(idir)%array(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
631 151684 : ALLOCATE (rho_set%drhob(idir)%array(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
632 : END DO
633 11668 : rho_set%owns%drho_spin = .TRUE.
634 11668 : rho_set%has%drho_spin = .FALSE.
635 : END IF
636 : !
637 : END SELECT
638 :
639 : ! tau part
640 175197 : IF (needs%tau) THEN
641 2728 : NULLIFY (rho_set%tau)
642 13640 : ALLOCATE (rho_set%tau(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
643 2728 : rho_set%owns%tau = .TRUE.
644 : END IF
645 175197 : IF (needs%tau_spin) THEN
646 34 : NULLIFY (rho_set%tau_a, rho_set%tau_b)
647 170 : ALLOCATE (rho_set%tau_a(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
648 136 : ALLOCATE (rho_set%tau_b(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
649 34 : rho_set%owns%tau_spin = .TRUE.
650 34 : rho_set%has%tau_spin = .FALSE.
651 : END IF
652 :
653 : ! Laplace part
654 175197 : IF (needs%laplace_rho) THEN
655 0 : NULLIFY (rho_set%laplace_rho)
656 0 : ALLOCATE (rho_set%laplace_rho(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
657 0 : rho_set%owns%laplace_rho = .TRUE.
658 : END IF
659 175197 : IF (needs%laplace_rho_spin) THEN
660 0 : NULLIFY (rho_set%laplace_rhoa)
661 0 : NULLIFY (rho_set%laplace_rhob)
662 0 : ALLOCATE (rho_set%laplace_rhoa(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
663 0 : ALLOCATE (rho_set%laplace_rhob(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
664 0 : rho_set%owns%laplace_rho_spin = .TRUE.
665 0 : rho_set%has%laplace_rho_spin = .TRUE.
666 : END IF
667 :
668 175197 : END SUBROUTINE xc_rho_set_atom_update
669 :
670 : ! **************************************************************************************************
671 : !> \brief ...
672 : !> \param rho_set ...
673 : !> \param lsd ...
674 : !> \param nspins ...
675 : !> \param needs ...
676 : !> \param rho ...
677 : !> \param drho ...
678 : !> \param tau ...
679 : !> \param na ...
680 : !> \param ir ...
681 : ! **************************************************************************************************
682 6177080 : SUBROUTINE fill_rho_set(rho_set, lsd, nspins, needs, rho, drho, tau, na, ir)
683 :
684 : TYPE(xc_rho_set_type), INTENT(INOUT) :: rho_set
685 : LOGICAL, INTENT(IN) :: lsd
686 : INTEGER, INTENT(IN) :: nspins
687 : TYPE(xc_rho_cflags_type), INTENT(IN) :: needs
688 : REAL(dp), DIMENSION(:, :, :), INTENT(IN) :: rho
689 : REAL(dp), DIMENSION(:, :, :, :), INTENT(IN) :: drho
690 : REAL(dp), DIMENSION(:, :, :), INTENT(IN) :: tau
691 : INTEGER, INTENT(IN) :: na, ir
692 :
693 : REAL(KIND=dp), PARAMETER :: f13 = (1.0_dp/3.0_dp)
694 :
695 : INTEGER :: ia, idir, my_nspins
696 : LOGICAL :: gradient_f, tddft_split
697 :
698 6177080 : my_nspins = nspins
699 6177080 : tddft_split = .FALSE.
700 6177080 : IF (lsd .AND. nspins == 1) THEN
701 90600 : my_nspins = 2
702 90600 : tddft_split = .TRUE.
703 : END IF
704 :
705 : ! some checks
706 6177080 : IF (lsd) THEN
707 : ELSE
708 5348100 : CPASSERT(SIZE(rho, 3) == 1)
709 : END IF
710 5348100 : SELECT CASE (my_nspins)
711 : CASE (1)
712 5348100 : CPASSERT(.NOT. needs%rho_spin)
713 5348100 : CPASSERT(.NOT. needs%drho_spin)
714 5348100 : CPASSERT(.NOT. needs%norm_drho_spin)
715 5348100 : CPASSERT(.NOT. needs%rho_spin_1_3)
716 : CASE (2)
717 : CASE default
718 6177080 : CPABORT("Unsupported number of spins")
719 : END SELECT
720 :
721 : gradient_f = (needs%drho_spin .OR. needs%norm_drho_spin .OR. &
722 6177080 : needs%drho .OR. needs%norm_drho)
723 :
724 5348100 : SELECT CASE (my_nspins)
725 : CASE (1)
726 : ! Give rho to 1/3
727 5348100 : IF (needs%rho_1_3) THEN
728 6735000 : DO ia = 1, na
729 6735000 : rho_set%rho_1_3(ia, ir, 1) = MAX(rho(ia, ir, 1), 0.0_dp)**f13
730 : END DO
731 132200 : rho_set%owns%rho_1_3 = .TRUE.
732 132200 : rho_set%has%rho_1_3 = .TRUE.
733 : END IF
734 : ! Give the density
735 5348100 : IF (needs%rho) THEN
736 272803500 : DO ia = 1, na
737 272803500 : rho_set%rho(ia, ir, 1) = rho(ia, ir, 1)
738 : END DO
739 5348100 : rho_set%owns%rho = .TRUE.
740 5348100 : rho_set%has%rho = .TRUE.
741 : END IF
742 : ! Give the norm of the gradient of the density
743 5348100 : IF (needs%norm_drho) THEN
744 176481900 : DO ia = 1, na
745 176481900 : rho_set%norm_drho(ia, ir, 1) = drho(4, ia, ir, 1)
746 : END DO
747 3456900 : rho_set%owns%norm_drho = .TRUE.
748 3456900 : rho_set%has%norm_drho = .TRUE.
749 : END IF
750 : ! Give the three components of the gradient of the density
751 5348100 : IF (needs%drho) THEN
752 13856400 : DO idir = 1, 3
753 533622600 : DO ia = 1, na
754 530158500 : rho_set%drho(idir)%array(ia, ir, 1) = drho(idir, ia, ir, 1)
755 : END DO
756 : END DO
757 3464100 : rho_set%owns%drho = .TRUE.
758 3464100 : rho_set%has%drho = .TRUE.
759 : END IF
760 : CASE (2)
761 : ! Give the total density
762 828980 : IF (needs%rho) THEN
763 : ! this should never be the case unless you use LDA functionals with LSD
764 0 : IF (.NOT. tddft_split) THEN
765 0 : DO ia = 1, na
766 0 : rho_set%rho(ia, ir, 1) = rho(ia, ir, 1) + rho(ia, ir, 2)
767 : END DO
768 : ELSE
769 0 : DO ia = 1, na
770 0 : rho_set%rho(ia, ir, 1) = rho(ia, ir, 1)
771 : END DO
772 : END IF
773 0 : rho_set%owns%rho = .TRUE.
774 0 : rho_set%has%rho = .TRUE.
775 : END IF
776 : ! Give the total density to 1/3
777 828980 : IF (needs%rho_1_3) THEN
778 0 : IF (.NOT. tddft_split) THEN
779 0 : DO ia = 1, na
780 0 : rho_set%rho_1_3(ia, ir, 1) = MAX(rho(ia, ir, 1) + rho(ia, ir, 2), 0.0_dp)**f13
781 : END DO
782 : ELSE
783 0 : DO ia = 1, na
784 0 : rho_set%rho_1_3(ia, ir, 1) = MAX(rho(ia, ir, 1), 0.0_dp)**f13
785 : END DO
786 : END IF
787 0 : rho_set%owns%rho_1_3 = .TRUE.
788 0 : rho_set%has%rho_1_3 = .TRUE.
789 : END IF
790 : ! Give the spin densities to 1/3
791 828980 : IF (needs%rho_spin_1_3) THEN
792 75480 : IF (.NOT. tddft_split) THEN
793 3837960 : DO ia = 1, na
794 3762480 : rho_set%rhoa_1_3(ia, ir, 1) = MAX(rho(ia, ir, 1), 0.0_dp)**f13
795 3837960 : rho_set%rhob_1_3(ia, ir, 1) = MAX(rho(ia, ir, 2), 0.0_dp)**f13
796 : END DO
797 : ELSE
798 0 : DO ia = 1, na
799 0 : rho_set%rhoa_1_3(ia, ir, 1) = MAX(0.5_dp*rho(ia, ir, 1), 0.0_dp)**f13
800 0 : rho_set%rhob_1_3(ia, ir, 1) = rho_set%rhoa_1_3(ia, ir, 1)
801 : END DO
802 : END IF
803 75480 : rho_set%owns%rho_spin_1_3 = .TRUE.
804 75480 : rho_set%has%rho_spin_1_3 = .TRUE.
805 : END IF
806 : ! Give the spin densities rhoa and rhob
807 828980 : IF (needs%rho_spin) THEN
808 828980 : IF (.NOT. tddft_split) THEN
809 37645860 : DO ia = 1, na
810 36907480 : rho_set%rhoa(ia, ir, 1) = rho(ia, ir, 1)
811 37645860 : rho_set%rhob(ia, ir, 1) = rho(ia, ir, 2)
812 : END DO
813 : ELSE
814 4620600 : DO ia = 1, na
815 4530000 : rho_set%rhoa(ia, ir, 1) = 0.5_dp*rho(ia, ir, 1)
816 4620600 : rho_set%rhob(ia, ir, 1) = rho_set%rhoa(ia, ir, 1)
817 : END DO
818 : END IF
819 828980 : rho_set%owns%rho_spin = .TRUE.
820 828980 : rho_set%has%rho_spin = .TRUE.
821 : END IF
822 : ! Give the norm of the gradient of the total density
823 828980 : IF (needs%norm_drho) THEN
824 483480 : IF (.NOT. tddft_split) THEN
825 21249360 : DO ia = 1, na
826 : rho_set%norm_drho(ia, ir, 1) = SQRT( &
827 : (drho(1, ia, ir, 1) + drho(1, ia, ir, 2))**2 + &
828 : (drho(2, ia, ir, 1) + drho(2, ia, ir, 2))**2 + &
829 21249360 : (drho(3, ia, ir, 1) + drho(3, ia, ir, 2))**2)
830 : END DO
831 : ELSE
832 3396600 : DO ia = 1, na
833 3396600 : rho_set%norm_drho(ia, ir, 1) = drho(4, ia, ir, 1)
834 : END DO
835 : END IF
836 483480 : rho_set%owns%norm_drho = .TRUE.
837 483480 : rho_set%has%norm_drho = .TRUE.
838 : END IF
839 : ! Give the norm of the gradient of rhoa and of rhob separatedly
840 828980 : IF (needs%norm_drho_spin) THEN
841 497880 : IF (.NOT. tddft_split) THEN
842 21983760 : DO ia = 1, na
843 21552480 : rho_set%norm_drhoa(ia, ir, 1) = drho(4, ia, ir, 1)
844 21983760 : rho_set%norm_drhob(ia, ir, 1) = drho(4, ia, ir, 2)
845 : END DO
846 : ELSE
847 3396600 : DO ia = 1, na
848 3330000 : rho_set%norm_drhoa(ia, ir, 1) = 0.5_dp*drho(4, ia, ir, 1)
849 3396600 : rho_set%norm_drhob(ia, ir, 1) = rho_set%norm_drhoa(ia, ir, 1)
850 : END DO
851 : END IF
852 497880 : rho_set%owns%norm_drho_spin = .TRUE.
853 497880 : rho_set%has%norm_drho_spin = .TRUE.
854 : END IF
855 : ! Give the components of the gradient for the total rho
856 828980 : IF (needs%drho) THEN
857 0 : IF (.NOT. tddft_split) THEN
858 0 : DO idir = 1, 3
859 0 : DO ia = 1, na
860 0 : rho_set%drho(idir)%array(ia, ir, 1) = drho(idir, ia, ir, 1) + drho(idir, ia, ir, 2)
861 : END DO
862 : END DO
863 : ELSE
864 0 : DO idir = 1, 3
865 0 : DO ia = 1, na
866 0 : rho_set%drho(idir)%array(ia, ir, 1) = drho(idir, ia, ir, 1)
867 : END DO
868 : END DO
869 : END IF
870 0 : rho_set%owns%drho = .TRUE.
871 0 : rho_set%has%drho = .TRUE.
872 : END IF
873 : ! Give the components of the gradient for rhoa and rhob
874 7006060 : IF (needs%drho_spin) THEN
875 499380 : IF (.NOT. tddft_split) THEN
876 1731120 : DO idir = 1, 3
877 66613560 : DO ia = 1, na
878 64882440 : rho_set%drhoa(idir)%array(ia, ir, 1) = drho(idir, ia, ir, 1)
879 66180780 : rho_set%drhob(idir)%array(ia, ir, 1) = drho(idir, ia, ir, 2)
880 : END DO
881 : END DO
882 : ELSE
883 266400 : DO idir = 1, 3
884 10256400 : DO ia = 1, na
885 9990000 : rho_set%drhoa(idir)%array(ia, ir, 1) = 0.5_dp*drho(idir, ia, ir, 1)
886 10189800 : rho_set%drhob(idir)%array(ia, ir, 1) = rho_set%drhoa(idir)%array(ia, ir, 1)
887 : END DO
888 : END DO
889 : END IF
890 499380 : rho_set%owns%drho_spin = .TRUE.
891 499380 : rho_set%has%drho_spin = .TRUE.
892 : END IF
893 : !
894 : END SELECT
895 :
896 : ! tau part
897 6177080 : IF (needs%tau .OR. needs%tau_spin) THEN
898 6177080 : CPASSERT(SIZE(tau, 3) == my_nspins)
899 : END IF
900 6177080 : IF (needs%tau) THEN
901 90000 : IF (my_nspins == 2) THEN
902 0 : DO ia = 1, na
903 0 : rho_set%tau(ia, ir, 1) = tau(ia, ir, 1) + tau(ia, ir, 2)
904 : END DO
905 0 : rho_set%owns%tau = .TRUE.
906 0 : rho_set%has%tau = .TRUE.
907 : ELSE
908 4647600 : DO ia = 1, na
909 4647600 : rho_set%tau(ia, ir, 1) = tau(ia, ir, 1)
910 : END DO
911 90000 : rho_set%owns%tau = .TRUE.
912 90000 : rho_set%has%tau = .TRUE.
913 : END IF
914 : END IF
915 6177080 : IF (needs%tau_spin) THEN
916 40800 : DO ia = 1, na
917 40000 : rho_set%tau_a(ia, ir, 1) = tau(ia, ir, 1)
918 40800 : rho_set%tau_b(ia, ir, 1) = tau(ia, ir, 2)
919 : END DO
920 800 : rho_set%owns%tau_spin = .TRUE.
921 800 : rho_set%has%tau_spin = .TRUE.
922 : END IF
923 :
924 6177080 : END SUBROUTINE fill_rho_set
925 :
926 : END MODULE xc_atom
|
