(git:0de0cc2)
xc_atom.F
Go to the documentation of this file.
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 
12  cp_sll_xc_deriv_type
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,&
20  USE xc_derivative_desc, ONLY: &
23  USE xc_derivative_set_types, ONLY: xc_derivative_set_type,&
26  xc_derivative_type
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 
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  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  REAL(dp), DIMENSION(:, :, :), POINTER :: deriv_data
97  REAL(kind=dp) :: drho_cutoff
98  TYPE(xc_derivative_type), POINTER :: deriv_att
99 
100  CALL timeset(routinen, handle)
101  my_only_energy = .false.
102  IF (PRESENT(energy_only)) my_only_energy = energy_only
103 
104  IF (PRESENT(adiabatic_rescale_factor)) THEN
105  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  my_epr_xc = .false.
112  IF (PRESENT(epr_xc)) my_epr_xc = epr_xc
113  my_deriv_order = deriv_order
114  IF (my_epr_xc) my_deriv_order = 2
115 
116  gradient_f = (needs%drho_spin .OR. needs%norm_drho_spin .OR. &
117  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  deriv_order=my_deriv_order)
125 
126  CALL xc_rho_set_get(rho_set, drho_cutoff=drho_cutoff)
127 
128  NULLIFY (deriv_data)
129 
130  IF (my_epr_xc) THEN
131  ! nabla v_xc (using the vxg arrays)
132  ! there's no point doing this when lsd = false
133  IF (lsd) THEN
134  deriv_att => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa])
135  IF (ASSOCIATED(deriv_att)) THEN
136  CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
137  DO ir = 1, nr
138  DO ia = 1, na
139  DO idir = 1, 3
140  vxg(idir, ia, ir, 1) = rho_set%drhoa(idir)%array(ia, ir, 1)* &
141  deriv_data(ia, ir, 1)
142  END DO !idir
143  END DO !ia
144  END DO !ir
145  NULLIFY (deriv_data)
146  END IF
147  deriv_att => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob])
148  IF (ASSOCIATED(deriv_att)) THEN
149  CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
150  DO ir = 1, nr
151  DO ia = 1, na
152  DO idir = 1, 3
153  vxg(idir, ia, ir, 2) = rho_set%drhob(idir)%array(ia, ir, 1)* &
154  deriv_data(ia, ir, 1)
155  END DO !idir
156  END DO !ia
157  END DO !ir
158  NULLIFY (deriv_data)
159  END IF
160  END IF
161  ! EXC energy ! is that needed for epr?
162  deriv_att => xc_dset_get_derivative(deriv_set, [INTEGER::])
163  exc = 0.0_dp
164  IF (ASSOCIATED(deriv_att)) THEN
165  CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
166  DO ir = 1, nr
167  DO ia = 1, na
168  exc = exc + deriv_data(ia, ir, 1)*w(ia, ir)
169  END DO
170  END DO
171  NULLIFY (deriv_data)
172  END IF
173  ELSE
174 ! EXC energy
175  deriv_att => xc_dset_get_derivative(deriv_set, [INTEGER::])
176  exc = 0.0_dp
177  IF (ASSOCIATED(deriv_att)) THEN
178  CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
179  DO ir = 1, nr
180  DO ia = 1, na
181  exc = exc + deriv_data(ia, ir, 1)*w(ia, ir)
182  END DO
183  END DO
184  NULLIFY (deriv_data)
185  END IF
186  ! Calculate the potential only if needed
187  IF (.NOT. my_only_energy) THEN
188 ! Derivative with respect to the density
189  IF (lsd) THEN
190  deriv_att => xc_dset_get_derivative(deriv_set, [deriv_rhoa])
191  IF (ASSOCIATED(deriv_att)) THEN
192  CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
193  vxc(:, :, 1) = deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
194  NULLIFY (deriv_data)
195  END IF
196  deriv_att => xc_dset_get_derivative(deriv_set, [deriv_rhob])
197  IF (ASSOCIATED(deriv_att)) THEN
198  CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
199  vxc(:, :, 2) = deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
200  NULLIFY (deriv_data)
201  END IF
202  deriv_att => xc_dset_get_derivative(deriv_set, [deriv_rho])
203  IF (ASSOCIATED(deriv_att)) THEN
204  CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
205  vxc(:, :, 1) = vxc(:, :, 1) + deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
206  vxc(:, :, 2) = vxc(:, :, 2) + deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
207  NULLIFY (deriv_data)
208  END IF
209  ELSE
210  deriv_att => xc_dset_get_derivative(deriv_set, [deriv_rho])
211  IF (ASSOCIATED(deriv_att)) THEN
212  CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
213  vxc(:, :, 1) = deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
214  NULLIFY (deriv_data)
215  END IF
216  END IF ! lsd
217 
218 ! Derivatives with respect to the gradient
219  IF (lsd) THEN
220  deriv_att => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa])
221  IF (ASSOCIATED(deriv_att)) THEN
222  CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
223  DO ir = 1, nr
224  DO ia = 1, na
225  DO idir = 1, 3
226  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  rho_set%norm_drhoa(ia, ir, 1)*my_adiabatic_rescale_factor
230  ELSE
231  vxg(idir, ia, ir, 1) = 0.0_dp
232  END IF
233  END DO
234  END DO
235  END DO
236  NULLIFY (deriv_data)
237  END IF
238  deriv_att => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob])
239  IF (ASSOCIATED(deriv_att)) THEN
240  CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
241  DO ir = 1, nr
242  DO ia = 1, na
243  DO idir = 1, 3
244  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  rho_set%norm_drhob(ia, ir, 1)*my_adiabatic_rescale_factor
248  ELSE
249  vxg(idir, ia, ir, 2) = 0.0_dp
250  END IF
251  END DO
252  END DO
253  END DO
254  NULLIFY (deriv_data)
255  END IF
256  ! Cross Terms
257  deriv_att => xc_dset_get_derivative(deriv_set, [deriv_norm_drho])
258  IF (ASSOCIATED(deriv_att)) THEN
259  CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
260  DO ir = 1, nr
261  DO ia = 1, na
262  DO idir = 1, 3
263  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  my_adiabatic_rescale_factor
270  END IF
271  END DO
272  END DO
273  END DO
274  NULLIFY (deriv_data)
275  END IF
276  ELSE
277  deriv_att => xc_dset_get_derivative(deriv_set, [deriv_norm_drho])
278  IF (ASSOCIATED(deriv_att)) THEN
279  CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
280  DO ir = 1, nr
281  DO ia = 1, na
282  IF (rho_set%norm_drho(ia, ir, 1) > drho_cutoff) THEN
283  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  rho_set%norm_drho(ia, ir, 1)*my_adiabatic_rescale_factor
287  END DO
288  ELSE
289  vxg(1:3, ia, ir, 1) = 0.0_dp
290  END IF
291  END DO
292  END DO
293  NULLIFY (deriv_data)
294  END IF
295  END IF ! lsd
296 ! Derivative with respect to tau
297  IF (lsd) THEN
298  deriv_att => xc_dset_get_derivative(deriv_set, [deriv_tau_a])
299  IF (ASSOCIATED(deriv_att)) THEN
300  CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
301  vtau(:, :, 1) = deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
302  NULLIFY (deriv_data)
303  END IF
304  deriv_att => xc_dset_get_derivative(deriv_set, [deriv_tau_b])
305  IF (ASSOCIATED(deriv_att)) THEN
306  CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
307  vtau(:, :, 2) = deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
308  NULLIFY (deriv_data)
309  END IF
310  deriv_att => xc_dset_get_derivative(deriv_set, [deriv_tau])
311  IF (ASSOCIATED(deriv_att)) THEN
312  CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
313  vtau(:, :, 1) = vtau(:, :, 1) + deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
314  vtau(:, :, 2) = vtau(:, :, 2) + deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
315  NULLIFY (deriv_data)
316  END IF
317  ELSE
318  deriv_att => xc_dset_get_derivative(deriv_set, [deriv_tau])
319  IF (ASSOCIATED(deriv_att)) THEN
320  CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
321  vtau(:, :, 1) = deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
322  NULLIFY (deriv_data)
323  END IF
324  END IF ! lsd
325  END IF ! only_energy
326  END IF ! epr_xc
327 
328  CALL timestop(handle)
329 
330  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  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  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  CALL timeset(routinen, handle)
373 
374  nspins = SIZE(vxc, 3)
375  lsd = (nspins == 2)
376  IF (ASSOCIATED(rho_set%rhoa)) THEN
377  lsd = .true.
378  END IF
379  my_fac_triplet = 1.0_dp
380  IF ((PRESENT(do_triplet)) .AND. (do_triplet)) my_fac_triplet = -1.0_dp
381 
382  CALL xc_rho_set_get(rho_set, drho_cutoff=drho_cutoff)
383  xc_fun_section => section_vals_get_subs_vals(xc_section, &
384  "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  deriv_order=2)
392 
393  CALL divide_by_norm_drho(deriv_set, rho_set, lsd)
394 
395  ! multiply by w
396  pos => deriv_set%derivs
397  DO WHILE (cp_sll_xc_deriv_next(pos, el_att=deriv_att))
398  deriv_att%deriv_data(:, :, 1) = w(:, :)*deriv_att%deriv_data(:, :, 1)
399  END DO
400 
401  NULLIFY (pw_pool)
402  ALLOCATE (vxc_pw(nspins))
403  DO ispin = 1, nspins
404  vxc_pw(ispin)%array => vxc(:, :, ispin:ispin)
405  END DO
406 
407  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  xc_section, gapw=.true., vxg=vxg, tddfpt_fac=my_fac_triplet)
411 
412  DEALLOCATE (vxc_pw)
413 
414  ! zero the derivative data for the next call
415  pos => deriv_set%derivs
416  DO WHILE (cp_sll_xc_deriv_next(pos, el_att=deriv_att))
417  deriv_att%deriv_data = 0.0_dp
418  END DO
419 
420  CALL timestop(handle)
421 
422  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  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  SELECT CASE (nspins)
448  CASE (1)
449 ! What is this for?
450  IF (needs%rho_1_3) THEN
451  NULLIFY (rho_set%rho_1_3)
452  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  rho_set%owns%rho_1_3 = .true.
454  rho_set%has%rho_1_3 = .false.
455  END IF
456 ! Allocate the storage space for the density
457  IF (needs%rho) THEN
458  NULLIFY (rho_set%rho)
459  ALLOCATE (rho_set%rho(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
460  rho_set%owns%rho = .true.
461  rho_set%has%rho = .false.
462  END IF
463 ! Allocate the storage space for the norm of the gradient of the density
464  IF (needs%norm_drho) THEN
465  NULLIFY (rho_set%norm_drho)
466  ALLOCATE (rho_set%norm_drho(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
467  rho_set%owns%norm_drho = .true.
468  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  IF (needs%drho) THEN
472  DO idir = 1, 3
473  NULLIFY (rho_set%drho(idir)%array)
474  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  rho_set%owns%drho = .true.
477  rho_set%has%drho = .false.
478  END IF
479  CASE (2)
480 ! Allocate the storage space for the total density
481  IF (needs%rho) THEN
482  ! this should never be the case unless you use LDA functionals with LSD
483  NULLIFY (rho_set%rho)
484  ALLOCATE (rho_set%rho(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
485  rho_set%owns%rho = .true.
486  rho_set%has%rho = .false.
487  END IF
488 ! What is this for?
489  IF (needs%rho_1_3) THEN
490  NULLIFY (rho_set%rho_1_3)
491  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  rho_set%owns%rho_1_3 = .true.
493  rho_set%has%rho_1_3 = .false.
494  END IF
495 ! What is this for?
496  IF (needs%rho_spin_1_3) THEN
497  NULLIFY (rho_set%rhoa_1_3, rho_set%rhob_1_3)
498  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  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  rho_set%owns%rho_spin_1_3 = .true.
501  rho_set%has%rho_spin_1_3 = .false.
502  END IF
503 ! Allocate the storage space for the spin densities rhoa and rhob
504  IF (needs%rho_spin) THEN
505  NULLIFY (rho_set%rhoa, rho_set%rhob)
506  ALLOCATE (rho_set%rhoa(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
507  ALLOCATE (rho_set%rhob(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
508  rho_set%owns%rho_spin = .true.
509  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  IF (needs%norm_drho) THEN
513  NULLIFY (rho_set%norm_drho)
514  ALLOCATE (rho_set%norm_drho(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
515  rho_set%owns%norm_drho = .true.
516  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  IF (needs%norm_drho_spin) THEN
520  NULLIFY (rho_set%norm_drhoa, rho_set%norm_drhob)
521  ALLOCATE (rho_set%norm_drhoa(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
522  ALLOCATE (rho_set%norm_drhob(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
523  rho_set%owns%norm_drho_spin = .true.
524  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  IF (needs%drho) THEN
528  DO idir = 1, 3
529  NULLIFY (rho_set%drho(idir)%array)
530  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  rho_set%owns%drho = .true.
533  rho_set%has%drho = .false.
534  END IF
535 ! Allocate the storage space for the components of the gradient for rhoa and rhob
536  IF (needs%drho_spin) THEN
537  DO idir = 1, 3
538  NULLIFY (rho_set%drhoa(idir)%array, rho_set%drhob(idir)%array)
539  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  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  rho_set%owns%drho_spin = .true.
543  rho_set%has%drho_spin = .false.
544  END IF
545 !
546  END SELECT
547 
548  ! tau part
549  IF (needs%tau) THEN
550  NULLIFY (rho_set%tau)
551  ALLOCATE (rho_set%tau(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
552  rho_set%owns%tau = .true.
553  END IF
554  IF (needs%tau_spin) THEN
555  NULLIFY (rho_set%tau_a, rho_set%tau_b)
556  ALLOCATE (rho_set%tau_a(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
557  ALLOCATE (rho_set%tau_b(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
558  rho_set%owns%tau_spin = .true.
559  rho_set%has%tau_spin = .false.
560  END IF
561 
562  ! Laplace part
563  IF (needs%laplace_rho) THEN
564  NULLIFY (rho_set%laplace_rho)
565  ALLOCATE (rho_set%laplace_rho(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
566  rho_set%owns%laplace_rho = .true.
567  END IF
568  IF (needs%laplace_rho_spin) THEN
569  NULLIFY (rho_set%laplace_rhoa)
570  NULLIFY (rho_set%laplace_rhob)
571  ALLOCATE (rho_set%laplace_rhoa(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
572  ALLOCATE (rho_set%laplace_rhob(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
573  rho_set%owns%laplace_rho_spin = .true.
574  rho_set%has%laplace_rho_spin = .true.
575  END IF
576 
577  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  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  my_nspins = nspins
608  tddft_split = .false.
609  IF (lsd .AND. nspins == 1) THEN
610  my_nspins = 2
611  tddft_split = .true.
612  END IF
613 
614  ! some checks
615  IF (lsd) THEN
616  ELSE
617  cpassert(SIZE(rho, 3) == 1)
618  END IF
619  SELECT CASE (my_nspins)
620  CASE (1)
621  cpassert(.NOT. needs%rho_spin)
622  cpassert(.NOT. needs%drho_spin)
623  cpassert(.NOT. needs%norm_drho_spin)
624  cpassert(.NOT. needs%rho_spin_1_3)
625  CASE (2)
626  CASE default
627  cpabort("Unsupported number of spins")
628  END SELECT
629 
630  gradient_f = (needs%drho_spin .OR. needs%norm_drho_spin .OR. &
631  needs%drho .OR. needs%norm_drho)
632 
633  SELECT CASE (my_nspins)
634  CASE (1)
635  ! Give rho to 1/3
636  IF (needs%rho_1_3) THEN
637  DO ia = 1, na
638  rho_set%rho_1_3(ia, ir, 1) = max(rho(ia, ir, 1), 0.0_dp)**f13
639  END DO
640  rho_set%owns%rho_1_3 = .true.
641  rho_set%has%rho_1_3 = .true.
642  END IF
643  ! Give the density
644  IF (needs%rho) THEN
645  DO ia = 1, na
646  rho_set%rho(ia, ir, 1) = rho(ia, ir, 1)
647  END DO
648  rho_set%owns%rho = .true.
649  rho_set%has%rho = .true.
650  END IF
651  ! Give the norm of the gradient of the density
652  IF (needs%norm_drho) THEN
653  DO ia = 1, na
654  rho_set%norm_drho(ia, ir, 1) = drho(4, ia, ir, 1)
655  END DO
656  rho_set%owns%norm_drho = .true.
657  rho_set%has%norm_drho = .true.
658  END IF
659  ! Give the three components of the gradient of the density
660  IF (needs%drho) THEN
661  DO idir = 1, 3
662  DO ia = 1, na
663  rho_set%drho(idir)%array(ia, ir, 1) = drho(idir, ia, ir, 1)
664  END DO
665  END DO
666  rho_set%owns%drho = .true.
667  rho_set%has%drho = .true.
668  END IF
669  CASE (2)
670  ! Give the total density
671  IF (needs%rho) THEN
672  ! this should never be the case unless you use LDA functionals with LSD
673  IF (.NOT. tddft_split) THEN
674  DO ia = 1, na
675  rho_set%rho(ia, ir, 1) = rho(ia, ir, 1) + rho(ia, ir, 2)
676  END DO
677  ELSE
678  DO ia = 1, na
679  rho_set%rho(ia, ir, 1) = rho(ia, ir, 1)
680  END DO
681  END IF
682  rho_set%owns%rho = .true.
683  rho_set%has%rho = .true.
684  END IF
685  ! Give the total density to 1/3
686  IF (needs%rho_1_3) THEN
687  IF (.NOT. tddft_split) THEN
688  DO ia = 1, na
689  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  DO ia = 1, na
693  rho_set%rho_1_3(ia, ir, 1) = max(rho(ia, ir, 1), 0.0_dp)**f13
694  END DO
695  END IF
696  rho_set%owns%rho_1_3 = .true.
697  rho_set%has%rho_1_3 = .true.
698  END IF
699  ! Give the spin densities to 1/3
700  IF (needs%rho_spin_1_3) THEN
701  IF (.NOT. tddft_split) THEN
702  DO ia = 1, na
703  rho_set%rhoa_1_3(ia, ir, 1) = max(rho(ia, ir, 1), 0.0_dp)**f13
704  rho_set%rhob_1_3(ia, ir, 1) = max(rho(ia, ir, 2), 0.0_dp)**f13
705  END DO
706  ELSE
707  DO ia = 1, na
708  rho_set%rhoa_1_3(ia, ir, 1) = max(0.5_dp*rho(ia, ir, 1), 0.0_dp)**f13
709  rho_set%rhob_1_3(ia, ir, 1) = rho_set%rhoa_1_3(ia, ir, 1)
710  END DO
711  END IF
712  rho_set%owns%rho_spin_1_3 = .true.
713  rho_set%has%rho_spin_1_3 = .true.
714  END IF
715  ! Give the spin densities rhoa and rhob
716  IF (needs%rho_spin) THEN
717  IF (.NOT. tddft_split) THEN
718  DO ia = 1, na
719  rho_set%rhoa(ia, ir, 1) = rho(ia, ir, 1)
720  rho_set%rhob(ia, ir, 1) = rho(ia, ir, 2)
721  END DO
722  ELSE
723  DO ia = 1, na
724  rho_set%rhoa(ia, ir, 1) = 0.5_dp*rho(ia, ir, 1)
725  rho_set%rhob(ia, ir, 1) = rho_set%rhoa(ia, ir, 1)
726  END DO
727  END IF
728  rho_set%owns%rho_spin = .true.
729  rho_set%has%rho_spin = .true.
730  END IF
731  ! Give the norm of the gradient of the total density
732  IF (needs%norm_drho) THEN
733  IF (.NOT. tddft_split) THEN
734  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  (drho(3, ia, ir, 1) + drho(3, ia, ir, 2))**2)
739  END DO
740  ELSE
741  DO ia = 1, na
742  rho_set%norm_drho(ia, ir, 1) = drho(4, ia, ir, 1)
743  END DO
744  END IF
745  rho_set%owns%norm_drho = .true.
746  rho_set%has%norm_drho = .true.
747  END IF
748  ! Give the norm of the gradient of rhoa and of rhob separatedly
749  IF (needs%norm_drho_spin) THEN
750  IF (.NOT. tddft_split) THEN
751  DO ia = 1, na
752  rho_set%norm_drhoa(ia, ir, 1) = drho(4, ia, ir, 1)
753  rho_set%norm_drhob(ia, ir, 1) = drho(4, ia, ir, 2)
754  END DO
755  ELSE
756  DO ia = 1, na
757  rho_set%norm_drhoa(ia, ir, 1) = 0.5_dp*drho(4, ia, ir, 1)
758  rho_set%norm_drhob(ia, ir, 1) = rho_set%norm_drhoa(ia, ir, 1)
759  END DO
760  END IF
761  rho_set%owns%norm_drho_spin = .true.
762  rho_set%has%norm_drho_spin = .true.
763  END IF
764  ! Give the components of the gradient for the total rho
765  IF (needs%drho) THEN
766  IF (.NOT. tddft_split) THEN
767  DO idir = 1, 3
768  DO ia = 1, na
769  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  DO idir = 1, 3
774  DO ia = 1, na
775  rho_set%drho(idir)%array(ia, ir, 1) = drho(idir, ia, ir, 1)
776  END DO
777  END DO
778  END IF
779  rho_set%owns%drho = .true.
780  rho_set%has%drho = .true.
781  END IF
782  ! Give the components of the gradient for rhoa and rhob
783  IF (needs%drho_spin) THEN
784  IF (.NOT. tddft_split) THEN
785  DO idir = 1, 3
786  DO ia = 1, na
787  rho_set%drhoa(idir)%array(ia, ir, 1) = drho(idir, ia, ir, 1)
788  rho_set%drhob(idir)%array(ia, ir, 1) = drho(idir, ia, ir, 2)
789  END DO
790  END DO
791  ELSE
792  DO idir = 1, 3
793  DO ia = 1, na
794  rho_set%drhoa(idir)%array(ia, ir, 1) = 0.5_dp*drho(idir, ia, ir, 1)
795  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  rho_set%owns%drho_spin = .true.
800  rho_set%has%drho_spin = .true.
801  END IF
802  !
803  END SELECT
804 
805  ! tau part
806  IF (needs%tau .OR. needs%tau_spin) THEN
807  cpassert(SIZE(tau, 3) == my_nspins)
808  END IF
809  IF (needs%tau) THEN
810  IF (my_nspins == 2) THEN
811  DO ia = 1, na
812  rho_set%tau(ia, ir, 1) = tau(ia, ir, 1) + tau(ia, ir, 2)
813  END DO
814  rho_set%owns%tau = .true.
815  rho_set%has%tau = .true.
816  ELSE
817  DO ia = 1, na
818  rho_set%tau(ia, ir, 1) = tau(ia, ir, 1)
819  END DO
820  rho_set%owns%tau = .true.
821  rho_set%has%tau = .true.
822  END IF
823  END IF
824  IF (needs%tau_spin) THEN
825  DO ia = 1, na
826  rho_set%tau_a(ia, ir, 1) = tau(ia, ir, 1)
827  rho_set%tau_b(ia, ir, 1) = tau(ia, ir, 2)
828  END DO
829  rho_set%owns%tau_spin = .true.
830  rho_set%has%tau_spin = .true.
831  END IF
832 
833  END SUBROUTINE fill_rho_set
834 
835 END 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 ...
Definition: pw_pool_types.F:24
Definition: xc_atom.F:9
subroutine, public xc_rho_set_atom_update(rho_set, needs, nspins, bo)
...
Definition: xc_atom.F:432
subroutine, public fill_rho_set(rho_set, lsd, nspins, needs, rho, drho, tau, na, ir)
...
Definition: xc_atom.F:592
subroutine, public xc_2nd_deriv_of_r(rho_set, rho1_set, xc_section, deriv_set, w, vxc, vxg, do_triplet)
...
Definition: xc_atom.F:345
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, epr_xc, adiabatic_rescale_factor)
...
Definition: xc_atom.F:65
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:5412
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)
Calculates the second derivative of E_xc at rho in the direction rho1 (if you see the second derivati...
Definition: xc.F:2635