Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 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 93324 : 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 46662 : REAL(dp), DIMENSION(:, :, :), POINTER :: deriv_data
97 : REAL(KIND=dp) :: drho_cutoff
98 : TYPE(xc_derivative_type), POINTER :: deriv_att
99 :
100 46662 : CALL timeset(routineN, handle)
101 46662 : my_only_energy = .FALSE.
102 46662 : IF (PRESENT(energy_only)) my_only_energy = energy_only
103 :
104 46662 : IF (PRESENT(adiabatic_rescale_factor)) THEN
105 46662 : 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 46662 : my_epr_xc = .FALSE.
112 46662 : IF (PRESENT(epr_xc)) my_epr_xc = epr_xc
113 46662 : my_deriv_order = deriv_order
114 46662 : IF (my_epr_xc) my_deriv_order = 2
115 :
116 : gradient_f = (needs%drho_spin .OR. needs%norm_drho_spin .OR. &
117 46662 : 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 46662 : deriv_order=my_deriv_order)
125 :
126 46662 : CALL xc_rho_set_get(rho_set, drho_cutoff=drho_cutoff)
127 :
128 46662 : NULLIFY (deriv_data)
129 :
130 46662 : 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 46632 : deriv_att => xc_dset_get_derivative(deriv_set, [INTEGER::])
176 46632 : exc = 0.0_dp
177 46632 : IF (ASSOCIATED(deriv_att)) THEN
178 46560 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
179 2711640 : DO ir = 1, nr
180 136134120 : DO ia = 1, na
181 136087560 : exc = exc + deriv_data(ia, ir, 1)*w(ia, ir)
182 : END DO
183 : END DO
184 46560 : NULLIFY (deriv_data)
185 : END IF
186 : ! Calculate the potential only if needed
187 46632 : IF (.NOT. my_only_energy) THEN
188 : ! Derivative with respect to the density
189 44568 : IF (lsd) THEN
190 7640 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_rhoa])
191 7640 : IF (ASSOCIATED(deriv_att)) THEN
192 7636 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
193 51245192 : vxc(:, :, 1) = deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
194 7636 : NULLIFY (deriv_data)
195 : END IF
196 7640 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_rhob])
197 7640 : IF (ASSOCIATED(deriv_att)) THEN
198 7636 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
199 51245192 : vxc(:, :, 2) = deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
200 7636 : NULLIFY (deriv_data)
201 : END IF
202 7640 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_rho])
203 7640 : 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 36928 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_rho])
211 36928 : IF (ASSOCIATED(deriv_att)) THEN
212 36860 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
213 209268520 : vxc(:, :, 1) = deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
214 36860 : NULLIFY (deriv_data)
215 : END IF
216 : END IF ! lsd
217 :
218 : ! Derivatives with respect to the gradient
219 44568 : IF (lsd) THEN
220 7640 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa])
221 7640 : IF (ASSOCIATED(deriv_att)) THEN
222 4374 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
223 238354 : DO ir = 1, nr
224 11925834 : DO ia = 1, na
225 46983900 : DO idir = 1, 3
226 46749920 : 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 29887218 : rho_set%norm_drhoa(ia, ir, 1)*my_adiabatic_rescale_factor
230 : ELSE
231 5175222 : vxg(idir, ia, ir, 1) = 0.0_dp
232 : END IF
233 : END DO
234 : END DO
235 : END DO
236 4374 : NULLIFY (deriv_data)
237 : END IF
238 7640 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob])
239 7640 : IF (ASSOCIATED(deriv_att)) THEN
240 4374 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
241 238354 : DO ir = 1, nr
242 11925834 : DO ia = 1, na
243 46983900 : DO idir = 1, 3
244 46749920 : 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 29388486 : rho_set%norm_drhob(ia, ir, 1)*my_adiabatic_rescale_factor
248 : ELSE
249 5673954 : vxg(idir, ia, ir, 2) = 0.0_dp
250 : END IF
251 : END DO
252 : END DO
253 : END DO
254 4374 : NULLIFY (deriv_data)
255 : END IF
256 : ! Cross Terms
257 7640 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_norm_drho])
258 7640 : IF (ASSOCIATED(deriv_att)) THEN
259 4326 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
260 235906 : DO ir = 1, nr
261 11803386 : DO ia = 1, na
262 46501500 : DO idir = 1, 3
263 46269920 : 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 88817004 : my_adiabatic_rescale_factor
270 : END IF
271 : END DO
272 : END DO
273 : END DO
274 4326 : NULLIFY (deriv_data)
275 : END IF
276 : ELSE
277 36928 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_norm_drho])
278 36928 : IF (ASSOCIATED(deriv_att)) THEN
279 20428 : CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
280 1061428 : DO ir = 1, nr
281 53291428 : DO ia = 1, na
282 53271000 : IF (rho_set%norm_drho(ia, ir, 1) > drho_cutoff) THEN
283 180056936 : 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 180056936 : rho_set%norm_drho(ia, ir, 1)*my_adiabatic_rescale_factor
287 : END DO
288 : ELSE
289 28863064 : vxg(1:3, ia, ir, 1) = 0.0_dp
290 : END IF
291 : END DO
292 : END DO
293 20428 : NULLIFY (deriv_data)
294 : END IF
295 : END IF ! lsd
296 : ! Derivative with respect to tau
297 44568 : IF (lsd) THEN
298 7640 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_tau_a])
299 7640 : 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 7640 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_tau_b])
305 7640 : 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 7640 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_tau])
311 7640 : 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 36928 : deriv_att => xc_dset_get_derivative(deriv_set, [deriv_tau])
319 36928 : 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 46662 : CALL timestop(handle)
329 :
330 46662 : 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 10060 : 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 10060 : 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 10060 : CALL timeset(routineN, handle)
373 :
374 10060 : nspins = SIZE(vxc, 3)
375 10060 : lsd = (nspins == 2)
376 10060 : IF (ASSOCIATED(rho_set%rhoa)) THEN
377 776 : lsd = .TRUE.
378 : END IF
379 10060 : my_fac_triplet = 1.0_dp
380 10060 : IF ((PRESENT(do_triplet)) .AND. (do_triplet)) my_fac_triplet = -1.0_dp
381 :
382 10060 : CALL xc_rho_set_get(rho_set, drho_cutoff=drho_cutoff)
383 : xc_fun_section => section_vals_get_subs_vals(xc_section, &
384 10060 : "XC_FUNCTIONAL")
385 :
386 : ! Calculate the derivatives
387 : CALL xc_functionals_eval(xc_fun_section, &
388 : lsd=lsd, &
389 : rho_set=rho_set, &
390 : deriv_set=deriv_set, &
391 10060 : deriv_order=2)
392 :
393 10060 : CALL divide_by_norm_drho(deriv_set, rho_set, lsd)
394 :
395 : ! multiply by w
396 10060 : pos => deriv_set%derivs
397 67236 : DO WHILE (cp_sll_xc_deriv_next(pos, el_att=deriv_att))
398 291722012 : deriv_att%deriv_data(:, :, 1) = w(:, :)*deriv_att%deriv_data(:, :, 1)
399 : END DO
400 :
401 10060 : NULLIFY (pw_pool)
402 40368 : ALLOCATE (vxc_pw(nspins))
403 20248 : DO ispin = 1, nspins
404 20248 : vxc_pw(ispin)%array => vxc(:, :, ispin:ispin)
405 : END DO
406 :
407 10060 : NULLIFY (vxc_tau_pw)
408 :
409 : CALL xc_calc_2nd_deriv_analytical(vxc_pw, vxc_tau_pw, deriv_set, rho_set, rho1_set, pw_pool, &
410 10060 : xc_section, gapw=.TRUE., vxg=vxg, tddfpt_fac=my_fac_triplet)
411 :
412 10060 : DEALLOCATE (vxc_pw)
413 :
414 : ! zero the derivative data for the next call
415 10060 : pos => deriv_set%derivs
416 67236 : DO WHILE (cp_sll_xc_deriv_next(pos, el_att=deriv_att))
417 145923212 : deriv_att%deriv_data = 0.0_dp
418 : END DO
419 :
420 10060 : CALL timestop(handle)
421 :
422 20120 : END SUBROUTINE xc_2nd_deriv_of_r
423 :
424 : ! **************************************************************************************************
425 : !> \brief ...
426 : !> \param rho_set ...
427 : !> \param needs ...
428 : !> \param nspins ...
429 : !> \param bo ...
430 : ! **************************************************************************************************
431 107367 : SUBROUTINE xc_rho_set_atom_update(rho_set, needs, nspins, bo)
432 :
433 : ! This routine allocates the storage arrays for rho and drho
434 : ! In calculate_vxc_atom this is called once for each atomic_kind,
435 : ! After the loop over all the atoms of the kind and over all the points
436 : ! of the radial grid for each atom, rho_set is deallocated.
437 : ! Within the same kind, at each new point on the radial grid, the rho_set
438 : ! arrays rho and drho are overwritten.
439 :
440 : TYPE(xc_rho_set_type), INTENT(INOUT) :: rho_set
441 : TYPE(xc_rho_cflags_type), INTENT(IN) :: needs
442 : INTEGER, INTENT(IN) :: nspins
443 : INTEGER, DIMENSION(2, 3), INTENT(IN) :: bo
444 :
445 : INTEGER :: idir
446 :
447 200583 : SELECT CASE (nspins)
448 : CASE (1)
449 : ! What is this for?
450 93216 : IF (needs%rho_1_3) THEN
451 4087 : NULLIFY (rho_set%rho_1_3)
452 20435 : ALLOCATE (rho_set%rho_1_3(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
453 4087 : rho_set%owns%rho_1_3 = .TRUE.
454 4087 : rho_set%has%rho_1_3 = .FALSE.
455 : END IF
456 : ! Allocate the storage space for the density
457 93216 : IF (needs%rho) THEN
458 93216 : NULLIFY (rho_set%rho)
459 466080 : ALLOCATE (rho_set%rho(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
460 93216 : rho_set%owns%rho = .TRUE.
461 93216 : rho_set%has%rho = .FALSE.
462 : END IF
463 : ! Allocate the storage space for the norm of the gradient of the density
464 93216 : IF (needs%norm_drho) THEN
465 66194 : NULLIFY (rho_set%norm_drho)
466 330970 : ALLOCATE (rho_set%norm_drho(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
467 66194 : rho_set%owns%norm_drho = .TRUE.
468 66194 : rho_set%has%norm_drho = .FALSE.
469 : END IF
470 : ! Allocate the storage space for the three components of the gradient of the density
471 93216 : IF (needs%drho) THEN
472 182512 : DO idir = 1, 3
473 136884 : NULLIFY (rho_set%drho(idir)%array)
474 730048 : ALLOCATE (rho_set%drho(idir)%array(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
475 : END DO
476 45628 : rho_set%owns%drho = .TRUE.
477 45628 : rho_set%has%drho = .FALSE.
478 : END IF
479 : CASE (2)
480 : ! Allocate the storage space for the total density
481 14151 : IF (needs%rho) THEN
482 : ! this should never be the case unless you use LDA functionals with LSD
483 0 : NULLIFY (rho_set%rho)
484 0 : ALLOCATE (rho_set%rho(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
485 0 : rho_set%owns%rho = .TRUE.
486 0 : rho_set%has%rho = .FALSE.
487 : END IF
488 : ! What is this for?
489 14151 : IF (needs%rho_1_3) THEN
490 0 : NULLIFY (rho_set%rho_1_3)
491 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)))
492 0 : rho_set%owns%rho_1_3 = .TRUE.
493 0 : rho_set%has%rho_1_3 = .FALSE.
494 : END IF
495 : ! What is this for?
496 14151 : IF (needs%rho_spin_1_3) THEN
497 2452 : NULLIFY (rho_set%rhoa_1_3, rho_set%rhob_1_3)
498 12260 : ALLOCATE (rho_set%rhoa_1_3(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
499 12260 : ALLOCATE (rho_set%rhob_1_3(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
500 2452 : rho_set%owns%rho_spin_1_3 = .TRUE.
501 2452 : rho_set%has%rho_spin_1_3 = .FALSE.
502 : END IF
503 : ! Allocate the storage space for the spin densities rhoa and rhob
504 14151 : IF (needs%rho_spin) THEN
505 14151 : NULLIFY (rho_set%rhoa, rho_set%rhob)
506 70755 : ALLOCATE (rho_set%rhoa(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
507 70755 : ALLOCATE (rho_set%rhob(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
508 14151 : rho_set%owns%rho_spin = .TRUE.
509 14151 : rho_set%has%rho_spin = .FALSE.
510 : END IF
511 : ! Allocate the storage space for the norm of the gradient of the total density
512 14151 : IF (needs%norm_drho) THEN
513 8749 : NULLIFY (rho_set%norm_drho)
514 43745 : ALLOCATE (rho_set%norm_drho(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
515 8749 : rho_set%owns%norm_drho = .TRUE.
516 8749 : rho_set%has%norm_drho = .FALSE.
517 : END IF
518 : ! Allocate the storage space for the norm of the gradient of rhoa and of rhob separatedly
519 14151 : IF (needs%norm_drho_spin) THEN
520 8797 : NULLIFY (rho_set%norm_drhoa, rho_set%norm_drhob)
521 43985 : ALLOCATE (rho_set%norm_drhoa(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
522 43985 : ALLOCATE (rho_set%norm_drhob(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
523 8797 : rho_set%owns%norm_drho_spin = .TRUE.
524 8797 : rho_set%has%norm_drho_spin = .FALSE.
525 : END IF
526 : ! Allocate the storage space for the components of the gradient for the total rho
527 14151 : IF (needs%drho) THEN
528 0 : DO idir = 1, 3
529 0 : NULLIFY (rho_set%drho(idir)%array)
530 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)))
531 : END DO
532 0 : rho_set%owns%drho = .TRUE.
533 0 : rho_set%has%drho = .FALSE.
534 : END IF
535 : ! Allocate the storage space for the components of the gradient for rhoa and rhob
536 121518 : IF (needs%drho_spin) THEN
537 31672 : DO idir = 1, 3
538 23754 : NULLIFY (rho_set%drhoa(idir)%array, rho_set%drhob(idir)%array)
539 118770 : ALLOCATE (rho_set%drhoa(idir)%array(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
540 126688 : ALLOCATE (rho_set%drhob(idir)%array(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
541 : END DO
542 7918 : rho_set%owns%drho_spin = .TRUE.
543 7918 : rho_set%has%drho_spin = .FALSE.
544 : END IF
545 : !
546 : END SELECT
547 :
548 : ! tau part
549 107367 : IF (needs%tau) THEN
550 996 : NULLIFY (rho_set%tau)
551 4980 : ALLOCATE (rho_set%tau(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
552 996 : rho_set%owns%tau = .TRUE.
553 : END IF
554 107367 : IF (needs%tau_spin) THEN
555 2 : NULLIFY (rho_set%tau_a, rho_set%tau_b)
556 10 : ALLOCATE (rho_set%tau_a(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
557 10 : ALLOCATE (rho_set%tau_b(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
558 2 : rho_set%owns%tau_spin = .TRUE.
559 2 : rho_set%has%tau_spin = .FALSE.
560 : END IF
561 :
562 : ! Laplace part
563 107367 : IF (needs%laplace_rho) THEN
564 0 : NULLIFY (rho_set%laplace_rho)
565 0 : ALLOCATE (rho_set%laplace_rho(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
566 0 : rho_set%owns%laplace_rho = .TRUE.
567 : END IF
568 107367 : IF (needs%laplace_rho_spin) THEN
569 0 : NULLIFY (rho_set%laplace_rhoa)
570 0 : NULLIFY (rho_set%laplace_rhob)
571 0 : ALLOCATE (rho_set%laplace_rhoa(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
572 0 : ALLOCATE (rho_set%laplace_rhob(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
573 0 : rho_set%owns%laplace_rho_spin = .TRUE.
574 0 : rho_set%has%laplace_rho_spin = .TRUE.
575 : END IF
576 :
577 107367 : END SUBROUTINE xc_rho_set_atom_update
578 :
579 : ! **************************************************************************************************
580 : !> \brief ...
581 : !> \param rho_set ...
582 : !> \param lsd ...
583 : !> \param nspins ...
584 : !> \param needs ...
585 : !> \param rho ...
586 : !> \param drho ...
587 : !> \param tau ...
588 : !> \param na ...
589 : !> \param ir ...
590 : ! **************************************************************************************************
591 3626980 : SUBROUTINE fill_rho_set(rho_set, lsd, nspins, needs, rho, drho, tau, na, ir)
592 :
593 : TYPE(xc_rho_set_type), INTENT(INOUT) :: rho_set
594 : LOGICAL, INTENT(IN) :: lsd
595 : INTEGER, INTENT(IN) :: nspins
596 : TYPE(xc_rho_cflags_type), INTENT(IN) :: needs
597 : REAL(dp), DIMENSION(:, :, :), INTENT(IN) :: rho
598 : REAL(dp), DIMENSION(:, :, :, :), INTENT(IN) :: drho
599 : REAL(dp), DIMENSION(:, :, :), INTENT(IN) :: tau
600 : INTEGER, INTENT(IN) :: na, ir
601 :
602 : REAL(KIND=dp), PARAMETER :: f13 = (1.0_dp/3.0_dp)
603 :
604 : INTEGER :: ia, idir, my_nspins
605 : LOGICAL :: gradient_f, tddft_split
606 :
607 3626980 : my_nspins = nspins
608 3626980 : tddft_split = .FALSE.
609 3626980 : IF (lsd .AND. nspins == 1) THEN
610 59400 : my_nspins = 2
611 59400 : tddft_split = .TRUE.
612 : END IF
613 :
614 : ! some checks
615 3626980 : IF (lsd) THEN
616 : ELSE
617 3002400 : CPASSERT(SIZE(rho, 3) == 1)
618 : END IF
619 3002400 : SELECT CASE (my_nspins)
620 : CASE (1)
621 3002400 : CPASSERT(.NOT. needs%rho_spin)
622 3002400 : CPASSERT(.NOT. needs%drho_spin)
623 3002400 : CPASSERT(.NOT. needs%norm_drho_spin)
624 3002400 : CPASSERT(.NOT. needs%rho_spin_1_3)
625 : CASE (2)
626 : CASE default
627 3626980 : CPABORT("Unsupported number of spins")
628 : END SELECT
629 :
630 : gradient_f = (needs%drho_spin .OR. needs%norm_drho_spin .OR. &
631 3626980 : needs%drho .OR. needs%norm_drho)
632 :
633 3002400 : SELECT CASE (my_nspins)
634 : CASE (1)
635 : ! Give rho to 1/3
636 3002400 : IF (needs%rho_1_3) THEN
637 6837000 : DO ia = 1, na
638 6837000 : rho_set%rho_1_3(ia, ir, 1) = MAX(rho(ia, ir, 1), 0.0_dp)**f13
639 : END DO
640 134200 : rho_set%owns%rho_1_3 = .TRUE.
641 134200 : rho_set%has%rho_1_3 = .TRUE.
642 : END IF
643 : ! Give the density
644 3002400 : IF (needs%rho) THEN
645 153302400 : DO ia = 1, na
646 153302400 : rho_set%rho(ia, ir, 1) = rho(ia, ir, 1)
647 : END DO
648 3002400 : rho_set%owns%rho = .TRUE.
649 3002400 : rho_set%has%rho = .TRUE.
650 : END IF
651 : ! Give the norm of the gradient of the density
652 3002400 : IF (needs%norm_drho) THEN
653 88695600 : DO ia = 1, na
654 88695600 : rho_set%norm_drho(ia, ir, 1) = drho(4, ia, ir, 1)
655 : END DO
656 1735600 : rho_set%owns%norm_drho = .TRUE.
657 1735600 : rho_set%has%norm_drho = .TRUE.
658 : END IF
659 : ! Give the three components of the gradient of the density
660 3002400 : IF (needs%drho) THEN
661 6942400 : DO idir = 1, 3
662 267822400 : DO ia = 1, na
663 266086800 : rho_set%drho(idir)%array(ia, ir, 1) = drho(idir, ia, ir, 1)
664 : END DO
665 : END DO
666 1735600 : rho_set%owns%drho = .TRUE.
667 1735600 : rho_set%has%drho = .TRUE.
668 : END IF
669 : CASE (2)
670 : ! Give the total density
671 624580 : IF (needs%rho) THEN
672 : ! this should never be the case unless you use LDA functionals with LSD
673 0 : IF (.NOT. tddft_split) THEN
674 0 : DO ia = 1, na
675 0 : rho_set%rho(ia, ir, 1) = rho(ia, ir, 1) + rho(ia, ir, 2)
676 : END DO
677 : ELSE
678 0 : DO ia = 1, na
679 0 : rho_set%rho(ia, ir, 1) = rho(ia, ir, 1)
680 : END DO
681 : END IF
682 0 : rho_set%owns%rho = .TRUE.
683 0 : rho_set%has%rho = .TRUE.
684 : END IF
685 : ! Give the total density to 1/3
686 624580 : IF (needs%rho_1_3) THEN
687 0 : IF (.NOT. tddft_split) THEN
688 0 : DO ia = 1, na
689 0 : rho_set%rho_1_3(ia, ir, 1) = MAX(rho(ia, ir, 1) + rho(ia, ir, 2), 0.0_dp)**f13
690 : END DO
691 : ELSE
692 0 : DO ia = 1, na
693 0 : rho_set%rho_1_3(ia, ir, 1) = MAX(rho(ia, ir, 1), 0.0_dp)**f13
694 : END DO
695 : END IF
696 0 : rho_set%owns%rho_1_3 = .TRUE.
697 0 : rho_set%has%rho_1_3 = .TRUE.
698 : END IF
699 : ! Give the spin densities to 1/3
700 624580 : IF (needs%rho_spin_1_3) THEN
701 75780 : IF (.NOT. tddft_split) THEN
702 3853260 : DO ia = 1, na
703 3777480 : rho_set%rhoa_1_3(ia, ir, 1) = MAX(rho(ia, ir, 1), 0.0_dp)**f13
704 3853260 : rho_set%rhob_1_3(ia, ir, 1) = MAX(rho(ia, ir, 2), 0.0_dp)**f13
705 : END DO
706 : ELSE
707 0 : DO ia = 1, na
708 0 : rho_set%rhoa_1_3(ia, ir, 1) = MAX(0.5_dp*rho(ia, ir, 1), 0.0_dp)**f13
709 0 : rho_set%rhob_1_3(ia, ir, 1) = rho_set%rhoa_1_3(ia, ir, 1)
710 : END DO
711 : END IF
712 75780 : rho_set%owns%rho_spin_1_3 = .TRUE.
713 75780 : rho_set%has%rho_spin_1_3 = .TRUE.
714 : END IF
715 : ! Give the spin densities rhoa and rhob
716 624580 : IF (needs%rho_spin) THEN
717 624580 : IF (.NOT. tddft_split) THEN
718 28812660 : DO ia = 1, na
719 28247480 : rho_set%rhoa(ia, ir, 1) = rho(ia, ir, 1)
720 28812660 : rho_set%rhob(ia, ir, 1) = rho(ia, ir, 2)
721 : END DO
722 : ELSE
723 3029400 : DO ia = 1, na
724 2970000 : rho_set%rhoa(ia, ir, 1) = 0.5_dp*rho(ia, ir, 1)
725 3029400 : rho_set%rhob(ia, ir, 1) = rho_set%rhoa(ia, ir, 1)
726 : END DO
727 : END IF
728 624580 : rho_set%owns%rho_spin = .TRUE.
729 624580 : rho_set%has%rho_spin = .TRUE.
730 : END IF
731 : ! Give the norm of the gradient of the total density
732 624580 : IF (needs%norm_drho) THEN
733 318380 : IF (.NOT. tddft_split) THEN
734 13992060 : DO ia = 1, na
735 : rho_set%norm_drho(ia, ir, 1) = SQRT( &
736 : (drho(1, ia, ir, 1) + drho(1, ia, ir, 2))**2 + &
737 : (drho(2, ia, ir, 1) + drho(2, ia, ir, 2))**2 + &
738 13992060 : (drho(3, ia, ir, 1) + drho(3, ia, ir, 2))**2)
739 : END DO
740 : ELSE
741 2233800 : DO ia = 1, na
742 2233800 : rho_set%norm_drho(ia, ir, 1) = drho(4, ia, ir, 1)
743 : END DO
744 : END IF
745 318380 : rho_set%owns%norm_drho = .TRUE.
746 318380 : rho_set%has%norm_drho = .TRUE.
747 : END IF
748 : ! Give the norm of the gradient of rhoa and of rhob separatedly
749 624580 : IF (needs%norm_drho_spin) THEN
750 320780 : IF (.NOT. tddft_split) THEN
751 14114460 : DO ia = 1, na
752 13837480 : rho_set%norm_drhoa(ia, ir, 1) = drho(4, ia, ir, 1)
753 14114460 : rho_set%norm_drhob(ia, ir, 1) = drho(4, ia, ir, 2)
754 : END DO
755 : ELSE
756 2233800 : DO ia = 1, na
757 2190000 : rho_set%norm_drhoa(ia, ir, 1) = 0.5_dp*drho(4, ia, ir, 1)
758 2233800 : rho_set%norm_drhob(ia, ir, 1) = rho_set%norm_drhoa(ia, ir, 1)
759 : END DO
760 : END IF
761 320780 : rho_set%owns%norm_drho_spin = .TRUE.
762 320780 : rho_set%has%norm_drho_spin = .TRUE.
763 : END IF
764 : ! Give the components of the gradient for the total rho
765 624580 : IF (needs%drho) THEN
766 0 : IF (.NOT. tddft_split) THEN
767 0 : DO idir = 1, 3
768 0 : DO ia = 1, na
769 0 : rho_set%drho(idir)%array(ia, ir, 1) = drho(idir, ia, ir, 1) + drho(idir, ia, ir, 2)
770 : END DO
771 : END DO
772 : ELSE
773 0 : DO idir = 1, 3
774 0 : DO ia = 1, na
775 0 : rho_set%drho(idir)%array(ia, ir, 1) = drho(idir, ia, ir, 1)
776 : END DO
777 : END DO
778 : END IF
779 0 : rho_set%owns%drho = .TRUE.
780 0 : rho_set%has%drho = .TRUE.
781 : END IF
782 : ! Give the components of the gradient for rhoa and rhob
783 4251560 : IF (needs%drho_spin) THEN
784 322280 : IF (.NOT. tddft_split) THEN
785 1113920 : DO idir = 1, 3
786 42851360 : DO ia = 1, na
787 41737440 : rho_set%drhoa(idir)%array(ia, ir, 1) = drho(idir, ia, ir, 1)
788 42572880 : rho_set%drhob(idir)%array(ia, ir, 1) = drho(idir, ia, ir, 2)
789 : END DO
790 : END DO
791 : ELSE
792 175200 : DO idir = 1, 3
793 6745200 : DO ia = 1, na
794 6570000 : rho_set%drhoa(idir)%array(ia, ir, 1) = 0.5_dp*drho(idir, ia, ir, 1)
795 6701400 : rho_set%drhob(idir)%array(ia, ir, 1) = rho_set%drhoa(idir)%array(ia, ir, 1)
796 : END DO
797 : END DO
798 : END IF
799 322280 : rho_set%owns%drho_spin = .TRUE.
800 322280 : rho_set%has%drho_spin = .TRUE.
801 : END IF
802 : !
803 : END SELECT
804 :
805 : ! tau part
806 3626980 : IF (needs%tau .OR. needs%tau_spin) THEN
807 3626980 : CPASSERT(SIZE(tau, 3) == my_nspins)
808 : END IF
809 3626980 : IF (needs%tau) THEN
810 33400 : IF (my_nspins == 2) THEN
811 0 : DO ia = 1, na
812 0 : rho_set%tau(ia, ir, 1) = tau(ia, ir, 1) + tau(ia, ir, 2)
813 : END DO
814 0 : rho_set%owns%tau = .TRUE.
815 0 : rho_set%has%tau = .TRUE.
816 : ELSE
817 1890600 : DO ia = 1, na
818 1890600 : rho_set%tau(ia, ir, 1) = tau(ia, ir, 1)
819 : END DO
820 33400 : rho_set%owns%tau = .TRUE.
821 33400 : rho_set%has%tau = .TRUE.
822 : END IF
823 : END IF
824 3626980 : IF (needs%tau_spin) THEN
825 0 : DO ia = 1, na
826 0 : rho_set%tau_a(ia, ir, 1) = tau(ia, ir, 1)
827 0 : rho_set%tau_b(ia, ir, 1) = tau(ia, ir, 2)
828 : END DO
829 0 : rho_set%owns%tau_spin = .TRUE.
830 0 : rho_set%has%tau_spin = .TRUE.
831 : END IF
832 :
833 3626980 : END SUBROUTINE fill_rho_set
834 :
835 : END MODULE xc_atom
|
