(git:d5c4d39)
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-2026 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(:, :), INTENT(IN) :: w
81 LOGICAL, INTENT(IN) :: lsd
82 INTEGER, INTENT(in) :: na, nr
83 REAL(dp) :: exc
84 REAL(dp), DIMENSION(:, :, :), POINTER :: vxc
85 REAL(dp), DIMENSION(:, :, :, :), POINTER :: vxg
86 REAL(dp), DIMENSION(:, :, :), POINTER :: vtau
87 LOGICAL, INTENT(IN), OPTIONAL :: energy_only
88 REAL(dp), INTENT(IN), OPTIONAL :: adiabatic_rescale_factor
89
90 CHARACTER(LEN=*), PARAMETER :: routinen = 'vxc_of_r_new'
91
92 INTEGER :: handle, ia, idir, ir
93 LOGICAL :: gradient_f, my_only_energy
94 REAL(dp) :: my_adiabatic_rescale_factor
95 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(:, :), INTENT(IN) :: w
303 LOGICAL, INTENT(IN) :: lsd
304 INTEGER, INTENT(in) :: na, nr
305 REAL(dp) :: exc
306 REAL(dp), DIMENSION(:, :, :), POINTER :: vxc
307 REAL(dp), DIMENSION(:, :, :, :), POINTER :: vxg
308 REAL(dp), DIMENSION(:, :, :), POINTER :: vtau
309
310 CHARACTER(LEN=*), PARAMETER :: routinen = 'vxc_of_r_epr'
311
312 INTEGER :: handle, ia, idir, ir, my_deriv_order
313 LOGICAL :: gradient_f
314 REAL(dp) :: my_adiabatic_rescale_factor
315 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 vtau ...
398!> \param do_triplet ...
399!> \param do_sf ...
400! **************************************************************************************************
401 SUBROUTINE xc_2nd_deriv_of_r(rho_set, rho1_set, xc_section, &
402 deriv_set, w, vxc, vxg, vtau, do_triplet, do_sf)
403
404! As input of this routine one gets rho and drho on a one dimensional grid.
405! The grid is the angular grid corresponding to a given point ir on the radial grid.
406! The derivatives are calculated on this one dimensional grid, the results are stored in
407! vxc(1:na,ir,ispin), vxg(1:na,ir,ispin), vxg_cross(1:na,ir,ispin)
408! Afterwords the arrays containing the derivatives are put to zero so that the routine
409! can safely be called for the next radial point ir
410
411 TYPE(xc_rho_set_type), INTENT(IN) :: rho_set, rho1_set
412 TYPE(section_vals_type), POINTER :: xc_section
413 TYPE(xc_derivative_set_type), INTENT(INOUT) :: deriv_set
414 REAL(dp), DIMENSION(:, :), INTENT(IN) :: w
415 REAL(dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: vxc
416 REAL(dp), DIMENSION(:, :, :, :), POINTER :: vxg
417 REAL(dp), CONTIGUOUS, DIMENSION(:, :, :), &
418 OPTIONAL, POINTER :: vtau
419 LOGICAL, INTENT(IN), OPTIONAL :: do_triplet, do_sf
420
421 CHARACTER(LEN=*), PARAMETER :: routinen = 'xc_2nd_deriv_of_r'
422
423 INTEGER :: handle, ispin, nspins
424 LOGICAL :: lsd, my_do_sf
425 REAL(dp) :: drho_cutoff, my_fac_triplet
426 TYPE(cp_sll_xc_deriv_type), POINTER :: pos
427 TYPE(pw_pool_type), POINTER :: pw_pool
428 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: vxc_pw, vxc_tau_pw
429 TYPE(section_vals_type), POINTER :: xc_fun_section
430 TYPE(xc_derivative_type), POINTER :: deriv_att
431
432 CALL timeset(routinen, handle)
433
434 nspins = SIZE(vxc, 3)
435 lsd = (nspins == 2)
436 IF (ASSOCIATED(rho_set%rhoa)) THEN
437 lsd = .true.
438 END IF
439 my_fac_triplet = 1.0_dp
440 IF (PRESENT(do_triplet)) THEN
441 IF (do_triplet) my_fac_triplet = -1.0_dp
442 END IF
443
444 my_do_sf = .false.
445 IF (PRESENT(do_sf)) my_do_sf = do_sf
446
447 CALL xc_rho_set_get(rho_set, drho_cutoff=drho_cutoff)
448 xc_fun_section => section_vals_get_subs_vals(xc_section, &
449 "XC_FUNCTIONAL")
450
451 ! Calculate the derivatives
452 CALL xc_functionals_eval(xc_fun_section, &
453 lsd=lsd, &
454 rho_set=rho_set, &
455 deriv_set=deriv_set, &
456 deriv_order=2)
457
458 CALL divide_by_norm_drho(deriv_set, rho_set, lsd)
459
460 ! multiply by w
461 pos => deriv_set%derivs
462 DO WHILE (cp_sll_xc_deriv_next(pos, el_att=deriv_att))
463 deriv_att%deriv_data(:, :, 1) = w(:, :)*deriv_att%deriv_data(:, :, 1)
464 END DO
465
466 NULLIFY (pw_pool)
467 ALLOCATE (vxc_pw(nspins))
468 DO ispin = 1, nspins
469 vxc_pw(ispin)%array => vxc(:, :, ispin:ispin)
470 END DO
471
472 NULLIFY (vxc_tau_pw)
473 IF (PRESENT(vtau)) THEN
474 IF (ASSOCIATED(vtau)) THEN
475 ALLOCATE (vxc_tau_pw(nspins))
476 DO ispin = 1, nspins
477 vxc_tau_pw(ispin)%array => vtau(:, :, ispin:ispin)
478 END DO
479 END IF
480 END IF
481
482 CALL xc_calc_2nd_deriv_analytical(vxc_pw, vxc_tau_pw, deriv_set, rho_set, rho1_set, pw_pool, &
483 xc_section, gapw=.true., vxg=vxg, &
484 tddfpt_fac=my_fac_triplet, spinflip=do_sf)
485
486 DEALLOCATE (vxc_pw)
487 IF (ASSOCIATED(vxc_tau_pw)) DEALLOCATE (vxc_tau_pw)
488
489 ! zero the derivative data for the next call
490 pos => deriv_set%derivs
491 DO WHILE (cp_sll_xc_deriv_next(pos, el_att=deriv_att))
492 deriv_att%deriv_data = 0.0_dp
493 END DO
494
495 CALL timestop(handle)
496
497 END SUBROUTINE xc_2nd_deriv_of_r
498
499! **************************************************************************************************
500!> \brief ...
501!> \param rho_set ...
502!> \param needs ...
503!> \param nspins ...
504!> \param bo ...
505! **************************************************************************************************
506 SUBROUTINE xc_rho_set_atom_update(rho_set, needs, nspins, bo)
507
508! This routine allocates the storage arrays for rho and drho
509! In calculate_vxc_atom this is called once for each atomic_kind,
510! After the loop over all the atoms of the kind and over all the points
511! of the radial grid for each atom, rho_set is deallocated.
512! Within the same kind, at each new point on the radial grid, the rho_set
513! arrays rho and drho are overwritten.
514
515 TYPE(xc_rho_set_type), INTENT(INOUT) :: rho_set
516 TYPE(xc_rho_cflags_type), INTENT(IN) :: needs
517 INTEGER, INTENT(IN) :: nspins
518 INTEGER, DIMENSION(2, 3), INTENT(IN) :: bo
519
520 INTEGER :: idir
521
522 SELECT CASE (nspins)
523 CASE (1)
524! What is this for?
525 IF (needs%rho_1_3) THEN
526 NULLIFY (rho_set%rho_1_3)
527 ALLOCATE (rho_set%rho_1_3(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
528 rho_set%owns%rho_1_3 = .true.
529 rho_set%has%rho_1_3 = .false.
530 END IF
531! Allocate the storage space for the density
532 IF (needs%rho) THEN
533 NULLIFY (rho_set%rho)
534 ALLOCATE (rho_set%rho(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
535 rho_set%owns%rho = .true.
536 rho_set%has%rho = .false.
537 END IF
538! Allocate the storage space for the norm of the gradient of the density
539 IF (needs%norm_drho) THEN
540 NULLIFY (rho_set%norm_drho)
541 ALLOCATE (rho_set%norm_drho(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
542 rho_set%owns%norm_drho = .true.
543 rho_set%has%norm_drho = .false.
544 END IF
545! Allocate the storage space for the three components of the gradient of the density
546 IF (needs%drho) THEN
547 DO idir = 1, 3
548 NULLIFY (rho_set%drho(idir)%array)
549 ALLOCATE (rho_set%drho(idir)%array(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
550 END DO
551 rho_set%owns%drho = .true.
552 rho_set%has%drho = .false.
553 END IF
554 CASE (2)
555! Allocate the storage space for the total density
556 IF (needs%rho) THEN
557 ! this should never be the case unless you use LDA functionals with LSD
558 NULLIFY (rho_set%rho)
559 ALLOCATE (rho_set%rho(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
560 rho_set%owns%rho = .true.
561 rho_set%has%rho = .false.
562 END IF
563! What is this for?
564 IF (needs%rho_1_3) THEN
565 NULLIFY (rho_set%rho_1_3)
566 ALLOCATE (rho_set%rho_1_3(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
567 rho_set%owns%rho_1_3 = .true.
568 rho_set%has%rho_1_3 = .false.
569 END IF
570! What is this for?
571 IF (needs%rho_spin_1_3) THEN
572 NULLIFY (rho_set%rhoa_1_3, rho_set%rhob_1_3)
573 ALLOCATE (rho_set%rhoa_1_3(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
574 ALLOCATE (rho_set%rhob_1_3(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
575 rho_set%owns%rho_spin_1_3 = .true.
576 rho_set%has%rho_spin_1_3 = .false.
577 END IF
578! Allocate the storage space for the spin densities rhoa and rhob
579 IF (needs%rho_spin) THEN
580 NULLIFY (rho_set%rhoa, rho_set%rhob)
581 ALLOCATE (rho_set%rhoa(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
582 ALLOCATE (rho_set%rhob(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
583 rho_set%owns%rho_spin = .true.
584 rho_set%has%rho_spin = .false.
585 END IF
586! Allocate the storage space for the norm of the gradient of the total density
587 IF (needs%norm_drho) THEN
588 NULLIFY (rho_set%norm_drho)
589 ALLOCATE (rho_set%norm_drho(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
590 rho_set%owns%norm_drho = .true.
591 rho_set%has%norm_drho = .false.
592 END IF
593! Allocate the storage space for the norm of the gradient of rhoa and of rhob separatedly
594 IF (needs%norm_drho_spin) THEN
595 NULLIFY (rho_set%norm_drhoa, rho_set%norm_drhob)
596 ALLOCATE (rho_set%norm_drhoa(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
597 ALLOCATE (rho_set%norm_drhob(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
598 rho_set%owns%norm_drho_spin = .true.
599 rho_set%has%norm_drho_spin = .false.
600 END IF
601! Allocate the storage space for the components of the gradient for the total rho
602 IF (needs%drho) THEN
603 DO idir = 1, 3
604 NULLIFY (rho_set%drho(idir)%array)
605 ALLOCATE (rho_set%drho(idir)%array(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
606 END DO
607 rho_set%owns%drho = .true.
608 rho_set%has%drho = .false.
609 END IF
610! Allocate the storage space for the components of the gradient for rhoa and rhob
611 IF (needs%drho_spin) THEN
612 DO idir = 1, 3
613 NULLIFY (rho_set%drhoa(idir)%array, rho_set%drhob(idir)%array)
614 ALLOCATE (rho_set%drhoa(idir)%array(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
615 ALLOCATE (rho_set%drhob(idir)%array(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
616 END DO
617 rho_set%owns%drho_spin = .true.
618 rho_set%has%drho_spin = .false.
619 END IF
620!
621 END SELECT
622
623 ! tau part
624 IF (needs%tau) THEN
625 NULLIFY (rho_set%tau)
626 ALLOCATE (rho_set%tau(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
627 rho_set%owns%tau = .true.
628 END IF
629 IF (needs%tau_spin) THEN
630 NULLIFY (rho_set%tau_a, rho_set%tau_b)
631 ALLOCATE (rho_set%tau_a(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
632 ALLOCATE (rho_set%tau_b(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
633 rho_set%owns%tau_spin = .true.
634 rho_set%has%tau_spin = .false.
635 END IF
636
637 ! Laplace part
638 IF (needs%laplace_rho) THEN
639 NULLIFY (rho_set%laplace_rho)
640 ALLOCATE (rho_set%laplace_rho(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
641 rho_set%owns%laplace_rho = .true.
642 END IF
643 IF (needs%laplace_rho_spin) THEN
644 NULLIFY (rho_set%laplace_rhoa)
645 NULLIFY (rho_set%laplace_rhob)
646 ALLOCATE (rho_set%laplace_rhoa(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
647 ALLOCATE (rho_set%laplace_rhob(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
648 rho_set%owns%laplace_rho_spin = .true.
649 rho_set%has%laplace_rho_spin = .true.
650 END IF
651
652 END SUBROUTINE xc_rho_set_atom_update
653
654! **************************************************************************************************
655!> \brief ...
656!> \param rho_set ...
657!> \param lsd ...
658!> \param nspins ...
659!> \param needs ...
660!> \param rho ...
661!> \param drho ...
662!> \param tau ...
663!> \param na ...
664!> \param ir ...
665! **************************************************************************************************
666 SUBROUTINE fill_rho_set(rho_set, lsd, nspins, needs, rho, drho, tau, na, ir)
667
668 TYPE(xc_rho_set_type), INTENT(INOUT) :: rho_set
669 LOGICAL, INTENT(IN) :: lsd
670 INTEGER, INTENT(IN) :: nspins
671 TYPE(xc_rho_cflags_type), INTENT(IN) :: needs
672 REAL(dp), DIMENSION(:, :, :), INTENT(IN) :: rho
673 REAL(dp), DIMENSION(:, :, :, :), INTENT(IN) :: drho
674 REAL(dp), DIMENSION(:, :, :), INTENT(IN) :: tau
675 INTEGER, INTENT(IN) :: na, ir
676
677 REAL(kind=dp), PARAMETER :: f13 = (1.0_dp/3.0_dp)
678
679 INTEGER :: ia, idir, my_nspins
680 LOGICAL :: gradient_f, tddft_split
681
682 my_nspins = nspins
683 tddft_split = .false.
684 IF (lsd .AND. nspins == 1) THEN
685 my_nspins = 2
686 tddft_split = .true.
687 END IF
688
689 ! some checks
690 IF (lsd) THEN
691 ELSE
692 cpassert(SIZE(rho, 3) == 1)
693 END IF
694 SELECT CASE (my_nspins)
695 CASE (1)
696 cpassert(.NOT. needs%rho_spin)
697 cpassert(.NOT. needs%drho_spin)
698 cpassert(.NOT. needs%norm_drho_spin)
699 cpassert(.NOT. needs%rho_spin_1_3)
700 CASE (2)
701 CASE default
702 cpabort("Unsupported number of spins")
703 END SELECT
704
705 gradient_f = (needs%drho_spin .OR. needs%norm_drho_spin .OR. &
706 needs%drho .OR. needs%norm_drho)
707
708 SELECT CASE (my_nspins)
709 CASE (1)
710 ! Give rho to 1/3
711 IF (needs%rho_1_3) THEN
712 DO ia = 1, na
713 rho_set%rho_1_3(ia, ir, 1) = max(rho(ia, ir, 1), 0.0_dp)**f13
714 END DO
715 rho_set%owns%rho_1_3 = .true.
716 rho_set%has%rho_1_3 = .true.
717 END IF
718 ! Give the density
719 IF (needs%rho) THEN
720 DO ia = 1, na
721 rho_set%rho(ia, ir, 1) = rho(ia, ir, 1)
722 END DO
723 rho_set%owns%rho = .true.
724 rho_set%has%rho = .true.
725 END IF
726 ! Give the norm of the gradient of the density
727 IF (needs%norm_drho) THEN
728 DO ia = 1, na
729 rho_set%norm_drho(ia, ir, 1) = drho(4, ia, ir, 1)
730 END DO
731 rho_set%owns%norm_drho = .true.
732 rho_set%has%norm_drho = .true.
733 END IF
734 ! Give the three components of the gradient of the density
735 IF (needs%drho) THEN
736 DO idir = 1, 3
737 DO ia = 1, na
738 rho_set%drho(idir)%array(ia, ir, 1) = drho(idir, ia, ir, 1)
739 END DO
740 END DO
741 rho_set%owns%drho = .true.
742 rho_set%has%drho = .true.
743 END IF
744 CASE (2)
745 ! Give the total density
746 IF (needs%rho) THEN
747 ! this should never be the case unless you use LDA functionals with LSD
748 IF (.NOT. tddft_split) THEN
749 DO ia = 1, na
750 rho_set%rho(ia, ir, 1) = rho(ia, ir, 1) + rho(ia, ir, 2)
751 END DO
752 ELSE
753 DO ia = 1, na
754 rho_set%rho(ia, ir, 1) = rho(ia, ir, 1)
755 END DO
756 END IF
757 rho_set%owns%rho = .true.
758 rho_set%has%rho = .true.
759 END IF
760 ! Give the total density to 1/3
761 IF (needs%rho_1_3) THEN
762 IF (.NOT. tddft_split) THEN
763 DO ia = 1, na
764 rho_set%rho_1_3(ia, ir, 1) = max(rho(ia, ir, 1) + rho(ia, ir, 2), 0.0_dp)**f13
765 END DO
766 ELSE
767 DO ia = 1, na
768 rho_set%rho_1_3(ia, ir, 1) = max(rho(ia, ir, 1), 0.0_dp)**f13
769 END DO
770 END IF
771 rho_set%owns%rho_1_3 = .true.
772 rho_set%has%rho_1_3 = .true.
773 END IF
774 ! Give the spin densities to 1/3
775 IF (needs%rho_spin_1_3) THEN
776 IF (.NOT. tddft_split) THEN
777 DO ia = 1, na
778 rho_set%rhoa_1_3(ia, ir, 1) = max(rho(ia, ir, 1), 0.0_dp)**f13
779 rho_set%rhob_1_3(ia, ir, 1) = max(rho(ia, ir, 2), 0.0_dp)**f13
780 END DO
781 ELSE
782 DO ia = 1, na
783 rho_set%rhoa_1_3(ia, ir, 1) = max(0.5_dp*rho(ia, ir, 1), 0.0_dp)**f13
784 rho_set%rhob_1_3(ia, ir, 1) = rho_set%rhoa_1_3(ia, ir, 1)
785 END DO
786 END IF
787 rho_set%owns%rho_spin_1_3 = .true.
788 rho_set%has%rho_spin_1_3 = .true.
789 END IF
790 ! Give the spin densities rhoa and rhob
791 IF (needs%rho_spin) THEN
792 IF (.NOT. tddft_split) THEN
793 DO ia = 1, na
794 rho_set%rhoa(ia, ir, 1) = rho(ia, ir, 1)
795 rho_set%rhob(ia, ir, 1) = rho(ia, ir, 2)
796 END DO
797 ELSE
798 DO ia = 1, na
799 rho_set%rhoa(ia, ir, 1) = 0.5_dp*rho(ia, ir, 1)
800 rho_set%rhob(ia, ir, 1) = rho_set%rhoa(ia, ir, 1)
801 END DO
802 END IF
803 rho_set%owns%rho_spin = .true.
804 rho_set%has%rho_spin = .true.
805 END IF
806 ! Give the norm of the gradient of the total density
807 IF (needs%norm_drho) THEN
808 IF (.NOT. tddft_split) THEN
809 DO ia = 1, na
810 rho_set%norm_drho(ia, ir, 1) = sqrt( &
811 (drho(1, ia, ir, 1) + drho(1, ia, ir, 2))**2 + &
812 (drho(2, ia, ir, 1) + drho(2, ia, ir, 2))**2 + &
813 (drho(3, ia, ir, 1) + drho(3, ia, ir, 2))**2)
814 END DO
815 ELSE
816 DO ia = 1, na
817 rho_set%norm_drho(ia, ir, 1) = drho(4, ia, ir, 1)
818 END DO
819 END IF
820 rho_set%owns%norm_drho = .true.
821 rho_set%has%norm_drho = .true.
822 END IF
823 ! Give the norm of the gradient of rhoa and of rhob separatedly
824 IF (needs%norm_drho_spin) THEN
825 IF (.NOT. tddft_split) THEN
826 DO ia = 1, na
827 rho_set%norm_drhoa(ia, ir, 1) = drho(4, ia, ir, 1)
828 rho_set%norm_drhob(ia, ir, 1) = drho(4, ia, ir, 2)
829 END DO
830 ELSE
831 DO ia = 1, na
832 rho_set%norm_drhoa(ia, ir, 1) = 0.5_dp*drho(4, ia, ir, 1)
833 rho_set%norm_drhob(ia, ir, 1) = rho_set%norm_drhoa(ia, ir, 1)
834 END DO
835 END IF
836 rho_set%owns%norm_drho_spin = .true.
837 rho_set%has%norm_drho_spin = .true.
838 END IF
839 ! Give the components of the gradient for the total rho
840 IF (needs%drho) THEN
841 IF (.NOT. tddft_split) THEN
842 DO idir = 1, 3
843 DO ia = 1, na
844 rho_set%drho(idir)%array(ia, ir, 1) = drho(idir, ia, ir, 1) + drho(idir, ia, ir, 2)
845 END DO
846 END DO
847 ELSE
848 DO idir = 1, 3
849 DO ia = 1, na
850 rho_set%drho(idir)%array(ia, ir, 1) = drho(idir, ia, ir, 1)
851 END DO
852 END DO
853 END IF
854 rho_set%owns%drho = .true.
855 rho_set%has%drho = .true.
856 END IF
857 ! Give the components of the gradient for rhoa and rhob
858 IF (needs%drho_spin) THEN
859 IF (.NOT. tddft_split) THEN
860 DO idir = 1, 3
861 DO ia = 1, na
862 rho_set%drhoa(idir)%array(ia, ir, 1) = drho(idir, ia, ir, 1)
863 rho_set%drhob(idir)%array(ia, ir, 1) = drho(idir, ia, ir, 2)
864 END DO
865 END DO
866 ELSE
867 DO idir = 1, 3
868 DO ia = 1, na
869 rho_set%drhoa(idir)%array(ia, ir, 1) = 0.5_dp*drho(idir, ia, ir, 1)
870 rho_set%drhob(idir)%array(ia, ir, 1) = rho_set%drhoa(idir)%array(ia, ir, 1)
871 END DO
872 END DO
873 END IF
874 rho_set%owns%drho_spin = .true.
875 rho_set%has%drho_spin = .true.
876 END IF
877 !
878 END SELECT
879
880 ! tau part
881 IF (needs%tau .OR. needs%tau_spin) THEN
882 cpassert(SIZE(tau, 3) == my_nspins)
883 END IF
884 IF (needs%tau) THEN
885 IF (my_nspins == 2) THEN
886 DO ia = 1, na
887 rho_set%tau(ia, ir, 1) = tau(ia, ir, 1) + tau(ia, ir, 2)
888 END DO
889 rho_set%owns%tau = .true.
890 rho_set%has%tau = .true.
891 ELSE
892 DO ia = 1, na
893 rho_set%tau(ia, ir, 1) = tau(ia, ir, 1)
894 END DO
895 rho_set%owns%tau = .true.
896 rho_set%has%tau = .true.
897 END IF
898 END IF
899 IF (needs%tau_spin) THEN
900 DO ia = 1, na
901 rho_set%tau_a(ia, ir, 1) = tau(ia, ir, 1)
902 rho_set%tau_b(ia, ir, 1) = tau(ia, ir, 2)
903 END DO
904 rho_set%owns%tau_spin = .true.
905 rho_set%has%tau_spin = .true.
906 END IF
907
908 END SUBROUTINE fill_rho_set
909
910END 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:507
subroutine, public xc_2nd_deriv_of_r(rho_set, rho1_set, xc_section, deriv_set, w, vxc, vxg, vtau, do_triplet, do_sf)
...
Definition xc_atom.F:403
subroutine, public fill_rho_set(rho_set, lsd, nspins, needs, rho, drho, tau, na, ir)
...
Definition xc_atom.F:667
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:5659
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:2055
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