(git:badb799)
Loading...
Searching...
No Matches
xc_atom.F
Go to the documentation of this file.
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! **************************************************************************************************
9MODULE xc_atom
10
15 USE kinds, ONLY: dp
17 USE pw_types, ONLY: pw_r3d_rs_type
18 USE xc, ONLY: divide_by_norm_drho,&
20 USE xc_derivative_desc, ONLY: &
31#include "../base/base_uses.f90"
32
33 IMPLICIT NONE
34
35 PRIVATE
36
37 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_atom'
38
40
41CONTAINS
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 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(:, :), POINTER :: 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 REAL(dp), DIMENSION(:, :, :), POINTER :: deriv_data
96 REAL(kind=dp) :: drho_cutoff
97 TYPE(xc_derivative_type), POINTER :: deriv_att
98
99 CALL timeset(routinen, handle)
100 my_only_energy = .false.
101 IF (PRESENT(energy_only)) my_only_energy = energy_only
102
103 IF (PRESENT(adiabatic_rescale_factor)) THEN
104 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 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 deriv_order=deriv_order)
118
119 CALL xc_rho_set_get(rho_set, drho_cutoff=drho_cutoff)
120
121 NULLIFY (deriv_data)
122
123 ! EXC energy
124 deriv_att => xc_dset_get_derivative(deriv_set, [INTEGER::])
125 exc = 0.0_dp
126 IF (ASSOCIATED(deriv_att)) THEN
127 CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
128 DO ir = 1, nr
129 DO ia = 1, na
130 exc = exc + deriv_data(ia, ir, 1)*w(ia, ir)
131 END DO
132 END DO
133 NULLIFY (deriv_data)
134 END IF
135 ! Calculate the potential only if needed
136 IF (.NOT. my_only_energy) THEN
137 ! Derivative with respect to the density
138 IF (lsd) THEN
139 deriv_att => xc_dset_get_derivative(deriv_set, [deriv_rhoa])
140 IF (ASSOCIATED(deriv_att)) THEN
141 CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
142 vxc(:, :, 1) = deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
143 NULLIFY (deriv_data)
144 END IF
145 deriv_att => xc_dset_get_derivative(deriv_set, [deriv_rhob])
146 IF (ASSOCIATED(deriv_att)) THEN
147 CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
148 vxc(:, :, 2) = deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
149 NULLIFY (deriv_data)
150 END IF
151 deriv_att => xc_dset_get_derivative(deriv_set, [deriv_rho])
152 IF (ASSOCIATED(deriv_att)) THEN
153 CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
154 vxc(:, :, 1) = vxc(:, :, 1) + deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
155 vxc(:, :, 2) = vxc(:, :, 2) + deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
156 NULLIFY (deriv_data)
157 END IF
158 ELSE
159 deriv_att => xc_dset_get_derivative(deriv_set, [deriv_rho])
160 IF (ASSOCIATED(deriv_att)) THEN
161 CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
162 vxc(:, :, 1) = deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
163 NULLIFY (deriv_data)
164 END IF
165 END IF ! lsd
166
167 ! Derivatives with respect to the gradient
168 IF (lsd) THEN
169 deriv_att => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa])
170 IF (ASSOCIATED(deriv_att)) THEN
171 CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
172 DO ir = 1, nr
173 DO ia = 1, na
174 DO idir = 1, 3
175 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 rho_set%norm_drhoa(ia, ir, 1)*my_adiabatic_rescale_factor
179 ELSE
180 vxg(idir, ia, ir, 1) = 0.0_dp
181 END IF
182 END DO
183 END DO
184 END DO
185 NULLIFY (deriv_data)
186 END IF
187 deriv_att => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob])
188 IF (ASSOCIATED(deriv_att)) THEN
189 CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
190 DO ir = 1, nr
191 DO ia = 1, na
192 DO idir = 1, 3
193 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 rho_set%norm_drhob(ia, ir, 1)*my_adiabatic_rescale_factor
197 ELSE
198 vxg(idir, ia, ir, 2) = 0.0_dp
199 END IF
200 END DO
201 END DO
202 END DO
203 NULLIFY (deriv_data)
204 END IF
205 ! Cross Terms
206 deriv_att => xc_dset_get_derivative(deriv_set, [deriv_norm_drho])
207 IF (ASSOCIATED(deriv_att)) THEN
208 CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
209 DO ir = 1, nr
210 DO ia = 1, na
211 DO idir = 1, 3
212 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 my_adiabatic_rescale_factor
219 END IF
220 END DO
221 END DO
222 END DO
223 NULLIFY (deriv_data)
224 END IF
225 ELSE
226 deriv_att => xc_dset_get_derivative(deriv_set, [deriv_norm_drho])
227 IF (ASSOCIATED(deriv_att)) THEN
228 CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
229 DO ir = 1, nr
230 DO ia = 1, na
231 IF (rho_set%norm_drho(ia, ir, 1) > drho_cutoff) THEN
232 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 rho_set%norm_drho(ia, ir, 1)*my_adiabatic_rescale_factor
236 END DO
237 ELSE
238 vxg(1:3, ia, ir, 1) = 0.0_dp
239 END IF
240 END DO
241 END DO
242 NULLIFY (deriv_data)
243 END IF
244 END IF ! lsd
245 ! Derivative with respect to tau
246 IF (lsd) THEN
247 deriv_att => xc_dset_get_derivative(deriv_set, [deriv_tau_a])
248 IF (ASSOCIATED(deriv_att)) THEN
249 CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
250 vtau(:, :, 1) = deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
251 NULLIFY (deriv_data)
252 END IF
253 deriv_att => xc_dset_get_derivative(deriv_set, [deriv_tau_b])
254 IF (ASSOCIATED(deriv_att)) THEN
255 CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
256 vtau(:, :, 2) = deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
257 NULLIFY (deriv_data)
258 END IF
259 deriv_att => xc_dset_get_derivative(deriv_set, [deriv_tau])
260 IF (ASSOCIATED(deriv_att)) THEN
261 CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
262 vtau(:, :, 1) = vtau(:, :, 1) + deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
263 vtau(:, :, 2) = vtau(:, :, 2) + deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
264 NULLIFY (deriv_data)
265 END IF
266 ELSE
267 deriv_att => xc_dset_get_derivative(deriv_set, [deriv_tau])
268 IF (ASSOCIATED(deriv_att)) THEN
269 CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
270 vtau(:, :, 1) = deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
271 NULLIFY (deriv_data)
272 END IF
273 END IF ! lsd
274 END IF ! only_energy
275
276 CALL timestop(handle)
277
278 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 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(:, :), POINTER :: 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 REAL(dp), DIMENSION(:, :, :), POINTER :: deriv_data
316 REAL(kind=dp) :: drho_cutoff
317 TYPE(xc_derivative_type), POINTER :: deriv_att
318
319 CALL timeset(routinen, handle)
320
321 mark_used(vxc)
322 mark_used(vtau)
323
324 my_adiabatic_rescale_factor = 1.0_dp
325 my_deriv_order = 2
326
327 gradient_f = (needs%drho_spin .OR. needs%norm_drho_spin .OR. &
328 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 deriv_order=my_deriv_order)
336
337 CALL xc_rho_set_get(rho_set, drho_cutoff=drho_cutoff)
338
339 NULLIFY (deriv_data)
340
341 ! nabla v_xc (using the vxg arrays)
342 ! there's no point doing this when lsd = false
343 IF (lsd) THEN
344 deriv_att => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa])
345 IF (ASSOCIATED(deriv_att)) THEN
346 CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
347 DO ir = 1, nr
348 DO ia = 1, na
349 DO idir = 1, 3
350 vxg(idir, ia, ir, 1) = rho_set%drhoa(idir)%array(ia, ir, 1)* &
351 deriv_data(ia, ir, 1)
352 END DO !idir
353 END DO !ia
354 END DO !ir
355 NULLIFY (deriv_data)
356 END IF
357 deriv_att => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob])
358 IF (ASSOCIATED(deriv_att)) THEN
359 CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
360 DO ir = 1, nr
361 DO ia = 1, na
362 DO idir = 1, 3
363 vxg(idir, ia, ir, 2) = rho_set%drhob(idir)%array(ia, ir, 1)* &
364 deriv_data(ia, ir, 1)
365 END DO !idir
366 END DO !ia
367 END DO !ir
368 NULLIFY (deriv_data)
369 END IF
370 END IF
371 ! EXC energy ! is that needed for epr?
372 deriv_att => xc_dset_get_derivative(deriv_set, [INTEGER::])
373 exc = 0.0_dp
374 IF (ASSOCIATED(deriv_att)) THEN
375 CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
376 DO ir = 1, nr
377 DO ia = 1, na
378 exc = exc + deriv_data(ia, ir, 1)*w(ia, ir)
379 END DO
380 END DO
381 NULLIFY (deriv_data)
382 END IF
383
384 CALL timestop(handle)
385
386 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 do_triplet ...
398!> \param do_sf ...
399! **************************************************************************************************
400 SUBROUTINE xc_2nd_deriv_of_r(rho_set, rho1_set, xc_section, &
401 deriv_set, w, vxc, vxg, do_triplet, do_sf)
402
403! As input of this routine one gets rho and drho on a one dimensional grid.
404! The grid is the angular grid corresponding to a given point ir on the radial grid.
405! The derivatives are calculated on this one dimensional grid, the results are stored in
406! vxc(1:na,ir,ispin), vxg(1:na,ir,ispin), vxg_cross(1:na,ir,ispin)
407! Afterwords the arrays containing the derivatives are put to zero so that the routine
408! can safely be called for the next radial point ir
409
410 TYPE(xc_rho_set_type), INTENT(IN) :: rho_set, rho1_set
411 TYPE(section_vals_type), POINTER :: xc_section
412 TYPE(xc_derivative_set_type), INTENT(INOUT) :: deriv_set
413 REAL(dp), DIMENSION(:, :), POINTER :: w
414 REAL(dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: vxc
415 REAL(dp), DIMENSION(:, :, :, :), POINTER :: vxg
416 LOGICAL, INTENT(IN), OPTIONAL :: do_triplet, do_sf
417
418 CHARACTER(LEN=*), PARAMETER :: routinen = 'xc_2nd_deriv_of_r'
419
420 INTEGER :: handle, ispin, nspins
421 LOGICAL :: lsd, my_do_sf
422 REAL(dp) :: drho_cutoff, my_fac_triplet
423 TYPE(cp_sll_xc_deriv_type), POINTER :: pos
424 TYPE(pw_pool_type), POINTER :: pw_pool
425 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: vxc_pw, vxc_tau_pw
426 TYPE(section_vals_type), POINTER :: xc_fun_section
427 TYPE(xc_derivative_type), POINTER :: deriv_att
428
429 CALL timeset(routinen, handle)
430
431 nspins = SIZE(vxc, 3)
432 lsd = (nspins == 2)
433 IF (ASSOCIATED(rho_set%rhoa)) THEN
434 lsd = .true.
435 END IF
436 my_fac_triplet = 1.0_dp
437 IF (PRESENT(do_triplet)) THEN
438 IF (do_triplet) my_fac_triplet = -1.0_dp
439 END IF
440
441 my_do_sf = .false.
442 IF (PRESENT(do_sf)) my_do_sf = do_sf
443
444 CALL xc_rho_set_get(rho_set, drho_cutoff=drho_cutoff)
445 xc_fun_section => section_vals_get_subs_vals(xc_section, &
446 "XC_FUNCTIONAL")
447
448 ! Calculate the derivatives
449 CALL xc_functionals_eval(xc_fun_section, &
450 lsd=lsd, &
451 rho_set=rho_set, &
452 deriv_set=deriv_set, &
453 deriv_order=2)
454
455 CALL divide_by_norm_drho(deriv_set, rho_set, lsd)
456
457 ! multiply by w
458 pos => deriv_set%derivs
459 DO WHILE (cp_sll_xc_deriv_next(pos, el_att=deriv_att))
460 deriv_att%deriv_data(:, :, 1) = w(:, :)*deriv_att%deriv_data(:, :, 1)
461 END DO
462
463 NULLIFY (pw_pool)
464 ALLOCATE (vxc_pw(nspins))
465 DO ispin = 1, nspins
466 vxc_pw(ispin)%array => vxc(:, :, ispin:ispin)
467 END DO
468
469 NULLIFY (vxc_tau_pw)
470
471 CALL xc_calc_2nd_deriv_analytical(vxc_pw, vxc_tau_pw, deriv_set, rho_set, rho1_set, pw_pool, &
472 xc_section, gapw=.true., vxg=vxg, tddfpt_fac=my_fac_triplet, spinflip=do_sf)
473
474 DEALLOCATE (vxc_pw)
475
476 ! zero the derivative data for the next call
477 pos => deriv_set%derivs
478 DO WHILE (cp_sll_xc_deriv_next(pos, el_att=deriv_att))
479 deriv_att%deriv_data = 0.0_dp
480 END DO
481
482 CALL timestop(handle)
483
484 END SUBROUTINE xc_2nd_deriv_of_r
485
486! **************************************************************************************************
487!> \brief ...
488!> \param rho_set ...
489!> \param needs ...
490!> \param nspins ...
491!> \param bo ...
492! **************************************************************************************************
493 SUBROUTINE xc_rho_set_atom_update(rho_set, needs, nspins, bo)
494
495! This routine allocates the storage arrays for rho and drho
496! In calculate_vxc_atom this is called once for each atomic_kind,
497! After the loop over all the atoms of the kind and over all the points
498! of the radial grid for each atom, rho_set is deallocated.
499! Within the same kind, at each new point on the radial grid, the rho_set
500! arrays rho and drho are overwritten.
501
502 TYPE(xc_rho_set_type), INTENT(INOUT) :: rho_set
503 TYPE(xc_rho_cflags_type), INTENT(IN) :: needs
504 INTEGER, INTENT(IN) :: nspins
505 INTEGER, DIMENSION(2, 3), INTENT(IN) :: bo
506
507 INTEGER :: idir
508
509 SELECT CASE (nspins)
510 CASE (1)
511! What is this for?
512 IF (needs%rho_1_3) THEN
513 NULLIFY (rho_set%rho_1_3)
514 ALLOCATE (rho_set%rho_1_3(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
515 rho_set%owns%rho_1_3 = .true.
516 rho_set%has%rho_1_3 = .false.
517 END IF
518! Allocate the storage space for the density
519 IF (needs%rho) THEN
520 NULLIFY (rho_set%rho)
521 ALLOCATE (rho_set%rho(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
522 rho_set%owns%rho = .true.
523 rho_set%has%rho = .false.
524 END IF
525! Allocate the storage space for the norm of the gradient of the density
526 IF (needs%norm_drho) THEN
527 NULLIFY (rho_set%norm_drho)
528 ALLOCATE (rho_set%norm_drho(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
529 rho_set%owns%norm_drho = .true.
530 rho_set%has%norm_drho = .false.
531 END IF
532! Allocate the storage space for the three components of the gradient of the density
533 IF (needs%drho) THEN
534 DO idir = 1, 3
535 NULLIFY (rho_set%drho(idir)%array)
536 ALLOCATE (rho_set%drho(idir)%array(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
537 END DO
538 rho_set%owns%drho = .true.
539 rho_set%has%drho = .false.
540 END IF
541 CASE (2)
542! Allocate the storage space for the total density
543 IF (needs%rho) THEN
544 ! this should never be the case unless you use LDA functionals with LSD
545 NULLIFY (rho_set%rho)
546 ALLOCATE (rho_set%rho(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
547 rho_set%owns%rho = .true.
548 rho_set%has%rho = .false.
549 END IF
550! What is this for?
551 IF (needs%rho_1_3) THEN
552 NULLIFY (rho_set%rho_1_3)
553 ALLOCATE (rho_set%rho_1_3(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
554 rho_set%owns%rho_1_3 = .true.
555 rho_set%has%rho_1_3 = .false.
556 END IF
557! What is this for?
558 IF (needs%rho_spin_1_3) THEN
559 NULLIFY (rho_set%rhoa_1_3, rho_set%rhob_1_3)
560 ALLOCATE (rho_set%rhoa_1_3(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
561 ALLOCATE (rho_set%rhob_1_3(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
562 rho_set%owns%rho_spin_1_3 = .true.
563 rho_set%has%rho_spin_1_3 = .false.
564 END IF
565! Allocate the storage space for the spin densities rhoa and rhob
566 IF (needs%rho_spin) THEN
567 NULLIFY (rho_set%rhoa, rho_set%rhob)
568 ALLOCATE (rho_set%rhoa(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
569 ALLOCATE (rho_set%rhob(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
570 rho_set%owns%rho_spin = .true.
571 rho_set%has%rho_spin = .false.
572 END IF
573! Allocate the storage space for the norm of the gradient of the total density
574 IF (needs%norm_drho) THEN
575 NULLIFY (rho_set%norm_drho)
576 ALLOCATE (rho_set%norm_drho(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
577 rho_set%owns%norm_drho = .true.
578 rho_set%has%norm_drho = .false.
579 END IF
580! Allocate the storage space for the norm of the gradient of rhoa and of rhob separatedly
581 IF (needs%norm_drho_spin) THEN
582 NULLIFY (rho_set%norm_drhoa, rho_set%norm_drhob)
583 ALLOCATE (rho_set%norm_drhoa(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
584 ALLOCATE (rho_set%norm_drhob(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
585 rho_set%owns%norm_drho_spin = .true.
586 rho_set%has%norm_drho_spin = .false.
587 END IF
588! Allocate the storage space for the components of the gradient for the total rho
589 IF (needs%drho) THEN
590 DO idir = 1, 3
591 NULLIFY (rho_set%drho(idir)%array)
592 ALLOCATE (rho_set%drho(idir)%array(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
593 END DO
594 rho_set%owns%drho = .true.
595 rho_set%has%drho = .false.
596 END IF
597! Allocate the storage space for the components of the gradient for rhoa and rhob
598 IF (needs%drho_spin) THEN
599 DO idir = 1, 3
600 NULLIFY (rho_set%drhoa(idir)%array, rho_set%drhob(idir)%array)
601 ALLOCATE (rho_set%drhoa(idir)%array(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
602 ALLOCATE (rho_set%drhob(idir)%array(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
603 END DO
604 rho_set%owns%drho_spin = .true.
605 rho_set%has%drho_spin = .false.
606 END IF
607!
608 END SELECT
609
610 ! tau part
611 IF (needs%tau) THEN
612 NULLIFY (rho_set%tau)
613 ALLOCATE (rho_set%tau(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
614 rho_set%owns%tau = .true.
615 END IF
616 IF (needs%tau_spin) THEN
617 NULLIFY (rho_set%tau_a, rho_set%tau_b)
618 ALLOCATE (rho_set%tau_a(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
619 ALLOCATE (rho_set%tau_b(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
620 rho_set%owns%tau_spin = .true.
621 rho_set%has%tau_spin = .false.
622 END IF
623
624 ! Laplace part
625 IF (needs%laplace_rho) THEN
626 NULLIFY (rho_set%laplace_rho)
627 ALLOCATE (rho_set%laplace_rho(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
628 rho_set%owns%laplace_rho = .true.
629 END IF
630 IF (needs%laplace_rho_spin) THEN
631 NULLIFY (rho_set%laplace_rhoa)
632 NULLIFY (rho_set%laplace_rhob)
633 ALLOCATE (rho_set%laplace_rhoa(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
634 ALLOCATE (rho_set%laplace_rhob(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
635 rho_set%owns%laplace_rho_spin = .true.
636 rho_set%has%laplace_rho_spin = .true.
637 END IF
638
639 END SUBROUTINE xc_rho_set_atom_update
640
641! **************************************************************************************************
642!> \brief ...
643!> \param rho_set ...
644!> \param lsd ...
645!> \param nspins ...
646!> \param needs ...
647!> \param rho ...
648!> \param drho ...
649!> \param tau ...
650!> \param na ...
651!> \param ir ...
652! **************************************************************************************************
653 SUBROUTINE fill_rho_set(rho_set, lsd, nspins, needs, rho, drho, tau, na, ir)
654
655 TYPE(xc_rho_set_type), INTENT(INOUT) :: rho_set
656 LOGICAL, INTENT(IN) :: lsd
657 INTEGER, INTENT(IN) :: nspins
658 TYPE(xc_rho_cflags_type), INTENT(IN) :: needs
659 REAL(dp), DIMENSION(:, :, :), INTENT(IN) :: rho
660 REAL(dp), DIMENSION(:, :, :, :), INTENT(IN) :: drho
661 REAL(dp), DIMENSION(:, :, :), INTENT(IN) :: tau
662 INTEGER, INTENT(IN) :: na, ir
663
664 REAL(kind=dp), PARAMETER :: f13 = (1.0_dp/3.0_dp)
665
666 INTEGER :: ia, idir, my_nspins
667 LOGICAL :: gradient_f, tddft_split
668
669 my_nspins = nspins
670 tddft_split = .false.
671 IF (lsd .AND. nspins == 1) THEN
672 my_nspins = 2
673 tddft_split = .true.
674 END IF
675
676 ! some checks
677 IF (lsd) THEN
678 ELSE
679 cpassert(SIZE(rho, 3) == 1)
680 END IF
681 SELECT CASE (my_nspins)
682 CASE (1)
683 cpassert(.NOT. needs%rho_spin)
684 cpassert(.NOT. needs%drho_spin)
685 cpassert(.NOT. needs%norm_drho_spin)
686 cpassert(.NOT. needs%rho_spin_1_3)
687 CASE (2)
688 CASE default
689 cpabort("Unsupported number of spins")
690 END SELECT
691
692 gradient_f = (needs%drho_spin .OR. needs%norm_drho_spin .OR. &
693 needs%drho .OR. needs%norm_drho)
694
695 SELECT CASE (my_nspins)
696 CASE (1)
697 ! Give rho to 1/3
698 IF (needs%rho_1_3) THEN
699 DO ia = 1, na
700 rho_set%rho_1_3(ia, ir, 1) = max(rho(ia, ir, 1), 0.0_dp)**f13
701 END DO
702 rho_set%owns%rho_1_3 = .true.
703 rho_set%has%rho_1_3 = .true.
704 END IF
705 ! Give the density
706 IF (needs%rho) THEN
707 DO ia = 1, na
708 rho_set%rho(ia, ir, 1) = rho(ia, ir, 1)
709 END DO
710 rho_set%owns%rho = .true.
711 rho_set%has%rho = .true.
712 END IF
713 ! Give the norm of the gradient of the density
714 IF (needs%norm_drho) THEN
715 DO ia = 1, na
716 rho_set%norm_drho(ia, ir, 1) = drho(4, ia, ir, 1)
717 END DO
718 rho_set%owns%norm_drho = .true.
719 rho_set%has%norm_drho = .true.
720 END IF
721 ! Give the three components of the gradient of the density
722 IF (needs%drho) THEN
723 DO idir = 1, 3
724 DO ia = 1, na
725 rho_set%drho(idir)%array(ia, ir, 1) = drho(idir, ia, ir, 1)
726 END DO
727 END DO
728 rho_set%owns%drho = .true.
729 rho_set%has%drho = .true.
730 END IF
731 CASE (2)
732 ! Give the total density
733 IF (needs%rho) THEN
734 ! this should never be the case unless you use LDA functionals with LSD
735 IF (.NOT. tddft_split) THEN
736 DO ia = 1, na
737 rho_set%rho(ia, ir, 1) = rho(ia, ir, 1) + rho(ia, ir, 2)
738 END DO
739 ELSE
740 DO ia = 1, na
741 rho_set%rho(ia, ir, 1) = rho(ia, ir, 1)
742 END DO
743 END IF
744 rho_set%owns%rho = .true.
745 rho_set%has%rho = .true.
746 END IF
747 ! Give the total density to 1/3
748 IF (needs%rho_1_3) THEN
749 IF (.NOT. tddft_split) THEN
750 DO ia = 1, na
751 rho_set%rho_1_3(ia, ir, 1) = max(rho(ia, ir, 1) + rho(ia, ir, 2), 0.0_dp)**f13
752 END DO
753 ELSE
754 DO ia = 1, na
755 rho_set%rho_1_3(ia, ir, 1) = max(rho(ia, ir, 1), 0.0_dp)**f13
756 END DO
757 END IF
758 rho_set%owns%rho_1_3 = .true.
759 rho_set%has%rho_1_3 = .true.
760 END IF
761 ! Give the spin densities to 1/3
762 IF (needs%rho_spin_1_3) THEN
763 IF (.NOT. tddft_split) THEN
764 DO ia = 1, na
765 rho_set%rhoa_1_3(ia, ir, 1) = max(rho(ia, ir, 1), 0.0_dp)**f13
766 rho_set%rhob_1_3(ia, ir, 1) = max(rho(ia, ir, 2), 0.0_dp)**f13
767 END DO
768 ELSE
769 DO ia = 1, na
770 rho_set%rhoa_1_3(ia, ir, 1) = max(0.5_dp*rho(ia, ir, 1), 0.0_dp)**f13
771 rho_set%rhob_1_3(ia, ir, 1) = rho_set%rhoa_1_3(ia, ir, 1)
772 END DO
773 END IF
774 rho_set%owns%rho_spin_1_3 = .true.
775 rho_set%has%rho_spin_1_3 = .true.
776 END IF
777 ! Give the spin densities rhoa and rhob
778 IF (needs%rho_spin) THEN
779 IF (.NOT. tddft_split) THEN
780 DO ia = 1, na
781 rho_set%rhoa(ia, ir, 1) = rho(ia, ir, 1)
782 rho_set%rhob(ia, ir, 1) = rho(ia, ir, 2)
783 END DO
784 ELSE
785 DO ia = 1, na
786 rho_set%rhoa(ia, ir, 1) = 0.5_dp*rho(ia, ir, 1)
787 rho_set%rhob(ia, ir, 1) = rho_set%rhoa(ia, ir, 1)
788 END DO
789 END IF
790 rho_set%owns%rho_spin = .true.
791 rho_set%has%rho_spin = .true.
792 END IF
793 ! Give the norm of the gradient of the total density
794 IF (needs%norm_drho) THEN
795 IF (.NOT. tddft_split) THEN
796 DO ia = 1, na
797 rho_set%norm_drho(ia, ir, 1) = sqrt( &
798 (drho(1, ia, ir, 1) + drho(1, ia, ir, 2))**2 + &
799 (drho(2, ia, ir, 1) + drho(2, ia, ir, 2))**2 + &
800 (drho(3, ia, ir, 1) + drho(3, ia, ir, 2))**2)
801 END DO
802 ELSE
803 DO ia = 1, na
804 rho_set%norm_drho(ia, ir, 1) = drho(4, ia, ir, 1)
805 END DO
806 END IF
807 rho_set%owns%norm_drho = .true.
808 rho_set%has%norm_drho = .true.
809 END IF
810 ! Give the norm of the gradient of rhoa and of rhob separatedly
811 IF (needs%norm_drho_spin) THEN
812 IF (.NOT. tddft_split) THEN
813 DO ia = 1, na
814 rho_set%norm_drhoa(ia, ir, 1) = drho(4, ia, ir, 1)
815 rho_set%norm_drhob(ia, ir, 1) = drho(4, ia, ir, 2)
816 END DO
817 ELSE
818 DO ia = 1, na
819 rho_set%norm_drhoa(ia, ir, 1) = 0.5_dp*drho(4, ia, ir, 1)
820 rho_set%norm_drhob(ia, ir, 1) = rho_set%norm_drhoa(ia, ir, 1)
821 END DO
822 END IF
823 rho_set%owns%norm_drho_spin = .true.
824 rho_set%has%norm_drho_spin = .true.
825 END IF
826 ! Give the components of the gradient for the total rho
827 IF (needs%drho) THEN
828 IF (.NOT. tddft_split) THEN
829 DO idir = 1, 3
830 DO ia = 1, na
831 rho_set%drho(idir)%array(ia, ir, 1) = drho(idir, ia, ir, 1) + drho(idir, ia, ir, 2)
832 END DO
833 END DO
834 ELSE
835 DO idir = 1, 3
836 DO ia = 1, na
837 rho_set%drho(idir)%array(ia, ir, 1) = drho(idir, ia, ir, 1)
838 END DO
839 END DO
840 END IF
841 rho_set%owns%drho = .true.
842 rho_set%has%drho = .true.
843 END IF
844 ! Give the components of the gradient for rhoa and rhob
845 IF (needs%drho_spin) THEN
846 IF (.NOT. tddft_split) THEN
847 DO idir = 1, 3
848 DO ia = 1, na
849 rho_set%drhoa(idir)%array(ia, ir, 1) = drho(idir, ia, ir, 1)
850 rho_set%drhob(idir)%array(ia, ir, 1) = drho(idir, ia, ir, 2)
851 END DO
852 END DO
853 ELSE
854 DO idir = 1, 3
855 DO ia = 1, na
856 rho_set%drhoa(idir)%array(ia, ir, 1) = 0.5_dp*drho(idir, ia, ir, 1)
857 rho_set%drhob(idir)%array(ia, ir, 1) = rho_set%drhoa(idir)%array(ia, ir, 1)
858 END DO
859 END DO
860 END IF
861 rho_set%owns%drho_spin = .true.
862 rho_set%has%drho_spin = .true.
863 END IF
864 !
865 END SELECT
866
867 ! tau part
868 IF (needs%tau .OR. needs%tau_spin) THEN
869 cpassert(SIZE(tau, 3) == my_nspins)
870 END IF
871 IF (needs%tau) THEN
872 IF (my_nspins == 2) THEN
873 DO ia = 1, na
874 rho_set%tau(ia, ir, 1) = tau(ia, ir, 1) + tau(ia, ir, 2)
875 END DO
876 rho_set%owns%tau = .true.
877 rho_set%has%tau = .true.
878 ELSE
879 DO ia = 1, na
880 rho_set%tau(ia, ir, 1) = tau(ia, ir, 1)
881 END DO
882 rho_set%owns%tau = .true.
883 rho_set%has%tau = .true.
884 END IF
885 END IF
886 IF (needs%tau_spin) THEN
887 DO ia = 1, na
888 rho_set%tau_a(ia, ir, 1) = tau(ia, ir, 1)
889 rho_set%tau_b(ia, ir, 1) = tau(ia, ir, 2)
890 END DO
891 rho_set%owns%tau_spin = .true.
892 rho_set%has%tau_spin = .true.
893 END IF
894
895 END SUBROUTINE fill_rho_set
896
897END MODULE xc_atom
logical function, public cp_sll_xc_deriv_next(iterator, el_att)
returns true if the actual element is valid (i.e. iterator ont at end) moves the iterator to the next...
objects that represent the structure of input sections and the data contained in an input section
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
subroutine, public vxc_of_r_epr(xc_fun_section, rho_set, deriv_set, needs, w, lsd, na, nr, exc, vxc, vxg, vtau)
Specific EPR version of vxc_of_r_new.
Definition xc_atom.F:297
subroutine, public vxc_of_r_new(xc_fun_section, rho_set, deriv_set, deriv_order, needs, w, lsd, na, nr, exc, vxc, vxg, vtau, energy_only, adiabatic_rescale_factor)
...
Definition xc_atom.F:64
subroutine, public xc_rho_set_atom_update(rho_set, needs, nspins, bo)
...
Definition xc_atom.F:494
subroutine, public fill_rho_set(rho_set, lsd, nspins, needs, rho, drho, tau, na, ir)
...
Definition xc_atom.F:654
subroutine, public xc_2nd_deriv_of_r(rho_set, rho1_set, xc_section, deriv_set, w, vxc, vxg, do_triplet, do_sf)
...
Definition xc_atom.F:402
Module with functions to handle derivative descriptors. derivative description are strings have the f...
integer, parameter, public deriv_norm_drho
integer, parameter, public deriv_norm_drhoa
integer, parameter, public deriv_rhob
integer, parameter, public deriv_rhoa
integer, parameter, public deriv_tau
integer, parameter, public deriv_tau_b
integer, parameter, public deriv_tau_a
integer, parameter, public deriv_rho
integer, parameter, public deriv_norm_drhob
represent a group ofunctional derivatives
type(xc_derivative_type) function, pointer, public xc_dset_get_derivative(derivative_set, description, allocate_deriv)
returns the requested xc_derivative
Provides types for the management of the xc-functionals and their derivatives.
subroutine, public xc_derivative_get(deriv, split_desc, order, deriv_data, accept_null_data)
returns various information on the given derivative
subroutine, public xc_functionals_eval(functionals, lsd, rho_set, deriv_set, deriv_order)
...
contains the structure
contains the structure
subroutine, public xc_rho_set_get(rho_set, can_return_null, rho, drho, norm_drho, rhoa, rhob, norm_drhoa, norm_drhob, rho_1_3, rhoa_1_3, rhob_1_3, laplace_rho, laplace_rhoa, laplace_rhob, drhoa, drhob, rho_cutoff, drho_cutoff, tau_cutoff, tau, tau_a, tau_b, local_bounds)
returns the various attributes of rho_set
Exchange and Correlation functional calculations.
Definition xc.F:17
subroutine, public divide_by_norm_drho(deriv_set, rho_set, lsd)
divides derivatives from deriv_set by norm_drho
Definition xc.F:6239
subroutine, public xc_calc_2nd_deriv_analytical(v_xc, v_xc_tau, deriv_set, rho_set, rho1_set, pw_pool, xc_section, gapw, vxg, tddfpt_fac, compute_virial, virial_xc, spinflip)
Calculates the second derivative of E_xc at rho in the direction rho1 (if you see the second derivati...
Definition xc.F:2650
represent a single linked list that stores pointers to the elements
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
A derivative set contains the different derivatives of a xc-functional in form of a linked list.
represent a derivative of a functional
contains a flag for each component of xc_rho_set, so that you can use it to tell which components you...
represent a density, with all the representation and data needed to perform a functional evaluation