Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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, 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 epr_xc ...
60 : !> \param adiabatic_rescale_factor ...
61 : ! **************************************************************************************************
62 93532 : SUBROUTINE vxc_of_r_new(xc_fun_section, rho_set, deriv_set, deriv_order, needs, w, &
63 : lsd, na, nr, exc, vxc, vxg, vtau, &
64 : energy_only, epr_xc, adiabatic_rescale_factor)
65 :
66 : ! This routine updates rho_set by giving to it the rho and drho that are needed.
67 : ! Since for the local densities rho1_h and rho1_s local grids are used it is not possible
68 : ! to call xc_rho_set_update.
69 : ! As input of this routine one gets rho and drho on a one dimensional grid.
70 : ! The grid is the angular grid corresponding to a given point ir_pnt on the radial grid.
71 : ! The derivatives are calculated on this one dimensional grid, the results are stored in
72 : ! exc, vxc(1:na,ir_pnt,ispin), vxg(1:na,ir_pnt,ispin), vxg_cross(1:na,ir_pnt,ispin)
73 : ! Afterwords the arrays containing the derivatives are put to zero so that the routine
74 : ! can safely be called for the next radial point ir_pnt
75 :
76 : TYPE(section_vals_type), POINTER :: xc_fun_section
77 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
78 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
79 : INTEGER, INTENT(in) :: deriv_order
80 : TYPE(xc_rho_cflags_type), INTENT(in) :: needs
81 : REAL(dp), DIMENSION(:, :), POINTER :: w
82 : LOGICAL, INTENT(IN) :: lsd
83 : INTEGER, INTENT(in) :: na, nr
84 : REAL(dp) :: exc
85 : REAL(dp), DIMENSION(:, :, :), POINTER :: vxc
86 : REAL(dp), DIMENSION(:, :, :, :), POINTER :: vxg
87 : REAL(dp), DIMENSION(:, :, :), POINTER :: vtau
88 : LOGICAL, INTENT(IN), OPTIONAL :: energy_only, epr_xc
89 : REAL(dp), INTENT(IN), OPTIONAL :: adiabatic_rescale_factor
90 :
91 : CHARACTER(LEN=*), PARAMETER :: routineN = 'vxc_of_r_new'
92 :
93 : INTEGER :: handle, ia, idir, ir, my_deriv_order
94 : LOGICAL :: gradient_f, my_epr_xc, my_only_energy
95 : REAL(dp) :: my_adiabatic_rescale_factor
96 46766 : REAL(dp), DIMENSION(:, :, :), POINTER :: deriv_data
97 : REAL(KIND=dp) :: drho_cutoff
98 : TYPE(xc_derivative_type), POINTER :: deriv_att
99 :
100 46766 : CALL timeset(routineN, handle)
101 46766 : my_only_energy = .FALSE.
102 46766 : IF (PRESENT(energy_only)) my_only_energy = energy_only
103 :
104 46766 : IF (PRESENT(adiabatic_rescale_factor)) THEN
105 46766 : my_adiabatic_rescale_factor = adiabatic_rescale_factor
106 : ELSE
107 : my_adiabatic_rescale_factor = 1.0_dp
108 : END IF
109 :
110 : ! needed for the epr routines
111 46766 : my_epr_xc = .FALSE.
112 46766 : IF (PRESENT(epr_xc)) my_epr_xc = epr_xc
113 46766 : my_deriv_order = deriv_order
114 46766 : IF (my_epr_xc) my_deriv_order = 2
115 :
116 : gradient_f = (needs%drho_spin .OR. needs%norm_drho_spin .OR. &
117 46766 : needs%drho .OR. needs%norm_drho)
118 :
119 : ! Calculate the derivatives
120 : CALL xc_functionals_eval(xc_fun_section, &
121 : lsd=lsd, &
122 : rho_set=rho_set, &
123 : deriv_set=deriv_set, &
124 46766 : deriv_order=my_deriv_order)
125 :
126 46766 : CALL xc_rho_set_get(rho_set, drho_cutoff=drho_cutoff)
127 :
128 46766 : NULLIFY (deriv_data)
129 :
130 46766 : IF (my_epr_xc) THEN
131 : ! nabla v_xc (using the vxg arrays)
132 : ! there's no point doing this when lsd = false
133 30 : IF (lsd) THEN
134 30 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa])
135 30 : IF (ASSOCIATED(deriv_att)) THEN
136 30 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
137 1530 : DO ir = 1, nr
138 76530 : DO ia = 1, na
139 301500 : DO idir = 1, 3
140 : vxg(idir, ia, ir, 1) = rho_set%drhoa(idir)%array(ia, ir, 1)* &
141 300000 : deriv_data(ia, ir, 1)
142 : END DO !idir
143 : END DO !ia
144 : END DO !ir
145 30 : NULLIFY (deriv_data)
146 : END IF
147 30 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob])
148 30 : IF (ASSOCIATED(deriv_att)) THEN
149 30 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
150 1530 : DO ir = 1, nr
151 76530 : DO ia = 1, na
152 301500 : DO idir = 1, 3
153 : vxg(idir, ia, ir, 2) = rho_set%drhob(idir)%array(ia, ir, 1)* &
154 300000 : deriv_data(ia, ir, 1)
155 : END DO !idir
156 : END DO !ia
157 : END DO !ir
158 30 : NULLIFY (deriv_data)
159 : END IF
160 : END IF
161 : ! EXC energy ! is that needed for epr?
162 30 : deriv_att => xc_dset_get_derivative(deriv_set, [INTEGER::])
163 30 : exc = 0.0_dp
164 30 : IF (ASSOCIATED(deriv_att)) THEN
165 30 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
166 1530 : DO ir = 1, nr
167 76530 : DO ia = 1, na
168 76500 : exc = exc + deriv_data(ia, ir, 1)*w(ia, ir)
169 : END DO
170 : END DO
171 30 : NULLIFY (deriv_data)
172 : END IF
173 : ELSE
174 : ! EXC energy
175 46736 : deriv_att => xc_dset_get_derivative(deriv_set, [INTEGER::])
176 46736 : exc = 0.0_dp
177 46736 : IF (ASSOCIATED(deriv_att)) THEN
178 46664 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
179 2716944 : DO ir = 1, nr
180 136399424 : DO ia = 1, na
181 136352760 : exc = exc + deriv_data(ia, ir, 1)*w(ia, ir)
182 : END DO
183 : END DO
184 46664 : NULLIFY (deriv_data)
185 : END IF
186 : ! Calculate the potential only if needed
187 46736 : IF (.NOT. my_only_energy) THEN
188 : ! Derivative with respect to the density
189 44688 : IF (lsd) THEN
190 7626 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_rhoa])
191 7626 : IF (ASSOCIATED(deriv_att)) THEN
192 7622 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
193 51173764 : vxc(:, :, 1) = deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
194 7622 : NULLIFY (deriv_data)
195 : END IF
196 7626 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_rhob])
197 7626 : IF (ASSOCIATED(deriv_att)) THEN
198 7622 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
199 51173764 : vxc(:, :, 2) = deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
200 7622 : NULLIFY (deriv_data)
201 : END IF
202 7626 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_rho])
203 7626 : IF (ASSOCIATED(deriv_att)) THEN
204 0 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
205 0 : vxc(:, :, 1) = vxc(:, :, 1) + deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
206 0 : vxc(:, :, 2) = vxc(:, :, 2) + deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
207 0 : NULLIFY (deriv_data)
208 : END IF
209 : ELSE
210 37062 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_rho])
211 37062 : IF (ASSOCIATED(deriv_att)) THEN
212 36994 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
213 209952188 : vxc(:, :, 1) = deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
214 36994 : NULLIFY (deriv_data)
215 : END IF
216 : END IF ! lsd
217 :
218 : ! Derivatives with respect to the gradient
219 44688 : IF (lsd) THEN
220 7626 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa])
221 7626 : IF (ASSOCIATED(deriv_att)) THEN
222 4348 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
223 237028 : DO ir = 1, nr
224 11859508 : DO ia = 1, na
225 46722600 : DO idir = 1, 3
226 46489920 : IF (rho_set%norm_drhoa(ia, ir, 1) > drho_cutoff) THEN
227 : vxg(idir, ia, ir, 1) = rho_set%drhoa(idir)%array(ia, ir, 1)* &
228 : deriv_data(ia, ir, 1)*w(ia, ir)/ &
229 29725674 : rho_set%norm_drhoa(ia, ir, 1)*my_adiabatic_rescale_factor
230 : ELSE
231 5141766 : vxg(idir, ia, ir, 1) = 0.0_dp
232 : END IF
233 : END DO
234 : END DO
235 : END DO
236 4348 : NULLIFY (deriv_data)
237 : END IF
238 7626 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob])
239 7626 : IF (ASSOCIATED(deriv_att)) THEN
240 4348 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
241 237028 : DO ir = 1, nr
242 11859508 : DO ia = 1, na
243 46722600 : DO idir = 1, 3
244 46489920 : IF (rho_set%norm_drhob(ia, ir, 1) > drho_cutoff) THEN
245 : vxg(idir, ia, ir, 2) = rho_set%drhob(idir)%array(ia, ir, 1)* &
246 : deriv_data(ia, ir, 1)*w(ia, ir)/ &
247 29226342 : rho_set%norm_drhob(ia, ir, 1)*my_adiabatic_rescale_factor
248 : ELSE
249 5641098 : vxg(idir, ia, ir, 2) = 0.0_dp
250 : END IF
251 : END DO
252 : END DO
253 : END DO
254 4348 : NULLIFY (deriv_data)
255 : END IF
256 : ! Cross Terms
257 7626 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_norm_drho])
258 7626 : IF (ASSOCIATED(deriv_att)) THEN
259 4300 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
260 234580 : DO ir = 1, nr
261 11737060 : DO ia = 1, na
262 46240200 : DO idir = 1, 3
263 46009920 : IF (rho_set%norm_drho(ia, ir, 1) > drho_cutoff) THEN
264 : vxg(idir, ia, ir, 1:2) = &
265 : vxg(idir, ia, ir, 1:2) + ( &
266 : rho_set%drhoa(idir)%array(ia, ir, 1) + &
267 : rho_set%drhob(idir)%array(ia, ir, 1))* &
268 : deriv_data(ia, ir, 1)*w(ia, ir)/rho_set%norm_drho(ia, ir, 1)* &
269 88328160 : my_adiabatic_rescale_factor
270 : END IF
271 : END DO
272 : END DO
273 : END DO
274 4300 : NULLIFY (deriv_data)
275 : END IF
276 : ELSE
277 37062 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_norm_drho])
278 37062 : IF (ASSOCIATED(deriv_att)) THEN
279 20698 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
280 1075198 : DO ir = 1, nr
281 53980198 : DO ia = 1, na
282 53959500 : IF (rho_set%norm_drho(ia, ir, 1) > drho_cutoff) THEN
283 182730936 : DO idir = 1, 3
284 : vxg(idir, ia, ir, 1) = rho_set%drho(idir)%array(ia, ir, 1)* &
285 : deriv_data(ia, ir, 1)*w(ia, ir)/ &
286 182730936 : rho_set%norm_drho(ia, ir, 1)*my_adiabatic_rescale_factor
287 : END DO
288 : ELSE
289 28889064 : vxg(1:3, ia, ir, 1) = 0.0_dp
290 : END IF
291 : END DO
292 : END DO
293 20698 : NULLIFY (deriv_data)
294 : END IF
295 : END IF ! lsd
296 : ! Derivative with respect to tau
297 44688 : IF (lsd) THEN
298 7626 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_tau_a])
299 7626 : IF (ASSOCIATED(deriv_att)) THEN
300 0 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
301 0 : vtau(:, :, 1) = deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
302 0 : NULLIFY (deriv_data)
303 : END IF
304 7626 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_tau_b])
305 7626 : IF (ASSOCIATED(deriv_att)) THEN
306 0 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
307 0 : vtau(:, :, 2) = deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
308 0 : NULLIFY (deriv_data)
309 : END IF
310 7626 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_tau])
311 7626 : IF (ASSOCIATED(deriv_att)) THEN
312 0 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
313 0 : vtau(:, :, 1) = vtau(:, :, 1) + deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
314 0 : vtau(:, :, 2) = vtau(:, :, 2) + deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
315 0 : NULLIFY (deriv_data)
316 : END IF
317 : ELSE
318 37062 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_tau])
319 37062 : IF (ASSOCIATED(deriv_att)) THEN
320 616 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
321 3782432 : vtau(:, :, 1) = deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
322 616 : NULLIFY (deriv_data)
323 : END IF
324 : END IF ! lsd
325 : END IF ! only_energy
326 : END IF ! epr_xc
327 :
328 46766 : CALL timestop(handle)
329 :
330 46766 : END SUBROUTINE vxc_of_r_new
331 :
332 : ! **************************************************************************************************
333 : !> \brief ...
334 : !> \param rho_set ...
335 : !> \param rho1_set ...
336 : !> \param xc_section ...
337 : !> \param deriv_set ...
338 : !> \param w ...
339 : !> \param vxc ...
340 : !> \param vxg ...
341 : !> \param do_triplet ...
342 : ! **************************************************************************************************
343 9278 : SUBROUTINE xc_2nd_deriv_of_r(rho_set, rho1_set, xc_section, &
344 : deriv_set, w, vxc, vxg, do_triplet)
345 :
346 : ! As input of this routine one gets rho and drho on a one dimensional grid.
347 : ! The grid is the angular grid corresponding to a given point ir on the radial grid.
348 : ! The derivatives are calculated on this one dimensional grid, the results are stored in
349 : ! vxc(1:na,ir,ispin), vxg(1:na,ir,ispin), vxg_cross(1:na,ir,ispin)
350 : ! Afterwords the arrays containing the derivatives are put to zero so that the routine
351 : ! can safely be called for the next radial point ir
352 :
353 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set, rho1_set
354 : TYPE(section_vals_type), POINTER :: xc_section
355 : TYPE(xc_derivative_set_type), INTENT(INOUT) :: deriv_set
356 : REAL(dp), DIMENSION(:, :), POINTER :: w
357 : REAL(dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: vxc
358 : REAL(dp), DIMENSION(:, :, :, :), POINTER :: vxg
359 : LOGICAL, INTENT(IN), OPTIONAL :: do_triplet
360 :
361 : CHARACTER(LEN=*), PARAMETER :: routineN = 'xc_2nd_deriv_of_r'
362 :
363 : INTEGER :: handle, ispin, nspins
364 : LOGICAL :: lsd
365 : REAL(dp) :: drho_cutoff, my_fac_triplet
366 : TYPE(cp_sll_xc_deriv_type), POINTER :: pos
367 : TYPE(pw_pool_type), POINTER :: pw_pool
368 9278 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: vxc_pw, vxc_tau_pw
369 : TYPE(section_vals_type), POINTER :: xc_fun_section
370 : TYPE(xc_derivative_type), POINTER :: deriv_att
371 :
372 9278 : CALL timeset(routineN, handle)
373 :
374 9278 : nspins = SIZE(vxc, 3)
375 9278 : lsd = (nspins == 2)
376 9278 : IF (ASSOCIATED(rho_set%rhoa)) THEN
377 662 : lsd = .TRUE.
378 : END IF
379 9278 : my_fac_triplet = 1.0_dp
380 9278 : IF (PRESENT(do_triplet)) THEN
381 9082 : IF (do_triplet) my_fac_triplet = -1.0_dp
382 : END IF
383 :
384 9278 : CALL xc_rho_set_get(rho_set, drho_cutoff=drho_cutoff)
385 : xc_fun_section => section_vals_get_subs_vals(xc_section, &
386 9278 : "XC_FUNCTIONAL")
387 :
388 : ! Calculate the derivatives
389 : CALL xc_functionals_eval(xc_fun_section, &
390 : lsd=lsd, &
391 : rho_set=rho_set, &
392 : deriv_set=deriv_set, &
393 9278 : deriv_order=2)
394 :
395 9278 : CALL divide_by_norm_drho(deriv_set, rho_set, lsd)
396 :
397 : ! multiply by w
398 9278 : pos => deriv_set%derivs
399 61864 : DO WHILE (cp_sll_xc_deriv_next(pos, el_att=deriv_att))
400 268303050 : deriv_att%deriv_data(:, :, 1) = w(:, :)*deriv_att%deriv_data(:, :, 1)
401 : END DO
402 :
403 9278 : NULLIFY (pw_pool)
404 37240 : ALLOCATE (vxc_pw(nspins))
405 18684 : DO ispin = 1, nspins
406 18684 : vxc_pw(ispin)%array => vxc(:, :, ispin:ispin)
407 : END DO
408 :
409 9278 : NULLIFY (vxc_tau_pw)
410 :
411 : CALL xc_calc_2nd_deriv_analytical(vxc_pw, vxc_tau_pw, deriv_set, rho_set, rho1_set, pw_pool, &
412 9278 : xc_section, gapw=.TRUE., vxg=vxg, tddfpt_fac=my_fac_triplet)
413 :
414 9278 : DEALLOCATE (vxc_pw)
415 :
416 : ! zero the derivative data for the next call
417 9278 : pos => deriv_set%derivs
418 61864 : DO WHILE (cp_sll_xc_deriv_next(pos, el_att=deriv_att))
419 134208750 : deriv_att%deriv_data = 0.0_dp
420 : END DO
421 :
422 9278 : CALL timestop(handle)
423 :
424 18556 : END SUBROUTINE xc_2nd_deriv_of_r
425 :
426 : ! **************************************************************************************************
427 : !> \brief ...
428 : !> \param rho_set ...
429 : !> \param needs ...
430 : !> \param nspins ...
431 : !> \param bo ...
432 : ! **************************************************************************************************
433 105926 : SUBROUTINE xc_rho_set_atom_update(rho_set, needs, nspins, bo)
434 :
435 : ! This routine allocates the storage arrays for rho and drho
436 : ! In calculate_vxc_atom this is called once for each atomic_kind,
437 : ! After the loop over all the atoms of the kind and over all the points
438 : ! of the radial grid for each atom, rho_set is deallocated.
439 : ! Within the same kind, at each new point on the radial grid, the rho_set
440 : ! arrays rho and drho are overwritten.
441 :
442 : TYPE(xc_rho_set_type), INTENT(INOUT) :: rho_set
443 : TYPE(xc_rho_cflags_type), INTENT(IN) :: needs
444 : INTEGER, INTENT(IN) :: nspins
445 : INTEGER, DIMENSION(2, 3), INTENT(IN) :: bo
446 :
447 : INTEGER :: idir
448 :
449 198053 : SELECT CASE (nspins)
450 : CASE (1)
451 : ! What is this for?
452 92127 : IF (needs%rho_1_3) THEN
453 4007 : NULLIFY (rho_set%rho_1_3)
454 20035 : ALLOCATE (rho_set%rho_1_3(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
455 4007 : rho_set%owns%rho_1_3 = .TRUE.
456 4007 : rho_set%has%rho_1_3 = .FALSE.
457 : END IF
458 : ! Allocate the storage space for the density
459 92127 : IF (needs%rho) THEN
460 92127 : NULLIFY (rho_set%rho)
461 460635 : ALLOCATE (rho_set%rho(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
462 92127 : rho_set%owns%rho = .TRUE.
463 92127 : rho_set%has%rho = .FALSE.
464 : END IF
465 : ! Allocate the storage space for the norm of the gradient of the density
466 92127 : IF (needs%norm_drho) THEN
467 65817 : NULLIFY (rho_set%norm_drho)
468 329085 : ALLOCATE (rho_set%norm_drho(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
469 65817 : rho_set%owns%norm_drho = .TRUE.
470 65817 : rho_set%has%norm_drho = .FALSE.
471 : END IF
472 : ! Allocate the storage space for the three components of the gradient of the density
473 92127 : IF (needs%drho) THEN
474 178192 : DO idir = 1, 3
475 133644 : NULLIFY (rho_set%drho(idir)%array)
476 712768 : ALLOCATE (rho_set%drho(idir)%array(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
477 : END DO
478 44548 : rho_set%owns%drho = .TRUE.
479 44548 : rho_set%has%drho = .FALSE.
480 : END IF
481 : CASE (2)
482 : ! Allocate the storage space for the total density
483 13799 : IF (needs%rho) THEN
484 : ! this should never be the case unless you use LDA functionals with LSD
485 0 : NULLIFY (rho_set%rho)
486 0 : ALLOCATE (rho_set%rho(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
487 0 : rho_set%owns%rho = .TRUE.
488 0 : rho_set%has%rho = .FALSE.
489 : END IF
490 : ! What is this for?
491 13799 : IF (needs%rho_1_3) THEN
492 0 : NULLIFY (rho_set%rho_1_3)
493 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)))
494 0 : rho_set%owns%rho_1_3 = .TRUE.
495 0 : rho_set%has%rho_1_3 = .FALSE.
496 : END IF
497 : ! What is this for?
498 13799 : IF (needs%rho_spin_1_3) THEN
499 2456 : NULLIFY (rho_set%rhoa_1_3, rho_set%rhob_1_3)
500 12280 : ALLOCATE (rho_set%rhoa_1_3(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
501 12280 : ALLOCATE (rho_set%rhob_1_3(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
502 2456 : rho_set%owns%rho_spin_1_3 = .TRUE.
503 2456 : rho_set%has%rho_spin_1_3 = .FALSE.
504 : END IF
505 : ! Allocate the storage space for the spin densities rhoa and rhob
506 13799 : IF (needs%rho_spin) THEN
507 13799 : NULLIFY (rho_set%rhoa, rho_set%rhob)
508 68995 : ALLOCATE (rho_set%rhoa(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
509 68995 : ALLOCATE (rho_set%rhob(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
510 13799 : rho_set%owns%rho_spin = .TRUE.
511 13799 : rho_set%has%rho_spin = .FALSE.
512 : END IF
513 : ! Allocate the storage space for the norm of the gradient of the total density
514 13799 : IF (needs%norm_drho) THEN
515 8469 : NULLIFY (rho_set%norm_drho)
516 42345 : ALLOCATE (rho_set%norm_drho(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
517 8469 : rho_set%owns%norm_drho = .TRUE.
518 8469 : rho_set%has%norm_drho = .FALSE.
519 : END IF
520 : ! Allocate the storage space for the norm of the gradient of rhoa and of rhob separatedly
521 13799 : IF (needs%norm_drho_spin) THEN
522 8517 : NULLIFY (rho_set%norm_drhoa, rho_set%norm_drhob)
523 42585 : ALLOCATE (rho_set%norm_drhoa(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
524 42585 : ALLOCATE (rho_set%norm_drhob(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
525 8517 : rho_set%owns%norm_drho_spin = .TRUE.
526 8517 : rho_set%has%norm_drho_spin = .FALSE.
527 : END IF
528 : ! Allocate the storage space for the components of the gradient for the total rho
529 13799 : IF (needs%drho) THEN
530 0 : DO idir = 1, 3
531 0 : NULLIFY (rho_set%drho(idir)%array)
532 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)))
533 : END DO
534 0 : rho_set%owns%drho = .TRUE.
535 0 : rho_set%has%drho = .FALSE.
536 : END IF
537 : ! Allocate the storage space for the components of the gradient for rhoa and rhob
538 119725 : IF (needs%drho_spin) THEN
539 30552 : DO idir = 1, 3
540 22914 : NULLIFY (rho_set%drhoa(idir)%array, rho_set%drhob(idir)%array)
541 114570 : ALLOCATE (rho_set%drhoa(idir)%array(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
542 122208 : ALLOCATE (rho_set%drhob(idir)%array(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
543 : END DO
544 7638 : rho_set%owns%drho_spin = .TRUE.
545 7638 : rho_set%has%drho_spin = .FALSE.
546 : END IF
547 : !
548 : END SELECT
549 :
550 : ! tau part
551 105926 : IF (needs%tau) THEN
552 996 : NULLIFY (rho_set%tau)
553 4980 : ALLOCATE (rho_set%tau(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
554 996 : rho_set%owns%tau = .TRUE.
555 : END IF
556 105926 : IF (needs%tau_spin) THEN
557 2 : NULLIFY (rho_set%tau_a, rho_set%tau_b)
558 10 : ALLOCATE (rho_set%tau_a(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
559 10 : ALLOCATE (rho_set%tau_b(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
560 2 : rho_set%owns%tau_spin = .TRUE.
561 2 : rho_set%has%tau_spin = .FALSE.
562 : END IF
563 :
564 : ! Laplace part
565 105926 : IF (needs%laplace_rho) THEN
566 0 : NULLIFY (rho_set%laplace_rho)
567 0 : ALLOCATE (rho_set%laplace_rho(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
568 0 : rho_set%owns%laplace_rho = .TRUE.
569 : END IF
570 105926 : IF (needs%laplace_rho_spin) THEN
571 0 : NULLIFY (rho_set%laplace_rhoa)
572 0 : NULLIFY (rho_set%laplace_rhob)
573 0 : ALLOCATE (rho_set%laplace_rhoa(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
574 0 : ALLOCATE (rho_set%laplace_rhob(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
575 0 : rho_set%owns%laplace_rho_spin = .TRUE.
576 0 : rho_set%has%laplace_rho_spin = .TRUE.
577 : END IF
578 :
579 105926 : END SUBROUTINE xc_rho_set_atom_update
580 :
581 : ! **************************************************************************************************
582 : !> \brief ...
583 : !> \param rho_set ...
584 : !> \param lsd ...
585 : !> \param nspins ...
586 : !> \param needs ...
587 : !> \param rho ...
588 : !> \param drho ...
589 : !> \param tau ...
590 : !> \param na ...
591 : !> \param ir ...
592 : ! **************************************************************************************************
593 3553980 : SUBROUTINE fill_rho_set(rho_set, lsd, nspins, needs, rho, drho, tau, na, ir)
594 :
595 : TYPE(xc_rho_set_type), INTENT(INOUT) :: rho_set
596 : LOGICAL, INTENT(IN) :: lsd
597 : INTEGER, INTENT(IN) :: nspins
598 : TYPE(xc_rho_cflags_type), INTENT(IN) :: needs
599 : REAL(dp), DIMENSION(:, :, :), INTENT(IN) :: rho
600 : REAL(dp), DIMENSION(:, :, :, :), INTENT(IN) :: drho
601 : REAL(dp), DIMENSION(:, :, :), INTENT(IN) :: tau
602 : INTEGER, INTENT(IN) :: na, ir
603 :
604 : REAL(KIND=dp), PARAMETER :: f13 = (1.0_dp/3.0_dp)
605 :
606 : INTEGER :: ia, idir, my_nspins
607 : LOGICAL :: gradient_f, tddft_split
608 :
609 3553980 : my_nspins = nspins
610 3553980 : tddft_split = .FALSE.
611 3553980 : IF (lsd .AND. nspins == 1) THEN
612 48000 : my_nspins = 2
613 48000 : tddft_split = .TRUE.
614 : END IF
615 :
616 : ! some checks
617 3553980 : IF (lsd) THEN
618 : ELSE
619 2942300 : CPASSERT(SIZE(rho, 3) == 1)
620 : END IF
621 2942300 : SELECT CASE (my_nspins)
622 : CASE (1)
623 2942300 : CPASSERT(.NOT. needs%rho_spin)
624 2942300 : CPASSERT(.NOT. needs%drho_spin)
625 2942300 : CPASSERT(.NOT. needs%norm_drho_spin)
626 2942300 : CPASSERT(.NOT. needs%rho_spin_1_3)
627 : CASE (2)
628 : CASE default
629 3553980 : CPABORT("Unsupported number of spins")
630 : END SELECT
631 :
632 : gradient_f = (needs%drho_spin .OR. needs%norm_drho_spin .OR. &
633 3553980 : needs%drho .OR. needs%norm_drho)
634 :
635 2942300 : SELECT CASE (my_nspins)
636 : CASE (1)
637 : ! Give rho to 1/3
638 2942300 : IF (needs%rho_1_3) THEN
639 6684000 : DO ia = 1, na
640 6684000 : rho_set%rho_1_3(ia, ir, 1) = MAX(rho(ia, ir, 1), 0.0_dp)**f13
641 : END DO
642 131200 : rho_set%owns%rho_1_3 = .TRUE.
643 131200 : rho_set%has%rho_1_3 = .TRUE.
644 : END IF
645 : ! Give the density
646 2942300 : IF (needs%rho) THEN
647 150237300 : DO ia = 1, na
648 150237300 : rho_set%rho(ia, ir, 1) = rho(ia, ir, 1)
649 : END DO
650 2942300 : rho_set%owns%rho = .TRUE.
651 2942300 : rho_set%has%rho = .TRUE.
652 : END IF
653 : ! Give the norm of the gradient of the density
654 2942300 : IF (needs%norm_drho) THEN
655 87017700 : DO ia = 1, na
656 87017700 : rho_set%norm_drho(ia, ir, 1) = drho(4, ia, ir, 1)
657 : END DO
658 1702700 : rho_set%owns%norm_drho = .TRUE.
659 1702700 : rho_set%has%norm_drho = .TRUE.
660 : END IF
661 : ! Give the three components of the gradient of the density
662 2942300 : IF (needs%drho) THEN
663 6810800 : DO idir = 1, 3
664 262755800 : DO ia = 1, na
665 261053100 : rho_set%drho(idir)%array(ia, ir, 1) = drho(idir, ia, ir, 1)
666 : END DO
667 : END DO
668 1702700 : rho_set%owns%drho = .TRUE.
669 1702700 : rho_set%has%drho = .TRUE.
670 : END IF
671 : CASE (2)
672 : ! Give the total density
673 611680 : IF (needs%rho) THEN
674 : ! this should never be the case unless you use LDA functionals with LSD
675 0 : IF (.NOT. tddft_split) THEN
676 0 : DO ia = 1, na
677 0 : rho_set%rho(ia, ir, 1) = rho(ia, ir, 1) + rho(ia, ir, 2)
678 : END DO
679 : ELSE
680 0 : DO ia = 1, na
681 0 : rho_set%rho(ia, ir, 1) = rho(ia, ir, 1)
682 : END DO
683 : END IF
684 0 : rho_set%owns%rho = .TRUE.
685 0 : rho_set%has%rho = .TRUE.
686 : END IF
687 : ! Give the total density to 1/3
688 611680 : IF (needs%rho_1_3) THEN
689 0 : IF (.NOT. tddft_split) THEN
690 0 : DO ia = 1, na
691 0 : rho_set%rho_1_3(ia, ir, 1) = MAX(rho(ia, ir, 1) + rho(ia, ir, 2), 0.0_dp)**f13
692 : END DO
693 : ELSE
694 0 : DO ia = 1, na
695 0 : rho_set%rho_1_3(ia, ir, 1) = MAX(rho(ia, ir, 1), 0.0_dp)**f13
696 : END DO
697 : END IF
698 0 : rho_set%owns%rho_1_3 = .TRUE.
699 0 : rho_set%has%rho_1_3 = .TRUE.
700 : END IF
701 : ! Give the spin densities to 1/3
702 611680 : IF (needs%rho_spin_1_3) THEN
703 75880 : IF (.NOT. tddft_split) THEN
704 3858360 : DO ia = 1, na
705 3782480 : rho_set%rhoa_1_3(ia, ir, 1) = MAX(rho(ia, ir, 1), 0.0_dp)**f13
706 3858360 : rho_set%rhob_1_3(ia, ir, 1) = MAX(rho(ia, ir, 2), 0.0_dp)**f13
707 : END DO
708 : ELSE
709 0 : DO ia = 1, na
710 0 : rho_set%rhoa_1_3(ia, ir, 1) = MAX(0.5_dp*rho(ia, ir, 1), 0.0_dp)**f13
711 0 : rho_set%rhob_1_3(ia, ir, 1) = rho_set%rhoa_1_3(ia, ir, 1)
712 : END DO
713 : END IF
714 75880 : rho_set%owns%rho_spin_1_3 = .TRUE.
715 75880 : rho_set%has%rho_spin_1_3 = .TRUE.
716 : END IF
717 : ! Give the spin densities rhoa and rhob
718 611680 : IF (needs%rho_spin) THEN
719 611680 : IF (.NOT. tddft_split) THEN
720 28736160 : DO ia = 1, na
721 28172480 : rho_set%rhoa(ia, ir, 1) = rho(ia, ir, 1)
722 28736160 : rho_set%rhob(ia, ir, 1) = rho(ia, ir, 2)
723 : END DO
724 : ELSE
725 2448000 : DO ia = 1, na
726 2400000 : rho_set%rhoa(ia, ir, 1) = 0.5_dp*rho(ia, ir, 1)
727 2448000 : rho_set%rhob(ia, ir, 1) = rho_set%rhoa(ia, ir, 1)
728 : END DO
729 : END IF
730 611680 : rho_set%owns%rho_spin = .TRUE.
731 611680 : rho_set%has%rho_spin = .TRUE.
732 : END IF
733 : ! Give the norm of the gradient of the total density
734 611680 : IF (needs%norm_drho) THEN
735 308480 : IF (.NOT. tddft_split) THEN
736 13884960 : DO ia = 1, na
737 : rho_set%norm_drho(ia, ir, 1) = SQRT( &
738 : (drho(1, ia, ir, 1) + drho(1, ia, ir, 2))**2 + &
739 : (drho(2, ia, ir, 1) + drho(2, ia, ir, 2))**2 + &
740 13884960 : (drho(3, ia, ir, 1) + drho(3, ia, ir, 2))**2)
741 : END DO
742 : ELSE
743 1836000 : DO ia = 1, na
744 1836000 : rho_set%norm_drho(ia, ir, 1) = drho(4, ia, ir, 1)
745 : END DO
746 : END IF
747 308480 : rho_set%owns%norm_drho = .TRUE.
748 308480 : rho_set%has%norm_drho = .TRUE.
749 : END IF
750 : ! Give the norm of the gradient of rhoa and of rhob separatedly
751 611680 : IF (needs%norm_drho_spin) THEN
752 310880 : IF (.NOT. tddft_split) THEN
753 14007360 : DO ia = 1, na
754 13732480 : rho_set%norm_drhoa(ia, ir, 1) = drho(4, ia, ir, 1)
755 14007360 : rho_set%norm_drhob(ia, ir, 1) = drho(4, ia, ir, 2)
756 : END DO
757 : ELSE
758 1836000 : DO ia = 1, na
759 1800000 : rho_set%norm_drhoa(ia, ir, 1) = 0.5_dp*drho(4, ia, ir, 1)
760 1836000 : rho_set%norm_drhob(ia, ir, 1) = rho_set%norm_drhoa(ia, ir, 1)
761 : END DO
762 : END IF
763 310880 : rho_set%owns%norm_drho_spin = .TRUE.
764 310880 : rho_set%has%norm_drho_spin = .TRUE.
765 : END IF
766 : ! Give the components of the gradient for the total rho
767 611680 : IF (needs%drho) THEN
768 0 : IF (.NOT. tddft_split) THEN
769 0 : DO idir = 1, 3
770 0 : DO ia = 1, na
771 0 : rho_set%drho(idir)%array(ia, ir, 1) = drho(idir, ia, ir, 1) + drho(idir, ia, ir, 2)
772 : END DO
773 : END DO
774 : ELSE
775 0 : DO idir = 1, 3
776 0 : DO ia = 1, na
777 0 : rho_set%drho(idir)%array(ia, ir, 1) = drho(idir, ia, ir, 1)
778 : END DO
779 : END DO
780 : END IF
781 0 : rho_set%owns%drho = .TRUE.
782 0 : rho_set%has%drho = .TRUE.
783 : END IF
784 : ! Give the components of the gradient for rhoa and rhob
785 4165660 : IF (needs%drho_spin) THEN
786 312380 : IF (.NOT. tddft_split) THEN
787 1105520 : DO idir = 1, 3
788 42527960 : DO ia = 1, na
789 41422440 : rho_set%drhoa(idir)%array(ia, ir, 1) = drho(idir, ia, ir, 1)
790 42251580 : rho_set%drhob(idir)%array(ia, ir, 1) = drho(idir, ia, ir, 2)
791 : END DO
792 : END DO
793 : ELSE
794 144000 : DO idir = 1, 3
795 5544000 : DO ia = 1, na
796 5400000 : rho_set%drhoa(idir)%array(ia, ir, 1) = 0.5_dp*drho(idir, ia, ir, 1)
797 5508000 : rho_set%drhob(idir)%array(ia, ir, 1) = rho_set%drhoa(idir)%array(ia, ir, 1)
798 : END DO
799 : END DO
800 : END IF
801 312380 : rho_set%owns%drho_spin = .TRUE.
802 312380 : rho_set%has%drho_spin = .TRUE.
803 : END IF
804 : !
805 : END SELECT
806 :
807 : ! tau part
808 3553980 : IF (needs%tau .OR. needs%tau_spin) THEN
809 3553980 : CPASSERT(SIZE(tau, 3) == my_nspins)
810 : END IF
811 3553980 : IF (needs%tau) THEN
812 33400 : IF (my_nspins == 2) THEN
813 0 : DO ia = 1, na
814 0 : rho_set%tau(ia, ir, 1) = tau(ia, ir, 1) + tau(ia, ir, 2)
815 : END DO
816 0 : rho_set%owns%tau = .TRUE.
817 0 : rho_set%has%tau = .TRUE.
818 : ELSE
819 1890600 : DO ia = 1, na
820 1890600 : rho_set%tau(ia, ir, 1) = tau(ia, ir, 1)
821 : END DO
822 33400 : rho_set%owns%tau = .TRUE.
823 33400 : rho_set%has%tau = .TRUE.
824 : END IF
825 : END IF
826 3553980 : IF (needs%tau_spin) THEN
827 0 : DO ia = 1, na
828 0 : rho_set%tau_a(ia, ir, 1) = tau(ia, ir, 1)
829 0 : rho_set%tau_b(ia, ir, 1) = tau(ia, ir, 2)
830 : END DO
831 0 : rho_set%owns%tau_spin = .TRUE.
832 0 : rho_set%has%tau_spin = .TRUE.
833 : END IF
834 :
835 3553980 : END SUBROUTINE fill_rho_set
836 :
837 : END MODULE xc_atom
|
