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 73096 : 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 73096 : REAL(dp), DIMENSION(:, :, :), POINTER :: deriv_data
96 : REAL(KIND=dp) :: drho_cutoff
97 : TYPE(xc_derivative_type), POINTER :: deriv_att
98 :
99 73096 : CALL timeset(routineN, handle)
100 73096 : my_only_energy = .FALSE.
101 73096 : IF (PRESENT(energy_only)) my_only_energy = energy_only
102 :
103 73096 : IF (PRESENT(adiabatic_rescale_factor)) THEN
104 73096 : 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 73096 : 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 73096 : deriv_order=deriv_order)
118 :
119 73096 : CALL xc_rho_set_get(rho_set, drho_cutoff=drho_cutoff)
120 :
121 73096 : NULLIFY (deriv_data)
122 :
123 : ! EXC energy
124 73096 : deriv_att => xc_dset_get_derivative(deriv_set, [INTEGER::])
125 73096 : exc = 0.0_dp
126 73096 : IF (ASSOCIATED(deriv_att)) THEN
127 73024 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
128 4065304 : DO ir = 1, nr
129 203847784 : DO ia = 1, na
130 203774760 : exc = exc + deriv_data(ia, ir, 1)*w(ia, ir)
131 : END DO
132 : END DO
133 73024 : NULLIFY (deriv_data)
134 : END IF
135 : ! Calculate the potential only if needed
136 73096 : IF (.NOT. my_only_energy) THEN
137 : ! Derivative with respect to the density
138 70044 : IF (lsd) THEN
139 9216 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_rhoa])
140 9216 : IF (ASSOCIATED(deriv_att)) THEN
141 9212 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
142 59530744 : vxc(:, :, 1) = deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
143 9212 : NULLIFY (deriv_data)
144 : END IF
145 9216 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_rhob])
146 9216 : IF (ASSOCIATED(deriv_att)) THEN
147 9212 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
148 59530744 : vxc(:, :, 2) = deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
149 9212 : NULLIFY (deriv_data)
150 : END IF
151 9216 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_rho])
152 9216 : 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 60828 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_rho])
160 60828 : IF (ASSOCIATED(deriv_att)) THEN
161 60760 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
162 331206320 : vxc(:, :, 1) = deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
163 60760 : NULLIFY (deriv_data)
164 : END IF
165 : END IF ! lsd
166 :
167 : ! Derivatives with respect to the gradient
168 70044 : IF (lsd) THEN
169 9216 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa])
170 9216 : IF (ASSOCIATED(deriv_att)) THEN
171 5702 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
172 306082 : DO ir = 1, nr
173 15313562 : DO ia = 1, na
174 60330300 : DO idir = 1, 3
175 60029920 : IF (rho_set%norm_drhoa(ia, ir, 1) > drho_cutoff) THEN
176 : vxg(idir, ia, ir, 1) = rho_set%drhoa(idir)%array(ia, ir, 1)* &
177 : deriv_data(ia, ir, 1)*w(ia, ir)/ &
178 38538702 : rho_set%norm_drhoa(ia, ir, 1)*my_adiabatic_rescale_factor
179 : ELSE
180 6483738 : vxg(idir, ia, ir, 1) = 0.0_dp
181 : END IF
182 : END DO
183 : END DO
184 : END DO
185 5702 : NULLIFY (deriv_data)
186 : END IF
187 9216 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob])
188 9216 : IF (ASSOCIATED(deriv_att)) THEN
189 5702 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
190 306082 : DO ir = 1, nr
191 15313562 : DO ia = 1, na
192 60330300 : DO idir = 1, 3
193 60029920 : IF (rho_set%norm_drhob(ia, ir, 1) > drho_cutoff) THEN
194 : vxg(idir, ia, ir, 2) = rho_set%drhob(idir)%array(ia, ir, 1)* &
195 : deriv_data(ia, ir, 1)*w(ia, ir)/ &
196 37883328 : rho_set%norm_drhob(ia, ir, 1)*my_adiabatic_rescale_factor
197 : ELSE
198 7139112 : vxg(idir, ia, ir, 2) = 0.0_dp
199 : END IF
200 : END DO
201 : END DO
202 : END DO
203 5702 : NULLIFY (deriv_data)
204 : END IF
205 : ! Cross Terms
206 9216 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_norm_drho])
207 9216 : IF (ASSOCIATED(deriv_att)) THEN
208 5414 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
209 291394 : DO ir = 1, nr
210 14578874 : DO ia = 1, na
211 57435900 : DO idir = 1, 3
212 57149920 : IF (rho_set%norm_drho(ia, ir, 1) > drho_cutoff) THEN
213 : vxg(idir, ia, ir, 1:2) = &
214 : vxg(idir, ia, ir, 1:2) + ( &
215 : rho_set%drhoa(idir)%array(ia, ir, 1) + &
216 : rho_set%drhob(idir)%array(ia, ir, 1))* &
217 : deriv_data(ia, ir, 1)*w(ia, ir)/rho_set%norm_drho(ia, ir, 1)* &
218 110301804 : my_adiabatic_rescale_factor
219 : END IF
220 : END DO
221 : END DO
222 : END DO
223 5414 : NULLIFY (deriv_data)
224 : END IF
225 : ELSE
226 60828 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_norm_drho])
227 60828 : IF (ASSOCIATED(deriv_att)) THEN
228 37858 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
229 1950358 : DO ir = 1, nr
230 97755358 : DO ia = 1, na
231 97717500 : IF (rho_set%norm_drho(ia, ir, 1) > drho_cutoff) THEN
232 331002208 : DO idir = 1, 3
233 : vxg(idir, ia, ir, 1) = rho_set%drho(idir)%array(ia, ir, 1)* &
234 : deriv_data(ia, ir, 1)*w(ia, ir)/ &
235 331002208 : rho_set%norm_drho(ia, ir, 1)*my_adiabatic_rescale_factor
236 : END DO
237 : ELSE
238 52217792 : vxg(1:3, ia, ir, 1) = 0.0_dp
239 : END IF
240 : END DO
241 : END DO
242 37858 : NULLIFY (deriv_data)
243 : END IF
244 : END IF ! lsd
245 : ! Derivative with respect to tau
246 70044 : IF (lsd) THEN
247 9216 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_tau_a])
248 9216 : IF (ASSOCIATED(deriv_att)) THEN
249 16 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
250 81632 : vtau(:, :, 1) = deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
251 16 : NULLIFY (deriv_data)
252 : END IF
253 9216 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_tau_b])
254 9216 : IF (ASSOCIATED(deriv_att)) THEN
255 16 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
256 81632 : vtau(:, :, 2) = deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
257 16 : NULLIFY (deriv_data)
258 : END IF
259 9216 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_tau])
260 9216 : IF (ASSOCIATED(deriv_att)) THEN
261 0 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
262 0 : vtau(:, :, 1) = vtau(:, :, 1) + deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
263 0 : vtau(:, :, 2) = vtau(:, :, 2) + deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
264 0 : NULLIFY (deriv_data)
265 : END IF
266 : ELSE
267 60828 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_tau])
268 60828 : IF (ASSOCIATED(deriv_att)) THEN
269 1682 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
270 9221164 : vtau(:, :, 1) = deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
271 1682 : NULLIFY (deriv_data)
272 : END IF
273 : END IF ! lsd
274 : END IF ! only_energy
275 :
276 73096 : CALL timestop(handle)
277 :
278 73096 : END SUBROUTINE vxc_of_r_new
279 :
280 : ! **************************************************************************************************
281 : !> \brief Specific EPR version of vxc_of_r_new
282 : !> \param xc_fun_section ...
283 : !> \param rho_set ...
284 : !> \param deriv_set ...
285 : !> \param needs ...
286 : !> \param w ...
287 : !> \param lsd ...
288 : !> \param na ...
289 : !> \param nr ...
290 : !> \param exc ...
291 : !> \param vxc ...
292 : !> \param vxg ...
293 : !> \param vtau ...
294 : ! **************************************************************************************************
295 30 : SUBROUTINE vxc_of_r_epr(xc_fun_section, rho_set, deriv_set, needs, w, &
296 : lsd, na, nr, exc, vxc, vxg, vtau)
297 :
298 : TYPE(section_vals_type), POINTER :: xc_fun_section
299 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
300 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
301 : TYPE(xc_rho_cflags_type), INTENT(IN) :: needs
302 : REAL(dp), DIMENSION(:, :), INTENT(IN) :: w
303 : LOGICAL, INTENT(IN) :: lsd
304 : INTEGER, INTENT(in) :: na, nr
305 : REAL(dp) :: exc
306 : REAL(dp), DIMENSION(:, :, :), POINTER :: vxc
307 : REAL(dp), DIMENSION(:, :, :, :), POINTER :: vxg
308 : REAL(dp), DIMENSION(:, :, :), POINTER :: vtau
309 :
310 : CHARACTER(LEN=*), PARAMETER :: routineN = 'vxc_of_r_epr'
311 :
312 : INTEGER :: handle, ia, idir, ir, my_deriv_order
313 : LOGICAL :: gradient_f
314 : REAL(dp) :: my_adiabatic_rescale_factor
315 30 : REAL(dp), DIMENSION(:, :, :), POINTER :: deriv_data
316 : REAL(KIND=dp) :: drho_cutoff
317 : TYPE(xc_derivative_type), POINTER :: deriv_att
318 :
319 30 : CALL timeset(routineN, handle)
320 :
321 : MARK_USED(vxc)
322 : MARK_USED(vtau)
323 :
324 30 : my_adiabatic_rescale_factor = 1.0_dp
325 30 : my_deriv_order = 2
326 :
327 : gradient_f = (needs%drho_spin .OR. needs%norm_drho_spin .OR. &
328 30 : needs%drho .OR. needs%norm_drho)
329 :
330 : ! Calculate the derivatives
331 : CALL xc_functionals_eval(xc_fun_section, &
332 : lsd=lsd, &
333 : rho_set=rho_set, &
334 : deriv_set=deriv_set, &
335 30 : deriv_order=my_deriv_order)
336 :
337 30 : CALL xc_rho_set_get(rho_set, drho_cutoff=drho_cutoff)
338 :
339 30 : NULLIFY (deriv_data)
340 :
341 : ! nabla v_xc (using the vxg arrays)
342 : ! there's no point doing this when lsd = false
343 30 : IF (lsd) THEN
344 30 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa])
345 30 : IF (ASSOCIATED(deriv_att)) THEN
346 30 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
347 1530 : DO ir = 1, nr
348 76530 : DO ia = 1, na
349 301500 : DO idir = 1, 3
350 : vxg(idir, ia, ir, 1) = rho_set%drhoa(idir)%array(ia, ir, 1)* &
351 300000 : deriv_data(ia, ir, 1)
352 : END DO !idir
353 : END DO !ia
354 : END DO !ir
355 30 : NULLIFY (deriv_data)
356 : END IF
357 30 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob])
358 30 : IF (ASSOCIATED(deriv_att)) THEN
359 30 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
360 1530 : DO ir = 1, nr
361 76530 : DO ia = 1, na
362 301500 : DO idir = 1, 3
363 : vxg(idir, ia, ir, 2) = rho_set%drhob(idir)%array(ia, ir, 1)* &
364 300000 : deriv_data(ia, ir, 1)
365 : END DO !idir
366 : END DO !ia
367 : END DO !ir
368 30 : NULLIFY (deriv_data)
369 : END IF
370 : END IF
371 : ! EXC energy ! is that needed for epr?
372 30 : deriv_att => xc_dset_get_derivative(deriv_set, [INTEGER::])
373 30 : exc = 0.0_dp
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 76500 : exc = exc + deriv_data(ia, ir, 1)*w(ia, ir)
379 : END DO
380 : END DO
381 30 : NULLIFY (deriv_data)
382 : END IF
383 :
384 30 : CALL timestop(handle)
385 :
386 30 : END SUBROUTINE vxc_of_r_epr
387 :
388 : ! **************************************************************************************************
389 : !> \brief ...
390 : !> \param rho_set ...
391 : !> \param rho1_set ...
392 : !> \param xc_section ...
393 : !> \param deriv_set ...
394 : !> \param w ...
395 : !> \param vxc ...
396 : !> \param vxg ...
397 : !> \param vtau ...
398 : !> \param do_triplet ...
399 : !> \param do_sf ...
400 : ! **************************************************************************************************
401 19724 : SUBROUTINE xc_2nd_deriv_of_r(rho_set, rho1_set, xc_section, &
402 19724 : deriv_set, w, vxc, vxg, vtau, do_triplet, do_sf)
403 :
404 : ! As input of this routine one gets rho and drho on a one dimensional grid.
405 : ! The grid is the angular grid corresponding to a given point ir on the radial grid.
406 : ! The derivatives are calculated on this one dimensional grid, the results are stored in
407 : ! vxc(1:na,ir,ispin), vxg(1:na,ir,ispin), vxg_cross(1:na,ir,ispin)
408 : ! Afterwords the arrays containing the derivatives are put to zero so that the routine
409 : ! can safely be called for the next radial point ir
410 :
411 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set, rho1_set
412 : TYPE(section_vals_type), POINTER :: xc_section
413 : TYPE(xc_derivative_set_type), INTENT(INOUT) :: deriv_set
414 : REAL(dp), DIMENSION(:, :), INTENT(IN) :: w
415 : REAL(dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: vxc
416 : REAL(dp), DIMENSION(:, :, :, :), POINTER :: vxg
417 : REAL(dp), CONTIGUOUS, DIMENSION(:, :, :), &
418 : OPTIONAL, POINTER :: vtau
419 : LOGICAL, INTENT(IN), OPTIONAL :: do_triplet, do_sf
420 :
421 : CHARACTER(LEN=*), PARAMETER :: routineN = 'xc_2nd_deriv_of_r'
422 :
423 : INTEGER :: handle, ispin, nspins
424 : LOGICAL :: lsd, my_do_sf
425 : REAL(dp) :: drho_cutoff, my_fac_triplet
426 : TYPE(cp_sll_xc_deriv_type), POINTER :: pos
427 : TYPE(pw_pool_type), POINTER :: pw_pool
428 19724 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: vxc_pw, vxc_tau_pw
429 : TYPE(section_vals_type), POINTER :: xc_fun_section
430 : TYPE(xc_derivative_type), POINTER :: deriv_att
431 :
432 19724 : CALL timeset(routineN, handle)
433 :
434 19724 : nspins = SIZE(vxc, 3)
435 19724 : lsd = (nspins == 2)
436 19724 : IF (ASSOCIATED(rho_set%rhoa)) THEN
437 1402 : lsd = .TRUE.
438 : END IF
439 19724 : my_fac_triplet = 1.0_dp
440 19724 : IF (PRESENT(do_triplet)) THEN
441 19362 : IF (do_triplet) my_fac_triplet = -1.0_dp
442 : END IF
443 :
444 19724 : my_do_sf = .FALSE.
445 19724 : IF (PRESENT(do_sf)) my_do_sf = do_sf
446 :
447 19724 : CALL xc_rho_set_get(rho_set, drho_cutoff=drho_cutoff)
448 : xc_fun_section => section_vals_get_subs_vals(xc_section, &
449 19724 : "XC_FUNCTIONAL")
450 :
451 : ! Calculate the derivatives
452 : CALL xc_functionals_eval(xc_fun_section, &
453 : lsd=lsd, &
454 : rho_set=rho_set, &
455 : deriv_set=deriv_set, &
456 19724 : deriv_order=2)
457 :
458 19724 : CALL divide_by_norm_drho(deriv_set, rho_set, lsd)
459 :
460 : ! multiply by w
461 19724 : pos => deriv_set%derivs
462 130192 : DO WHILE (cp_sll_xc_deriv_next(pos, el_att=deriv_att))
463 281823592 : deriv_att%deriv_data(:, :, 1) = w(:, :)*deriv_att%deriv_data(:, :, 1)
464 : END DO
465 :
466 19724 : NULLIFY (pw_pool)
467 79284 : ALLOCATE (vxc_pw(nspins))
468 39836 : DO ispin = 1, nspins
469 39836 : vxc_pw(ispin)%array => vxc(:, :, ispin:ispin)
470 : END DO
471 :
472 19724 : NULLIFY (vxc_tau_pw)
473 19724 : IF (PRESENT(vtau)) THEN
474 19724 : IF (ASSOCIATED(vtau)) THEN
475 0 : ALLOCATE (vxc_tau_pw(nspins))
476 0 : DO ispin = 1, nspins
477 0 : vxc_tau_pw(ispin)%array => vtau(:, :, ispin:ispin)
478 : END DO
479 : END IF
480 : END IF
481 :
482 : CALL xc_calc_2nd_deriv_analytical(vxc_pw, vxc_tau_pw, deriv_set, rho_set, rho1_set, pw_pool, &
483 : xc_section, gapw=.TRUE., vxg=vxg, &
484 19724 : tddfpt_fac=my_fac_triplet, spinflip=do_sf)
485 :
486 19724 : DEALLOCATE (vxc_pw)
487 19724 : IF (ASSOCIATED(vxc_tau_pw)) DEALLOCATE (vxc_tau_pw)
488 :
489 : ! zero the derivative data for the next call
490 19724 : pos => deriv_set%derivs
491 130192 : DO WHILE (cp_sll_xc_deriv_next(pos, el_att=deriv_att))
492 281934060 : deriv_att%deriv_data = 0.0_dp
493 : END DO
494 :
495 19724 : CALL timestop(handle)
496 :
497 39448 : END SUBROUTINE xc_2nd_deriv_of_r
498 :
499 : ! **************************************************************************************************
500 : !> \brief ...
501 : !> \param rho_set ...
502 : !> \param needs ...
503 : !> \param nspins ...
504 : !> \param bo ...
505 : ! **************************************************************************************************
506 165529 : SUBROUTINE xc_rho_set_atom_update(rho_set, needs, nspins, bo)
507 :
508 : ! This routine allocates the storage arrays for rho and drho
509 : ! In calculate_vxc_atom this is called once for each atomic_kind,
510 : ! After the loop over all the atoms of the kind and over all the points
511 : ! of the radial grid for each atom, rho_set is deallocated.
512 : ! Within the same kind, at each new point on the radial grid, the rho_set
513 : ! arrays rho and drho are overwritten.
514 :
515 : TYPE(xc_rho_set_type), INTENT(INOUT) :: rho_set
516 : TYPE(xc_rho_cflags_type), INTENT(IN) :: needs
517 : INTEGER, INTENT(IN) :: nspins
518 : INTEGER, DIMENSION(2, 3), INTENT(IN) :: bo
519 :
520 : INTEGER :: idir
521 :
522 313897 : SELECT CASE (nspins)
523 : CASE (1)
524 : ! What is this for?
525 148368 : IF (needs%rho_1_3) THEN
526 4055 : NULLIFY (rho_set%rho_1_3)
527 20275 : ALLOCATE (rho_set%rho_1_3(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
528 4055 : rho_set%owns%rho_1_3 = .TRUE.
529 4055 : rho_set%has%rho_1_3 = .FALSE.
530 : END IF
531 : ! Allocate the storage space for the density
532 148368 : IF (needs%rho) THEN
533 148368 : NULLIFY (rho_set%rho)
534 741840 : ALLOCATE (rho_set%rho(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
535 148368 : rho_set%owns%rho = .TRUE.
536 148368 : rho_set%has%rho = .FALSE.
537 : END IF
538 : ! Allocate the storage space for the norm of the gradient of the density
539 148368 : IF (needs%norm_drho) THEN
540 107186 : NULLIFY (rho_set%norm_drho)
541 535930 : ALLOCATE (rho_set%norm_drho(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
542 107186 : rho_set%owns%norm_drho = .TRUE.
543 107186 : rho_set%has%norm_drho = .FALSE.
544 : END IF
545 : ! Allocate the storage space for the three components of the gradient of the density
546 148368 : IF (needs%drho) THEN
547 342064 : DO idir = 1, 3
548 256548 : NULLIFY (rho_set%drho(idir)%array)
549 1368256 : ALLOCATE (rho_set%drho(idir)%array(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
550 : END DO
551 85516 : rho_set%owns%drho = .TRUE.
552 85516 : rho_set%has%drho = .FALSE.
553 : END IF
554 : CASE (2)
555 : ! Allocate the storage space for the total density
556 17161 : IF (needs%rho) THEN
557 : ! this should never be the case unless you use LDA functionals with LSD
558 0 : NULLIFY (rho_set%rho)
559 0 : ALLOCATE (rho_set%rho(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
560 0 : rho_set%owns%rho = .TRUE.
561 0 : rho_set%has%rho = .FALSE.
562 : END IF
563 : ! What is this for?
564 17161 : IF (needs%rho_1_3) THEN
565 0 : NULLIFY (rho_set%rho_1_3)
566 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)))
567 0 : rho_set%owns%rho_1_3 = .TRUE.
568 0 : rho_set%has%rho_1_3 = .FALSE.
569 : END IF
570 : ! What is this for?
571 17161 : IF (needs%rho_spin_1_3) THEN
572 2440 : NULLIFY (rho_set%rhoa_1_3, rho_set%rhob_1_3)
573 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)))
574 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)))
575 2440 : rho_set%owns%rho_spin_1_3 = .TRUE.
576 2440 : rho_set%has%rho_spin_1_3 = .FALSE.
577 : END IF
578 : ! Allocate the storage space for the spin densities rhoa and rhob
579 17161 : IF (needs%rho_spin) THEN
580 17161 : NULLIFY (rho_set%rhoa, rho_set%rhob)
581 85805 : ALLOCATE (rho_set%rhoa(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
582 68644 : ALLOCATE (rho_set%rhob(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
583 17161 : rho_set%owns%rho_spin = .TRUE.
584 17161 : rho_set%has%rho_spin = .FALSE.
585 : END IF
586 : ! Allocate the storage space for the norm of the gradient of the total density
587 17161 : IF (needs%norm_drho) THEN
588 10847 : NULLIFY (rho_set%norm_drho)
589 54235 : ALLOCATE (rho_set%norm_drho(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
590 10847 : rho_set%owns%norm_drho = .TRUE.
591 10847 : rho_set%has%norm_drho = .FALSE.
592 : END IF
593 : ! Allocate the storage space for the norm of the gradient of rhoa and of rhob separatedly
594 17161 : IF (needs%norm_drho_spin) THEN
595 11135 : NULLIFY (rho_set%norm_drhoa, rho_set%norm_drhob)
596 55675 : ALLOCATE (rho_set%norm_drhoa(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
597 44540 : ALLOCATE (rho_set%norm_drhob(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
598 11135 : rho_set%owns%norm_drho_spin = .TRUE.
599 11135 : rho_set%has%norm_drho_spin = .FALSE.
600 : END IF
601 : ! Allocate the storage space for the components of the gradient for the total rho
602 17161 : IF (needs%drho) THEN
603 0 : DO idir = 1, 3
604 0 : NULLIFY (rho_set%drho(idir)%array)
605 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)))
606 : END DO
607 0 : rho_set%owns%drho = .TRUE.
608 0 : rho_set%has%drho = .FALSE.
609 : END IF
610 : ! Allocate the storage space for the components of the gradient for rhoa and rhob
611 182690 : IF (needs%drho_spin) THEN
612 41024 : DO idir = 1, 3
613 30768 : NULLIFY (rho_set%drhoa(idir)%array, rho_set%drhob(idir)%array)
614 153840 : ALLOCATE (rho_set%drhoa(idir)%array(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
615 133328 : ALLOCATE (rho_set%drhob(idir)%array(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
616 : END DO
617 10256 : rho_set%owns%drho_spin = .TRUE.
618 10256 : rho_set%has%drho_spin = .FALSE.
619 : END IF
620 : !
621 : END SELECT
622 :
623 : ! tau part
624 165529 : IF (needs%tau) THEN
625 2512 : NULLIFY (rho_set%tau)
626 12560 : ALLOCATE (rho_set%tau(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
627 2512 : rho_set%owns%tau = .TRUE.
628 : END IF
629 165529 : IF (needs%tau_spin) THEN
630 34 : NULLIFY (rho_set%tau_a, rho_set%tau_b)
631 170 : ALLOCATE (rho_set%tau_a(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
632 136 : ALLOCATE (rho_set%tau_b(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
633 34 : rho_set%owns%tau_spin = .TRUE.
634 34 : rho_set%has%tau_spin = .FALSE.
635 : END IF
636 :
637 : ! Laplace part
638 165529 : IF (needs%laplace_rho) THEN
639 0 : NULLIFY (rho_set%laplace_rho)
640 0 : ALLOCATE (rho_set%laplace_rho(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
641 0 : rho_set%owns%laplace_rho = .TRUE.
642 : END IF
643 165529 : IF (needs%laplace_rho_spin) THEN
644 0 : NULLIFY (rho_set%laplace_rhoa)
645 0 : NULLIFY (rho_set%laplace_rhob)
646 0 : ALLOCATE (rho_set%laplace_rhoa(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
647 0 : ALLOCATE (rho_set%laplace_rhob(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
648 0 : rho_set%owns%laplace_rho_spin = .TRUE.
649 0 : rho_set%has%laplace_rho_spin = .TRUE.
650 : END IF
651 :
652 165529 : END SUBROUTINE xc_rho_set_atom_update
653 :
654 : ! **************************************************************************************************
655 : !> \brief ...
656 : !> \param rho_set ...
657 : !> \param lsd ...
658 : !> \param nspins ...
659 : !> \param needs ...
660 : !> \param rho ...
661 : !> \param drho ...
662 : !> \param tau ...
663 : !> \param na ...
664 : !> \param ir ...
665 : ! **************************************************************************************************
666 5859980 : SUBROUTINE fill_rho_set(rho_set, lsd, nspins, needs, rho, drho, tau, na, ir)
667 :
668 : TYPE(xc_rho_set_type), INTENT(INOUT) :: rho_set
669 : LOGICAL, INTENT(IN) :: lsd
670 : INTEGER, INTENT(IN) :: nspins
671 : TYPE(xc_rho_cflags_type), INTENT(IN) :: needs
672 : REAL(dp), DIMENSION(:, :, :), INTENT(IN) :: rho
673 : REAL(dp), DIMENSION(:, :, :, :), INTENT(IN) :: drho
674 : REAL(dp), DIMENSION(:, :, :), INTENT(IN) :: tau
675 : INTEGER, INTENT(IN) :: na, ir
676 :
677 : REAL(KIND=dp), PARAMETER :: f13 = (1.0_dp/3.0_dp)
678 :
679 : INTEGER :: ia, idir, my_nspins
680 : LOGICAL :: gradient_f, tddft_split
681 :
682 5859980 : my_nspins = nspins
683 5859980 : tddft_split = .FALSE.
684 5859980 : IF (lsd .AND. nspins == 1) THEN
685 90600 : my_nspins = 2
686 90600 : tddft_split = .TRUE.
687 : END IF
688 :
689 : ! some checks
690 5859980 : IF (lsd) THEN
691 : ELSE
692 5079300 : CPASSERT(SIZE(rho, 3) == 1)
693 : END IF
694 5079300 : SELECT CASE (my_nspins)
695 : CASE (1)
696 5079300 : CPASSERT(.NOT. needs%rho_spin)
697 5079300 : CPASSERT(.NOT. needs%drho_spin)
698 5079300 : CPASSERT(.NOT. needs%norm_drho_spin)
699 5079300 : CPASSERT(.NOT. needs%rho_spin_1_3)
700 : CASE (2)
701 : CASE default
702 5859980 : CPABORT("Unsupported number of spins")
703 : END SELECT
704 :
705 : gradient_f = (needs%drho_spin .OR. needs%norm_drho_spin .OR. &
706 5859980 : needs%drho .OR. needs%norm_drho)
707 :
708 5079300 : SELECT CASE (my_nspins)
709 : CASE (1)
710 : ! Give rho to 1/3
711 5079300 : IF (needs%rho_1_3) THEN
712 6765600 : DO ia = 1, na
713 6765600 : rho_set%rho_1_3(ia, ir, 1) = MAX(rho(ia, ir, 1), 0.0_dp)**f13
714 : END DO
715 132800 : rho_set%owns%rho_1_3 = .TRUE.
716 132800 : rho_set%has%rho_1_3 = .TRUE.
717 : END IF
718 : ! Give the density
719 5079300 : IF (needs%rho) THEN
720 259224300 : DO ia = 1, na
721 259224300 : rho_set%rho(ia, ir, 1) = rho(ia, ir, 1)
722 : END DO
723 5079300 : rho_set%owns%rho = .TRUE.
724 5079300 : rho_set%has%rho = .TRUE.
725 : END IF
726 : ! Give the norm of the gradient of the density
727 5079300 : IF (needs%norm_drho) THEN
728 164971200 : DO ia = 1, na
729 164971200 : rho_set%norm_drho(ia, ir, 1) = drho(4, ia, ir, 1)
730 : END DO
731 3231200 : rho_set%owns%norm_drho = .TRUE.
732 3231200 : rho_set%has%norm_drho = .TRUE.
733 : END IF
734 : ! Give the three components of the gradient of the density
735 5079300 : IF (needs%drho) THEN
736 12924800 : DO idir = 1, 3
737 498144800 : DO ia = 1, na
738 494913600 : rho_set%drho(idir)%array(ia, ir, 1) = drho(idir, ia, ir, 1)
739 : END DO
740 : END DO
741 3231200 : rho_set%owns%drho = .TRUE.
742 3231200 : rho_set%has%drho = .TRUE.
743 : END IF
744 : CASE (2)
745 : ! Give the total density
746 780680 : IF (needs%rho) THEN
747 : ! this should never be the case unless you use LDA functionals with LSD
748 0 : IF (.NOT. tddft_split) THEN
749 0 : DO ia = 1, na
750 0 : rho_set%rho(ia, ir, 1) = rho(ia, ir, 1) + rho(ia, ir, 2)
751 : END DO
752 : ELSE
753 0 : DO ia = 1, na
754 0 : rho_set%rho(ia, ir, 1) = rho(ia, ir, 1)
755 : END DO
756 : END IF
757 0 : rho_set%owns%rho = .TRUE.
758 0 : rho_set%has%rho = .TRUE.
759 : END IF
760 : ! Give the total density to 1/3
761 780680 : IF (needs%rho_1_3) THEN
762 0 : IF (.NOT. tddft_split) THEN
763 0 : DO ia = 1, na
764 0 : rho_set%rho_1_3(ia, ir, 1) = MAX(rho(ia, ir, 1) + rho(ia, ir, 2), 0.0_dp)**f13
765 : END DO
766 : ELSE
767 0 : DO ia = 1, na
768 0 : rho_set%rho_1_3(ia, ir, 1) = MAX(rho(ia, ir, 1), 0.0_dp)**f13
769 : END DO
770 : END IF
771 0 : rho_set%owns%rho_1_3 = .TRUE.
772 0 : rho_set%has%rho_1_3 = .TRUE.
773 : END IF
774 : ! Give the spin densities to 1/3
775 780680 : IF (needs%rho_spin_1_3) THEN
776 75480 : IF (.NOT. tddft_split) THEN
777 3837960 : DO ia = 1, na
778 3762480 : rho_set%rhoa_1_3(ia, ir, 1) = MAX(rho(ia, ir, 1), 0.0_dp)**f13
779 3837960 : rho_set%rhob_1_3(ia, ir, 1) = MAX(rho(ia, ir, 2), 0.0_dp)**f13
780 : END DO
781 : ELSE
782 0 : DO ia = 1, na
783 0 : rho_set%rhoa_1_3(ia, ir, 1) = MAX(0.5_dp*rho(ia, ir, 1), 0.0_dp)**f13
784 0 : rho_set%rhob_1_3(ia, ir, 1) = rho_set%rhoa_1_3(ia, ir, 1)
785 : END DO
786 : END IF
787 75480 : rho_set%owns%rho_spin_1_3 = .TRUE.
788 75480 : rho_set%has%rho_spin_1_3 = .TRUE.
789 : END IF
790 : ! Give the spin densities rhoa and rhob
791 780680 : IF (needs%rho_spin) THEN
792 780680 : IF (.NOT. tddft_split) THEN
793 35182560 : DO ia = 1, na
794 34492480 : rho_set%rhoa(ia, ir, 1) = rho(ia, ir, 1)
795 35182560 : rho_set%rhob(ia, ir, 1) = rho(ia, ir, 2)
796 : END DO
797 : ELSE
798 4620600 : DO ia = 1, na
799 4530000 : rho_set%rhoa(ia, ir, 1) = 0.5_dp*rho(ia, ir, 1)
800 4620600 : rho_set%rhob(ia, ir, 1) = rho_set%rhoa(ia, ir, 1)
801 : END DO
802 : END IF
803 780680 : rho_set%owns%rho_spin = .TRUE.
804 780680 : rho_set%has%rho_spin = .TRUE.
805 : END IF
806 : ! Give the norm of the gradient of the total density
807 780680 : IF (needs%norm_drho) THEN
808 434880 : IF (.NOT. tddft_split) THEN
809 18770760 : DO ia = 1, na
810 : rho_set%norm_drho(ia, ir, 1) = SQRT( &
811 : (drho(1, ia, ir, 1) + drho(1, ia, ir, 2))**2 + &
812 : (drho(2, ia, ir, 1) + drho(2, ia, ir, 2))**2 + &
813 18770760 : (drho(3, ia, ir, 1) + drho(3, ia, ir, 2))**2)
814 : END DO
815 : ELSE
816 3396600 : DO ia = 1, na
817 3396600 : rho_set%norm_drho(ia, ir, 1) = drho(4, ia, ir, 1)
818 : END DO
819 : END IF
820 434880 : rho_set%owns%norm_drho = .TRUE.
821 434880 : rho_set%has%norm_drho = .TRUE.
822 : END IF
823 : ! Give the norm of the gradient of rhoa and of rhob separatedly
824 780680 : IF (needs%norm_drho_spin) THEN
825 449280 : IF (.NOT. tddft_split) THEN
826 19505160 : DO ia = 1, na
827 19122480 : rho_set%norm_drhoa(ia, ir, 1) = drho(4, ia, ir, 1)
828 19505160 : rho_set%norm_drhob(ia, ir, 1) = drho(4, ia, ir, 2)
829 : END DO
830 : ELSE
831 3396600 : DO ia = 1, na
832 3330000 : rho_set%norm_drhoa(ia, ir, 1) = 0.5_dp*drho(4, ia, ir, 1)
833 3396600 : rho_set%norm_drhob(ia, ir, 1) = rho_set%norm_drhoa(ia, ir, 1)
834 : END DO
835 : END IF
836 449280 : rho_set%owns%norm_drho_spin = .TRUE.
837 449280 : rho_set%has%norm_drho_spin = .TRUE.
838 : END IF
839 : ! Give the components of the gradient for the total rho
840 780680 : IF (needs%drho) THEN
841 0 : IF (.NOT. tddft_split) THEN
842 0 : DO idir = 1, 3
843 0 : DO ia = 1, na
844 0 : rho_set%drho(idir)%array(ia, ir, 1) = drho(idir, ia, ir, 1) + drho(idir, ia, ir, 2)
845 : END DO
846 : END DO
847 : ELSE
848 0 : DO idir = 1, 3
849 0 : DO ia = 1, na
850 0 : rho_set%drho(idir)%array(ia, ir, 1) = drho(idir, ia, ir, 1)
851 : END DO
852 : END DO
853 : END IF
854 0 : rho_set%owns%drho = .TRUE.
855 0 : rho_set%has%drho = .TRUE.
856 : END IF
857 : ! Give the components of the gradient for rhoa and rhob
858 6640660 : IF (needs%drho_spin) THEN
859 450780 : IF (.NOT. tddft_split) THEN
860 1536720 : DO idir = 1, 3
861 59129160 : DO ia = 1, na
862 57592440 : rho_set%drhoa(idir)%array(ia, ir, 1) = drho(idir, ia, ir, 1)
863 58744980 : rho_set%drhob(idir)%array(ia, ir, 1) = drho(idir, ia, ir, 2)
864 : END DO
865 : END DO
866 : ELSE
867 266400 : DO idir = 1, 3
868 10256400 : DO ia = 1, na
869 9990000 : rho_set%drhoa(idir)%array(ia, ir, 1) = 0.5_dp*drho(idir, ia, ir, 1)
870 10189800 : rho_set%drhob(idir)%array(ia, ir, 1) = rho_set%drhoa(idir)%array(ia, ir, 1)
871 : END DO
872 : END DO
873 : END IF
874 450780 : rho_set%owns%drho_spin = .TRUE.
875 450780 : rho_set%has%drho_spin = .TRUE.
876 : END IF
877 : !
878 : END SELECT
879 :
880 : ! tau part
881 5859980 : IF (needs%tau .OR. needs%tau_spin) THEN
882 5859980 : CPASSERT(SIZE(tau, 3) == my_nspins)
883 : END IF
884 5859980 : IF (needs%tau) THEN
885 86700 : IF (my_nspins == 2) THEN
886 0 : DO ia = 1, na
887 0 : rho_set%tau(ia, ir, 1) = tau(ia, ir, 1) + tau(ia, ir, 2)
888 : END DO
889 0 : rho_set%owns%tau = .TRUE.
890 0 : rho_set%has%tau = .TRUE.
891 : ELSE
892 4608900 : DO ia = 1, na
893 4608900 : rho_set%tau(ia, ir, 1) = tau(ia, ir, 1)
894 : END DO
895 86700 : rho_set%owns%tau = .TRUE.
896 86700 : rho_set%has%tau = .TRUE.
897 : END IF
898 : END IF
899 5859980 : IF (needs%tau_spin) THEN
900 40800 : DO ia = 1, na
901 40000 : rho_set%tau_a(ia, ir, 1) = tau(ia, ir, 1)
902 40800 : rho_set%tau_b(ia, ir, 1) = tau(ia, ir, 2)
903 : END DO
904 800 : rho_set%owns%tau_spin = .TRUE.
905 800 : rho_set%has%tau_spin = .TRUE.
906 : END IF
907 :
908 5859980 : END SUBROUTINE fill_rho_set
909 :
910 : END MODULE xc_atom
|
