(git:34ef472)
qs_vxc_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 !> \brief routines that build the integrals of the Vxc potential calculated
10 !> for the atomic density in the basis set of spherical primitives
11 ! **************************************************************************************************
12 MODULE qs_vxc_atom
13  USE atomic_kind_types, ONLY: atomic_kind_type,&
14  get_atomic_kind
15  USE basis_set_types, ONLY: get_gto_basis_set,&
16  gto_basis_set_type
17  USE cp_control_types, ONLY: dft_control_type
18  USE input_constants, ONLY: tddfpt_excitations,&
19  tddfpt_triplet,&
20  xc_none
21  USE input_section_types, ONLY: section_vals_get_subs_vals,&
22  section_vals_type,&
23  section_vals_val_get
24  USE kinds, ONLY: dp
25  USE memory_utilities, ONLY: reallocate
26  USE message_passing, ONLY: mp_para_env_type
27  USE orbital_pointers, ONLY: indso,&
28  nsoset
29  USE paw_basis_types, ONLY: get_paw_basis_info
30  USE qs_environment_types, ONLY: get_qs_env,&
31  qs_environment_type
32  USE qs_grid_atom, ONLY: grid_atom_type
33  USE qs_harmonics_atom, ONLY: get_none0_cg_list,&
34  harmonics_atom_type
35  USE qs_kind_types, ONLY: get_qs_kind,&
36  has_nlcc,&
37  qs_kind_type
38  USE qs_linres_types, ONLY: nablavks_atom_type
39  USE qs_rho_atom_types, ONLY: get_rho_atom,&
40  rho_atom_coeff,&
41  rho_atom_type
42  USE util, ONLY: get_limit
43  USE xc_atom, ONLY: fill_rho_set,&
44  vxc_of_r_new,&
45  xc_2nd_deriv_of_r,&
46  xc_rho_set_atom_update
47  USE xc_derivative_set_types, ONLY: xc_derivative_set_type,&
48  xc_dset_create,&
49  xc_dset_release,&
50  xc_dset_zero_all
51  USE xc_derivatives, ONLY: xc_functionals_get_needs
52  USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
53  USE xc_rho_set_types, ONLY: xc_rho_set_create,&
54  xc_rho_set_release,&
55  xc_rho_set_type
56 #include "./base/base_uses.f90"
57 
58  IMPLICIT NONE
59 
60  PRIVATE
61 
62  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_vxc_atom'
63 
64  PUBLIC :: calculate_vxc_atom, &
65  calculate_xc_2nd_deriv_atom, &
66  calc_rho_angular, &
67  calculate_gfxc_atom, &
68  gfxc_atom_diff, &
69  gavxcgb_nogc
70 
71 CONTAINS
72 
73 ! **************************************************************************************************
74 !> \brief ...
75 !> \param qs_env ...
76 !> \param energy_only ...
77 !> \param exc1 the on-body ex energy contribution
78 !> \param gradient_atom_set ...
79 !> \param adiabatic_rescale_factor ...
80 !> \param kind_set_external provides a non-default kind_set to use
81 !> \param rho_atom_set_external provides a non-default atomic density set to use
82 !> \param xc_section_external provides an external non-default XC
83 ! **************************************************************************************************
84  SUBROUTINE calculate_vxc_atom(qs_env, energy_only, exc1, gradient_atom_set, &
85  adiabatic_rescale_factor, kind_set_external, &
86  rho_atom_set_external, xc_section_external)
87 
88  TYPE(qs_environment_type), POINTER :: qs_env
89  LOGICAL, INTENT(IN) :: energy_only
90  REAL(dp), INTENT(INOUT) :: exc1
91  TYPE(nablavks_atom_type), DIMENSION(:), OPTIONAL, &
92  POINTER :: gradient_atom_set
93  REAL(dp), INTENT(IN), OPTIONAL :: adiabatic_rescale_factor
94  TYPE(qs_kind_type), DIMENSION(:), OPTIONAL, &
95  POINTER :: kind_set_external
96  TYPE(rho_atom_type), DIMENSION(:), OPTIONAL, &
97  POINTER :: rho_atom_set_external
98  TYPE(section_vals_type), OPTIONAL, POINTER :: xc_section_external
99 
100  CHARACTER(LEN=*), PARAMETER :: routinen = 'calculate_vxc_atom'
101 
102  INTEGER :: bo(2), handle, ia, iat, iatom, idir, &
103  ikind, ir, ispin, myfun, na, natom, &
104  nr, nspins, num_pe
105  INTEGER, DIMENSION(2, 3) :: bounds
106  INTEGER, DIMENSION(:), POINTER :: atom_list
107  LOGICAL :: donlcc, epr_xc, gradient_f, lsd, nlcc, &
108  paw_atom, tau_f
109  REAL(dp) :: density_cut, exc_h, exc_s, gradient_cut, &
110  my_adiabatic_rescale_factor, tau_cut
111  REAL(dp), DIMENSION(1, 1, 1) :: tau_d
112  REAL(dp), DIMENSION(1, 1, 1, 1) :: rho_d
113  REAL(dp), DIMENSION(:, :), POINTER :: rho_nlcc, weight
114  REAL(dp), DIMENSION(:, :, :), POINTER :: rho_h, rho_s, tau_h, tau_s, vtau_h, &
115  vtau_s, vxc_h, vxc_s
116  REAL(dp), DIMENSION(:, :, :, :), POINTER :: drho_h, drho_s, vxg_h, vxg_s
117  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
118  TYPE(dft_control_type), POINTER :: dft_control
119  TYPE(grid_atom_type), POINTER :: grid_atom
120  TYPE(gto_basis_set_type), POINTER :: basis_1c
121  TYPE(harmonics_atom_type), POINTER :: harmonics
122  TYPE(mp_para_env_type), POINTER :: para_env
123  TYPE(qs_kind_type), DIMENSION(:), POINTER :: my_kind_set
124  TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: dr_h, dr_s, int_hh, int_ss, r_h, r_s
125  TYPE(rho_atom_coeff), DIMENSION(:, :), POINTER :: r_h_d, r_s_d
126  TYPE(rho_atom_type), DIMENSION(:), POINTER :: my_rho_atom_set
127  TYPE(rho_atom_type), POINTER :: rho_atom
128  TYPE(section_vals_type), POINTER :: input, my_xc_section, xc_fun_section
129  TYPE(xc_derivative_set_type) :: deriv_set
130  TYPE(xc_rho_cflags_type) :: needs
131  TYPE(xc_rho_set_type) :: rho_set_h, rho_set_s
132 
133 ! -------------------------------------------------------------------------
134 
135  CALL timeset(routinen, handle)
136 
137  NULLIFY (atom_list)
138  NULLIFY (my_kind_set)
139  NULLIFY (atomic_kind_set)
140  NULLIFY (grid_atom)
141  NULLIFY (harmonics)
142  NULLIFY (input)
143  NULLIFY (para_env)
144  NULLIFY (rho_atom)
145  NULLIFY (my_rho_atom_set)
146  NULLIFY (rho_nlcc)
147 
148  epr_xc = .false.
149  IF (PRESENT(gradient_atom_set)) THEN
150  epr_xc = .true.
151  END IF
152 
153  IF (PRESENT(adiabatic_rescale_factor)) THEN
154  my_adiabatic_rescale_factor = adiabatic_rescale_factor
155  ELSE
156  my_adiabatic_rescale_factor = 1.0_dp
157  END IF
158 
159  CALL get_qs_env(qs_env=qs_env, &
160  dft_control=dft_control, &
161  para_env=para_env, &
162  atomic_kind_set=atomic_kind_set, &
163  qs_kind_set=my_kind_set, &
164  input=input, &
165  rho_atom_set=my_rho_atom_set)
166 
167  IF (PRESENT(kind_set_external)) my_kind_set => kind_set_external
168  IF (PRESENT(rho_atom_set_external)) my_rho_atom_set => rho_atom_set_external
169 
170  nlcc = has_nlcc(my_kind_set)
171 
172  IF (epr_xc) THEN
173  my_xc_section => section_vals_get_subs_vals(input, &
174  "PROPERTIES%LINRES%EPR%PRINT%G_TENSOR%XC")
175  ELSE
176  my_xc_section => section_vals_get_subs_vals(input, "DFT%XC")
177  END IF
178 
179  IF (PRESENT(xc_section_external)) my_xc_section => xc_section_external
180 
181  xc_fun_section => section_vals_get_subs_vals(my_xc_section, "XC_FUNCTIONAL")
182  CALL section_vals_val_get(xc_fun_section, "_SECTION_PARAMETERS_", &
183  i_val=myfun)
184 
185  IF (myfun == xc_none) THEN
186  exc1 = 0.0_dp
187  my_rho_atom_set(:)%exc_h = 0.0_dp
188  my_rho_atom_set(:)%exc_s = 0.0_dp
189  ELSE
190  CALL section_vals_val_get(my_xc_section, "DENSITY_CUTOFF", &
191  r_val=density_cut)
192  CALL section_vals_val_get(my_xc_section, "GRADIENT_CUTOFF", &
193  r_val=gradient_cut)
194  CALL section_vals_val_get(my_xc_section, "TAU_CUTOFF", &
195  r_val=tau_cut)
196 
197  lsd = dft_control%lsd
198  nspins = dft_control%nspins
199  needs = xc_functionals_get_needs(xc_fun_section, &
200  lsd=lsd, &
201  calc_potential=.true.)
202 
203  ! whatever the xc, if epr_xc, drho_spin is needed
204  IF (epr_xc) needs%drho_spin = .true.
205 
206  gradient_f = (needs%drho .OR. needs%drho_spin)
207  tau_f = (needs%tau .OR. needs%tau_spin)
208 
209  ! Initialize energy contribution from the one center XC terms to zero
210  exc1 = 0.0_dp
211 
212  ! Nullify some pointers for work-arrays
213  NULLIFY (rho_h, drho_h, rho_s, drho_s, weight)
214  NULLIFY (vxc_h, vxc_s, vxg_h, vxg_s)
215  NULLIFY (tau_h, tau_s)
216  NULLIFY (vtau_h, vtau_s)
217 
218  ! Here starts the loop over all the atoms
219 
220  DO ikind = 1, SIZE(atomic_kind_set)
221  CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
222  CALL get_qs_kind(my_kind_set(ikind), paw_atom=paw_atom, &
223  harmonics=harmonics, grid_atom=grid_atom)
224  CALL get_qs_kind(my_kind_set(ikind), basis_set=basis_1c, basis_type="GAPW_1C")
225 
226  IF (.NOT. paw_atom) cycle
227 
228  nr = grid_atom%nr
229  na = grid_atom%ng_sphere
230 
231  ! Prepare the structures needed to calculate and store the xc derivatives
232 
233  ! Array dimension: here anly one dimensional arrays are used,
234  ! i.e. only the first column of deriv_data is read.
235  ! The other to dimensions are set to size equal 1
236  bounds(1:2, 1:3) = 1
237  bounds(2, 1) = na
238  bounds(2, 2) = nr
239 
240  ! create a place where to put the derivatives
241  CALL xc_dset_create(deriv_set, local_bounds=bounds)
242  ! create the place where to store the argument for the functionals
243  CALL xc_rho_set_create(rho_set_h, bounds, rho_cutoff=density_cut, &
244  drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
245  CALL xc_rho_set_create(rho_set_s, bounds, rho_cutoff=density_cut, &
246  drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
247 
248  ! allocate the required 3d arrays where to store rho and drho
249  CALL xc_rho_set_atom_update(rho_set_h, needs, nspins, bounds)
250  CALL xc_rho_set_atom_update(rho_set_s, needs, nspins, bounds)
251 
252  CALL reallocate(rho_h, 1, na, 1, nr, 1, nspins)
253  CALL reallocate(rho_s, 1, na, 1, nr, 1, nspins)
254  weight => grid_atom%weight
255  CALL reallocate(vxc_h, 1, na, 1, nr, 1, nspins)
256  CALL reallocate(vxc_s, 1, na, 1, nr, 1, nspins)
257  !
258  IF (gradient_f) THEN
259  CALL reallocate(drho_h, 1, 4, 1, na, 1, nr, 1, nspins)
260  CALL reallocate(drho_s, 1, 4, 1, na, 1, nr, 1, nspins)
261  CALL reallocate(vxg_h, 1, 3, 1, na, 1, nr, 1, nspins)
262  CALL reallocate(vxg_s, 1, 3, 1, na, 1, nr, 1, nspins)
263  END IF
264 
265  IF (tau_f) THEN
266  CALL reallocate(tau_h, 1, na, 1, nr, 1, nspins)
267  CALL reallocate(tau_s, 1, na, 1, nr, 1, nspins)
268  CALL reallocate(vtau_h, 1, na, 1, nr, 1, nspins)
269  CALL reallocate(vtau_s, 1, na, 1, nr, 1, nspins)
270  END IF
271 
272  ! NLCC: prepare rho and drho of the core charge for this KIND
273  donlcc = .false.
274  IF (nlcc) THEN
275  NULLIFY (rho_nlcc)
276  rho_nlcc => my_kind_set(ikind)%nlcc_pot
277  IF (ASSOCIATED(rho_nlcc)) donlcc = .true.
278  END IF
279 
280  ! Distribute the atoms of this kind
281 
282  num_pe = para_env%num_pe
283  bo = get_limit(natom, para_env%num_pe, para_env%mepos)
284 
285  DO iat = bo(1), bo(2)
286  iatom = atom_list(iat)
287 
288  my_rho_atom_set(iatom)%exc_h = 0.0_dp
289  my_rho_atom_set(iatom)%exc_s = 0.0_dp
290 
291  rho_atom => my_rho_atom_set(iatom)
292  rho_h = 0.0_dp
293  rho_s = 0.0_dp
294  IF (gradient_f) THEN
295  NULLIFY (r_h, r_s, dr_h, dr_s, r_h_d, r_s_d)
296  CALL get_rho_atom(rho_atom=rho_atom, rho_rad_h=r_h, &
297  rho_rad_s=r_s, drho_rad_h=dr_h, &
298  drho_rad_s=dr_s, rho_rad_h_d=r_h_d, &
299  rho_rad_s_d=r_s_d)
300  drho_h = 0.0_dp
301  drho_s = 0.0_dp
302  ELSE
303  NULLIFY (r_h, r_s)
304  CALL get_rho_atom(rho_atom=rho_atom, rho_rad_h=r_h, rho_rad_s=r_s)
305  rho_d = 0.0_dp
306  END IF
307  IF (tau_f) THEN
308  !compute tau on the grid all at once
309  CALL calc_tau_atom(tau_h, tau_s, rho_atom, my_kind_set(ikind), nspins)
310  ELSE
311  tau_d = 0.0_dp
312  END IF
313 
314  DO ir = 1, nr
315  CALL calc_rho_angular(grid_atom, harmonics, nspins, gradient_f, &
316  ir, r_h, r_s, rho_h, rho_s, dr_h, dr_s, &
317  r_h_d, r_s_d, drho_h, drho_s)
318  IF (donlcc) THEN
319  CALL calc_rho_nlcc(grid_atom, nspins, gradient_f, &
320  ir, rho_nlcc(:, 1), rho_h, rho_s, rho_nlcc(:, 2), drho_h, drho_s)
321  END IF
322  END DO
323  DO ir = 1, nr
324  IF (tau_f) THEN
325  CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, drho_h, tau_h, na, ir)
326  CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, drho_s, tau_s, na, ir)
327  ELSE IF (gradient_f) THEN
328  CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, drho_h, tau_d, na, ir)
329  CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, drho_s, tau_d, na, ir)
330  ELSE
331  CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, rho_d, tau_d, na, ir)
332  CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, rho_d, tau_d, na, ir)
333  END IF
334  END DO
335 
336  !-------------------!
337  ! hard atom density !
338  !-------------------!
339  CALL xc_dset_zero_all(deriv_set)
340  CALL vxc_of_r_new(xc_fun_section, rho_set_h, deriv_set, 1, needs, weight, &
341  lsd, na, nr, exc_h, vxc_h, vxg_h, vtau_h, energy_only=energy_only, &
342  epr_xc=epr_xc, adiabatic_rescale_factor=my_adiabatic_rescale_factor)
343  rho_atom%exc_h = rho_atom%exc_h + exc_h
344 
345  !-------------------!
346  ! soft atom density !
347  !-------------------!
348  CALL xc_dset_zero_all(deriv_set)
349  CALL vxc_of_r_new(xc_fun_section, rho_set_s, deriv_set, 1, needs, weight, &
350  lsd, na, nr, exc_s, vxc_s, vxg_s, vtau_s, energy_only=energy_only, &
351  epr_xc=epr_xc, adiabatic_rescale_factor=my_adiabatic_rescale_factor)
352  rho_atom%exc_s = rho_atom%exc_s + exc_s
353 
354  IF (epr_xc) THEN
355  DO ispin = 1, nspins
356  DO idir = 1, 3
357  DO ir = 1, nr
358  DO ia = 1, na
359  gradient_atom_set(iatom)%nablavks_vec_rad_h(idir, ispin)%r_coef(ir, ia) = &
360  gradient_atom_set(iatom)%nablavks_vec_rad_h(idir, ispin)%r_coef(ir, ia) &
361  + vxg_h(idir, ia, ir, ispin)
362  gradient_atom_set(iatom)%nablavks_vec_rad_s(idir, ispin)%r_coef(ir, ia) = &
363  gradient_atom_set(iatom)%nablavks_vec_rad_s(idir, ispin)%r_coef(ir, ia) &
364  + vxg_s(idir, ia, ir, ispin)
365  END DO ! ia
366  END DO ! ir
367  END DO ! idir
368  END DO ! ispin
369  END IF
370 
371  ! Add contributions to the exc energy
372 
373  exc1 = exc1 + rho_atom%exc_h - rho_atom%exc_s
374 
375  ! Integration to get the matrix elements relative to the vxc_atom
376  ! here the products with the primitives is done: gaVxcgb
377  ! internal transformation to get the integral in cartesian Gaussians
378 
379  IF (.NOT. energy_only) THEN
380  NULLIFY (int_hh, int_ss)
381  CALL get_rho_atom(rho_atom=rho_atom, ga_vlocal_gb_h=int_hh, ga_vlocal_gb_s=int_ss)
382  IF (gradient_f) THEN
383  CALL gavxcgb_gc(vxc_h, vxc_s, vxg_h, vxg_s, int_hh, int_ss, &
384  grid_atom, basis_1c, harmonics, nspins)
385  ELSE
386  CALL gavxcgb_nogc(vxc_h, vxc_s, int_hh, int_ss, &
387  grid_atom, basis_1c, harmonics, nspins)
388  END IF
389  IF (tau_f) THEN
390  CALL dgavtaudgb(vtau_h, vtau_s, int_hh, int_ss, &
391  grid_atom, basis_1c, harmonics, nspins)
392  END IF
393  END IF ! energy_only
394  NULLIFY (r_h, r_s, dr_h, dr_s)
395  END DO ! iat
396 
397  ! Release the xc structure used to store the xc derivatives
398  CALL xc_dset_release(deriv_set)
399  CALL xc_rho_set_release(rho_set_h)
400  CALL xc_rho_set_release(rho_set_s)
401 
402  END DO ! ikind
403 
404  CALL para_env%sum(exc1)
405 
406  IF (ASSOCIATED(rho_h)) DEALLOCATE (rho_h)
407  IF (ASSOCIATED(rho_s)) DEALLOCATE (rho_s)
408  IF (ASSOCIATED(vxc_h)) DEALLOCATE (vxc_h)
409  IF (ASSOCIATED(vxc_s)) DEALLOCATE (vxc_s)
410 
411  IF (gradient_f) THEN
412  IF (ASSOCIATED(drho_h)) DEALLOCATE (drho_h)
413  IF (ASSOCIATED(drho_s)) DEALLOCATE (drho_s)
414  IF (ASSOCIATED(vxg_h)) DEALLOCATE (vxg_h)
415  IF (ASSOCIATED(vxg_s)) DEALLOCATE (vxg_s)
416  END IF
417 
418  IF (tau_f) THEN
419  IF (ASSOCIATED(tau_h)) DEALLOCATE (tau_h)
420  IF (ASSOCIATED(tau_s)) DEALLOCATE (tau_s)
421  IF (ASSOCIATED(vtau_h)) DEALLOCATE (vtau_h)
422  IF (ASSOCIATED(vtau_s)) DEALLOCATE (vtau_s)
423  END IF
424 
425  END IF !xc_none
426 
427  CALL timestop(handle)
428 
429  END SUBROUTINE calculate_vxc_atom
430 
431 ! **************************************************************************************************
432 !> \brief ...
433 !> \param rho_atom_set ...
434 !> \param rho1_atom_set ...
435 !> \param qs_env ...
436 !> \param xc_section ...
437 !> \param para_env ...
438 !> \param do_tddft Initial implementation of TDDFT. Control parameters are read directly from
439 !> 'DFT' input section
440 !> \param do_tddfpt2 New implementation of TDDFT.
441 !> \param do_triplet ...
442 !> \param kind_set_external ...
443 ! **************************************************************************************************
444  SUBROUTINE calculate_xc_2nd_deriv_atom(rho_atom_set, rho1_atom_set, qs_env, xc_section, para_env, &
445  do_tddft, do_tddfpt2, do_triplet, kind_set_external)
446 
447  TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set, rho1_atom_set
448  TYPE(qs_environment_type), POINTER :: qs_env
449  TYPE(section_vals_type), POINTER :: xc_section
450  TYPE(mp_para_env_type), INTENT(IN) :: para_env
451  LOGICAL, INTENT(IN), OPTIONAL :: do_tddft, do_tddfpt2, do_triplet
452  TYPE(qs_kind_type), DIMENSION(:), OPTIONAL, &
453  POINTER :: kind_set_external
454 
455  CHARACTER(LEN=*), PARAMETER :: routinen = 'calculate_xc_2nd_deriv_atom'
456 
457  INTEGER :: atom, excitations, handle, iatom, ikind, &
458  ir, na, natom, nr, nspins, res_etype
459  INTEGER, DIMENSION(2) :: local_loop_limit
460  INTEGER, DIMENSION(2, 3) :: bounds
461  INTEGER, DIMENSION(:), POINTER :: atom_list
462  LOGICAL :: gradient_functional, lsd, lsd_singlets, &
463  my_tddft, paw_atom, scale_rho, tau_f
464  REAL(kind=dp) :: density_cut, gradient_cut, rtot, tau_cut
465  REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
466  POINTER :: vxc_h, vxc_s
467  REAL(kind=dp), DIMENSION(1, 1, 1) :: rtau
468  REAL(kind=dp), DIMENSION(1, 1, 1, 1) :: rrho
469  REAL(kind=dp), DIMENSION(:, :), POINTER :: weight
470  REAL(kind=dp), DIMENSION(:, :, :), POINTER :: rho1_h, rho1_s, rho_h, rho_s, tau1_h, &
471  tau1_s, tau_h, tau_s
472  REAL(kind=dp), DIMENSION(:, :, :, :), POINTER :: drho1_h, drho1_s, drho_h, drho_s, vxg_h, &
473  vxg_s
474  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
475  TYPE(grid_atom_type), POINTER :: grid_atom
476  TYPE(gto_basis_set_type), POINTER :: basis_1c
477  TYPE(harmonics_atom_type), POINTER :: harmonics
478  TYPE(qs_kind_type), DIMENSION(:), POINTER :: my_kind_set, qs_kind_set
479  TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: dr1_h, dr1_s, dr_h, dr_s, int_hh, &
480  int_ss, r1_h, r1_s, r_h, r_s
481  TYPE(rho_atom_coeff), DIMENSION(:, :), POINTER :: r1_h_d, r1_s_d, r_h_d, r_s_d
482  TYPE(rho_atom_type), POINTER :: rho1_atom, rho_atom
483  TYPE(section_vals_type), POINTER :: input, xc_fun_section
484  TYPE(xc_derivative_set_type) :: deriv_set
485  TYPE(xc_rho_cflags_type) :: needs
486  TYPE(xc_rho_set_type) :: rho1_set_h, rho1_set_s, rho_set_h, &
487  rho_set_s
488 
489 ! -------------------------------------------------------------------------
490 
491  CALL timeset(routinen, handle)
492 
493  NULLIFY (qs_kind_set)
494  NULLIFY (rho_h, rho_s, drho_h, drho_s, weight)
495  NULLIFY (rho1_h, rho1_s, drho1_h, drho1_s)
496  NULLIFY (vxc_h, vxc_s, vxg_h, vxg_s)
497  NULLIFY (tau_h, tau_s, tau1_h, tau1_s)
498 
499  CALL get_qs_env(qs_env=qs_env, &
500  input=input, &
501  qs_kind_set=qs_kind_set, &
502  atomic_kind_set=atomic_kind_set)
503 
504  IF (PRESENT(kind_set_external)) THEN
505  my_kind_set => kind_set_external
506  ELSE
507  my_kind_set => qs_kind_set
508  END IF
509 
510  my_tddft = .false.
511  IF (PRESENT(do_tddft)) my_tddft = do_tddft
512  CALL section_vals_val_get(input, "DFT%LSD", l_val=lsd)
513  CALL section_vals_val_get(xc_section, "DENSITY_CUTOFF", &
514  r_val=density_cut)
515  CALL section_vals_val_get(xc_section, "GRADIENT_CUTOFF", &
516  r_val=gradient_cut)
517  CALL section_vals_val_get(xc_section, "TAU_CUTOFF", &
518  r_val=tau_cut)
519  IF (my_tddft) THEN
520  CALL section_vals_val_get(input, "DFT%EXCITATIONS", &
521  i_val=excitations)
522  CALL section_vals_val_get(input, "DFT%TDDFPT%LSD_SINGLETS", &
523  l_val=lsd_singlets)
524  CALL section_vals_val_get(input, "DFT%TDDFPT%RES_ETYPE", &
525  i_val=res_etype)
526  END IF
527 
528  xc_fun_section => section_vals_get_subs_vals(xc_section, &
529  "XC_FUNCTIONAL")
530  IF (lsd) THEN
531  nspins = 2
532  ELSE
533  nspins = 1
534  END IF
535 
536  scale_rho = .false.
537  IF (my_tddft) THEN
538  IF (excitations == tddfpt_excitations) THEN
539  IF (nspins == 1 .AND. (lsd_singlets .OR. res_etype == tddfpt_triplet)) THEN
540  lsd = .true.
541  END IF
542  END IF
543  ELSEIF (PRESENT(do_tddfpt2) .AND. PRESENT(do_triplet)) THEN
544  IF (nspins == 1 .AND. do_triplet) THEN
545  lsd = .true.
546  scale_rho = .true.
547  END IF
548  ELSEIF (PRESENT(do_triplet)) THEN
549  IF (nspins == 1 .AND. do_triplet) lsd = .true.
550  END IF
551 
552  needs = xc_functionals_get_needs(xc_fun_section, lsd=lsd, &
553  calc_potential=.true.)
554  gradient_functional = needs%drho .OR. needs%drho_spin
555  tau_f = (needs%tau .OR. needs%tau_spin)
556  IF (tau_f) THEN
557  cpabort("Tau functionals not implemented for GAPW 2nd derivatives")
558  ELSE
559  rtau = 0.0_dp
560  END IF
561 
562  ! Here starts the loop over all the atoms
563  DO ikind = 1, SIZE(atomic_kind_set)
564 
565  NULLIFY (atom_list, harmonics, grid_atom)
566  CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
567  CALL get_qs_kind(my_kind_set(ikind), paw_atom=paw_atom, &
568  harmonics=harmonics, grid_atom=grid_atom)
569  CALL get_qs_kind(my_kind_set(ikind), basis_set=basis_1c, basis_type="GAPW_1C")
570  IF (.NOT. paw_atom) cycle
571 
572  nr = grid_atom%nr
573  na = grid_atom%ng_sphere
574 
575  ! Array dimension: here anly one dimensional arrays are used,
576  ! i.e. only the first column of deriv_data is read.
577  ! The other to dimensions are set to size equal 1.
578  bounds(1:2, 1:3) = 1
579  bounds(2, 1) = na
580  bounds(2, 2) = nr
581 
582  CALL xc_dset_create(deriv_set, local_bounds=bounds)
583  CALL xc_rho_set_create(rho_set_h, bounds, rho_cutoff=density_cut, &
584  drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
585  CALL xc_rho_set_create(rho_set_s, bounds, rho_cutoff=density_cut, &
586  drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
587  CALL xc_rho_set_create(rho1_set_h, bounds, rho_cutoff=density_cut, &
588  drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
589  CALL xc_rho_set_create(rho1_set_s, bounds, rho_cutoff=density_cut, &
590  drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
591 
592  ! allocate the required 3d arrays where to store rho and drho
593  IF (nspins == 1 .AND. .NOT. lsd) THEN
594  CALL xc_rho_set_atom_update(rho_set_h, needs, 1, bounds)
595  CALL xc_rho_set_atom_update(rho1_set_h, needs, 1, bounds)
596  CALL xc_rho_set_atom_update(rho_set_s, needs, 1, bounds)
597  CALL xc_rho_set_atom_update(rho1_set_s, needs, 1, bounds)
598  ELSE
599  CALL xc_rho_set_atom_update(rho_set_h, needs, 2, bounds)
600  CALL xc_rho_set_atom_update(rho1_set_h, needs, 2, bounds)
601  CALL xc_rho_set_atom_update(rho_set_s, needs, 2, bounds)
602  CALL xc_rho_set_atom_update(rho1_set_s, needs, 2, bounds)
603  END IF
604 
605  ALLOCATE (rho_h(1:na, 1:nr, 1:nspins), rho1_h(1:na, 1:nr, 1:nspins), &
606  rho_s(1:na, 1:nr, 1:nspins), rho1_s(1:na, 1:nr, 1:nspins))
607 
608  ALLOCATE (vxc_h(1:na, 1:nr, 1:nspins), vxc_s(1:na, 1:nr, 1:nspins))
609  vxc_h = 0.0_dp
610  vxc_s = 0.0_dp
611 
612  weight => grid_atom%weight
613 
614  IF (gradient_functional) THEN
615  ALLOCATE (drho_h(1:4, 1:na, 1:nr, 1:nspins), drho1_h(1:4, 1:na, 1:nr, 1:nspins), &
616  drho_s(1:4, 1:na, 1:nr, 1:nspins), drho1_s(1:4, 1:na, 1:nr, 1:nspins))
617  ALLOCATE (vxg_h(1:3, 1:na, 1:nr, 1:nspins), vxg_s(1:3, 1:na, 1:nr, 1:nspins))
618  ELSE
619  ALLOCATE (drho_h(1, 1, 1, 1), drho1_h(1, 1, 1, 1), &
620  drho_s(1, 1, 1, 1), drho1_s(1, 1, 1, 1))
621  ALLOCATE (vxg_h(1, 1, 1, 1), vxg_s(1, 1, 1, 1))
622  rrho = 0.0_dp
623  END IF
624  vxg_h = 0.0_dp
625  vxg_s = 0.0_dp
626 
627  ! parallelization
628  local_loop_limit = get_limit(natom, para_env%num_pe, para_env%mepos)
629 
630  DO iatom = local_loop_limit(1), local_loop_limit(2) !1,natom
631  atom = atom_list(iatom)
632 
633  rho_atom_set(atom)%exc_h = 0.0_dp
634  rho_atom_set(atom)%exc_s = 0.0_dp
635  rho1_atom_set(atom)%exc_h = 0.0_dp
636  rho1_atom_set(atom)%exc_s = 0.0_dp
637 
638  rho_atom => rho_atom_set(atom)
639  rho1_atom => rho1_atom_set(atom)
640  NULLIFY (r_h, r_s, dr_h, dr_s, r_h_d, r_s_d)
641  NULLIFY (r1_h, r1_s, dr1_h, dr1_s, r1_h_d, r1_s_d)
642  rho_h = 0.0_dp
643  rho_s = 0.0_dp
644  rho1_h = 0.0_dp
645  rho1_s = 0.0_dp
646  IF (gradient_functional) THEN
647  CALL get_rho_atom(rho_atom=rho_atom, &
648  rho_rad_h=r_h, rho_rad_s=r_s, &
649  drho_rad_h=dr_h, drho_rad_s=dr_s, &
650  rho_rad_h_d=r_h_d, rho_rad_s_d=r_s_d)
651  CALL get_rho_atom(rho_atom=rho1_atom, &
652  rho_rad_h=r1_h, rho_rad_s=r1_s, &
653  drho_rad_h=dr1_h, drho_rad_s=dr1_s, &
654  rho_rad_h_d=r1_h_d, rho_rad_s_d=r1_s_d)
655  drho_h = 0.0_dp; drho_s = 0.0_dp
656  drho1_h = 0.0_dp; drho1_s = 0.0_dp
657  ELSE
658  CALL get_rho_atom(rho_atom=rho_atom, &
659  rho_rad_h=r_h, rho_rad_s=r_s)
660  CALL get_rho_atom(rho_atom=rho1_atom, &
661  rho_rad_h=r1_h, rho_rad_s=r1_s)
662  END IF
663 
664  rtot = 0.0_dp
665 
666  DO ir = 1, nr
667  CALL calc_rho_angular(grid_atom, harmonics, nspins, gradient_functional, &
668  ir, r_h, r_s, rho_h, rho_s, dr_h, dr_s, r_h_d, r_s_d, &
669  drho_h, drho_s)
670  CALL calc_rho_angular(grid_atom, harmonics, nspins, gradient_functional, &
671  ir, r1_h, r1_s, rho1_h, rho1_s, dr1_h, dr1_s, r1_h_d, r1_s_d, &
672  drho1_h, drho1_s)
673  END DO
674  IF (scale_rho) THEN
675  rho_h = 2.0_dp*rho_h
676  rho_s = 2.0_dp*rho_s
677  IF (gradient_functional) THEN
678  drho_h = 2.0_dp*drho_h
679  drho_s = 2.0_dp*drho_s
680  END IF
681  END IF
682 
683  DO ir = 1, nr
684  IF (tau_f) THEN
685  CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, drho_h, tau_h, na, ir)
686  CALL fill_rho_set(rho1_set_h, lsd, nspins, needs, rho1_h, drho1_h, tau1_h, na, ir)
687  CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, drho_s, tau_s, na, ir)
688  CALL fill_rho_set(rho1_set_s, lsd, nspins, needs, rho1_s, drho1_s, tau1_s, na, ir)
689  ELSE IF (gradient_functional) THEN
690  CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, drho_h, rtau, na, ir)
691  CALL fill_rho_set(rho1_set_h, lsd, nspins, needs, rho1_h, drho1_h, rtau, na, ir)
692  CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, drho_s, rtau, na, ir)
693  CALL fill_rho_set(rho1_set_s, lsd, nspins, needs, rho1_s, drho1_s, rtau, na, ir)
694  ELSE
695  CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, rrho, rtau, na, ir)
696  CALL fill_rho_set(rho1_set_h, lsd, nspins, needs, rho1_h, rrho, rtau, na, ir)
697  CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, rrho, rtau, na, ir)
698  CALL fill_rho_set(rho1_set_s, lsd, nspins, needs, rho1_s, rrho, rtau, na, ir)
699  END IF
700  END DO
701 
702  CALL xc_2nd_deriv_of_r(xc_section=xc_section, &
703  rho_set=rho_set_h, rho1_set=rho1_set_h, &
704  deriv_set=deriv_set, &
705  w=weight, vxc=vxc_h, vxg=vxg_h, do_triplet=do_triplet)
706  CALL xc_2nd_deriv_of_r(xc_section=xc_section, &
707  rho_set=rho_set_s, rho1_set=rho1_set_s, &
708  deriv_set=deriv_set, &
709  w=weight, vxc=vxc_s, vxg=vxg_s, do_triplet=do_triplet)
710 
711  CALL get_rho_atom(rho_atom=rho1_atom, ga_vlocal_gb_h=int_hh, ga_vlocal_gb_s=int_ss)
712  IF (gradient_functional) THEN
713  CALL gavxcgb_gc(vxc_h, vxc_s, vxg_h, vxg_s, int_hh, int_ss, &
714  grid_atom, basis_1c, harmonics, nspins)
715  ELSE
716  CALL gavxcgb_nogc(vxc_h, vxc_s, int_hh, int_ss, &
717  grid_atom, basis_1c, harmonics, nspins)
718  END IF
719 
720  NULLIFY (r_h, r_s, dr_h, dr_s)
721 
722  END DO
723 
724  ! some cleanup
725  NULLIFY (weight)
726  DEALLOCATE (rho_h, rho_s, rho1_h, rho1_s, vxc_h, vxc_s)
727  DEALLOCATE (drho_h, drho_s, vxg_h, vxg_s)
728  DEALLOCATE (drho1_h, drho1_s)
729 
730  CALL xc_dset_release(deriv_set)
731  CALL xc_rho_set_release(rho_set_h)
732  CALL xc_rho_set_release(rho1_set_h)
733  CALL xc_rho_set_release(rho_set_s)
734  CALL xc_rho_set_release(rho1_set_s)
735 
736  END DO
737 
738  CALL timestop(handle)
739 
740  END SUBROUTINE calculate_xc_2nd_deriv_atom
741 
742 ! **************************************************************************************************
743 !> \brief ...
744 !> \param qs_env ...
745 !> \param rho0_atom_set ...
746 !> \param rho1_atom_set ...
747 !> \param rho2_atom_set ...
748 !> \param kind_set ...
749 !> \param xc_section ...
750 !> \param is_triplet ...
751 !> \param accuracy ...
752 ! **************************************************************************************************
753  SUBROUTINE calculate_gfxc_atom(qs_env, rho0_atom_set, rho1_atom_set, rho2_atom_set, &
754  kind_set, xc_section, is_triplet, accuracy)
755 
756  TYPE(qs_environment_type), POINTER :: qs_env
757  TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho0_atom_set, rho1_atom_set, &
758  rho2_atom_set
759  TYPE(qs_kind_type), DIMENSION(:), POINTER :: kind_set
760  TYPE(section_vals_type), OPTIONAL, POINTER :: xc_section
761  LOGICAL, INTENT(IN) :: is_triplet
762  INTEGER, INTENT(IN) :: accuracy
763 
764  CHARACTER(LEN=*), PARAMETER :: routinen = 'calculate_gfxc_atom'
765  REAL(kind=dp), PARAMETER :: epsrho = 5.e-4_dp
766 
767  INTEGER :: bo(2), handle, iat, iatom, ikind, ir, &
768  istep, mspins, myfun, na, natom, nf, &
769  nr, ns, nspins, nstep, num_pe
770  INTEGER, DIMENSION(2, 3) :: bounds
771  INTEGER, DIMENSION(:), POINTER :: atom_list
772  LOGICAL :: donlcc, gradient_f, lsd, nlcc, paw_atom, &
773  tau_f
774  REAL(dp) :: beta, density_cut, exc_h, exc_s, &
775  gradient_cut, oeps1, oeps2, tau_cut
776  REAL(dp), DIMENSION(1, 1, 1) :: tau_d
777  REAL(dp), DIMENSION(1, 1, 1, 1) :: rho_d
778  REAL(dp), DIMENSION(:, :), POINTER :: rho_nlcc, weight
779  REAL(dp), DIMENSION(:, :, :), POINTER :: rho0_h, rho0_s, rho1_h, rho1_s, rho_h, &
780  rho_s, tau0_h, tau0_s, tau1_h, tau1_s, &
781  tau_h, tau_s, vtau_h, vtau_s, vxc_h, &
782  vxc_s
783  REAL(dp), DIMENSION(:, :, :, :), POINTER :: drho0_h, drho0_s, drho1_h, drho1_s, &
784  drho_h, drho_s, vxg_h, vxg_s
785  REAL(kind=dp), DIMENSION(-4:4) :: ak, bl
786  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
787  TYPE(dft_control_type), POINTER :: dft_control
788  TYPE(grid_atom_type), POINTER :: grid_atom
789  TYPE(gto_basis_set_type), POINTER :: basis_1c
790  TYPE(harmonics_atom_type), POINTER :: harmonics
791  TYPE(mp_para_env_type), POINTER :: para_env
792  TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: dr_h, dr_s, fint_hh, fint_ss, int_hh, &
793  int_ss, r_h, r_s
794  TYPE(rho_atom_coeff), DIMENSION(:, :), POINTER :: r_h_d, r_s_d
795  TYPE(rho_atom_type), POINTER :: rho0_atom, rho1_atom, rho2_atom
796  TYPE(section_vals_type), POINTER :: xc_fun_section
797  TYPE(xc_derivative_set_type) :: deriv_set
798  TYPE(xc_rho_cflags_type) :: needs
799  TYPE(xc_rho_set_type) :: rho_set_h, rho_set_s
800 
801  CALL timeset(routinen, handle)
802 
803  ak = 0.0_dp
804  bl = 0.0_dp
805  SELECT CASE (accuracy)
806  CASE (:4)
807  nstep = 2
808  ak(-2:2) = (/1.0_dp, -8.0_dp, 0.0_dp, 8.0_dp, -1.0_dp/)/12.0_dp
809  bl(-2:2) = (/-1.0_dp, 16.0_dp, -30.0_dp, 16.0_dp, -1.0_dp/)/12.0_dp
810  CASE (5:7)
811  nstep = 3
812  ak(-3:3) = (/-1.0_dp, 9.0_dp, -45.0_dp, 0.0_dp, 45.0_dp, -9.0_dp, 1.0_dp/)/60.0_dp
813  bl(-3:3) = (/2.0_dp, -27.0_dp, 270.0_dp, -490.0_dp, 270.0_dp, -27.0_dp, 2.0_dp/)/180.0_dp
814  CASE (8:)
815  nstep = 4
816  ak(-4:4) = (/1.0_dp, -32.0_dp/3.0_dp, 56.0_dp, -224.0_dp, 0.0_dp, &
817  224.0_dp, -56.0_dp, 32.0_dp/3.0_dp, -1.0_dp/)/280.0_dp
818  bl(-4:4) = (/-1.0_dp, 128.0_dp/9.0_dp, -112.0_dp, 896.0_dp, -14350.0_dp/9.0_dp, &
819  896.0_dp, -112.0_dp, 128.0_dp/9.0_dp, -1.0_dp/)/560.0_dp
820  END SELECT
821  oeps1 = 1.0_dp/epsrho
822  oeps2 = 1.0_dp/(epsrho**2)
823 
824  CALL get_qs_env(qs_env=qs_env, &
825  dft_control=dft_control, &
826  para_env=para_env, &
827  atomic_kind_set=atomic_kind_set)
828 
829  xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
830  CALL section_vals_val_get(xc_fun_section, "_SECTION_PARAMETERS_", i_val=myfun)
831 
832  IF (myfun == xc_none) THEN
833  ! no action needed?
834  ELSE
835  CALL section_vals_val_get(xc_section, "DENSITY_CUTOFF", r_val=density_cut)
836  CALL section_vals_val_get(xc_section, "GRADIENT_CUTOFF", r_val=gradient_cut)
837  CALL section_vals_val_get(xc_section, "TAU_CUTOFF", r_val=tau_cut)
838 
839  nlcc = has_nlcc(kind_set)
840  lsd = dft_control%lsd
841  nspins = dft_control%nspins
842  mspins = nspins
843  IF (is_triplet) THEN
844  cpassert(nspins == 1)
845  lsd = .true.
846  mspins = 2
847  END IF
848  needs = xc_functionals_get_needs(xc_fun_section, lsd=lsd, calc_potential=.true.)
849  gradient_f = (needs%drho .OR. needs%drho_spin)
850  tau_f = (needs%tau .OR. needs%tau_spin)
851 
852  ! Here starts the loop over all the atoms
853  DO ikind = 1, SIZE(atomic_kind_set)
854  CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
855  CALL get_qs_kind(kind_set(ikind), paw_atom=paw_atom, &
856  harmonics=harmonics, grid_atom=grid_atom)
857  CALL get_qs_kind(kind_set(ikind), basis_set=basis_1c, basis_type="GAPW_1C")
858 
859  IF (.NOT. paw_atom) cycle
860 
861  nr = grid_atom%nr
862  na = grid_atom%ng_sphere
863 
864  ! Prepare the structures needed to calculate and store the xc derivatives
865 
866  ! Array dimension: here anly one dimensional arrays are used,
867  ! i.e. only the first column of deriv_data is read.
868  ! The other to dimensions are set to size equal 1
869  bounds(1:2, 1:3) = 1
870  bounds(2, 1) = na
871  bounds(2, 2) = nr
872 
873  ! create a place where to put the derivatives
874  CALL xc_dset_create(deriv_set, local_bounds=bounds)
875  ! create the place where to store the argument for the functionals
876  CALL xc_rho_set_create(rho_set_h, bounds, rho_cutoff=density_cut, &
877  drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
878  CALL xc_rho_set_create(rho_set_s, bounds, rho_cutoff=density_cut, &
879  drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
880 
881  ! allocate the required 3d arrays where to store rho and drho
882  CALL xc_rho_set_atom_update(rho_set_h, needs, mspins, bounds)
883  CALL xc_rho_set_atom_update(rho_set_s, needs, mspins, bounds)
884 
885  weight => grid_atom%weight
886 
887  ALLOCATE (rho_h(na, nr, mspins), rho_s(na, nr, mspins), &
888  rho0_h(na, nr, nspins), rho0_s(na, nr, nspins), &
889  rho1_h(na, nr, nspins), rho1_s(na, nr, nspins))
890  ALLOCATE (vxc_h(na, nr, mspins), vxc_s(na, nr, mspins))
891  IF (gradient_f) THEN
892  ALLOCATE (drho_h(4, na, nr, mspins), drho_s(4, na, nr, mspins), &
893  drho0_h(4, na, nr, nspins), drho0_s(4, na, nr, nspins), &
894  drho1_h(4, na, nr, nspins), drho1_s(4, na, nr, nspins))
895  ALLOCATE (vxg_h(3, na, nr, mspins), vxg_s(3, na, nr, mspins))
896  END IF
897  IF (tau_f) THEN
898  ALLOCATE (tau_h(na, nr, mspins), tau_s(na, nr, mspins), &
899  tau0_h(na, nr, nspins), tau0_s(na, nr, nspins), &
900  tau1_h(na, nr, nspins), tau1_s(na, nr, nspins))
901  ALLOCATE (vtau_h(na, nr, mspins), vtau_s(na, nr, mspins))
902  END IF
903  !
904  ! NLCC: prepare rho and drho of the core charge for this KIND
905  donlcc = .false.
906  IF (nlcc) THEN
907  NULLIFY (rho_nlcc)
908  rho_nlcc => kind_set(ikind)%nlcc_pot
909  IF (ASSOCIATED(rho_nlcc)) donlcc = .true.
910  END IF
911 
912  ! Distribute the atoms of this kind
913  num_pe = para_env%num_pe
914  bo = get_limit(natom, num_pe, para_env%mepos)
915 
916  DO iat = bo(1), bo(2)
917  iatom = atom_list(iat)
918  !
919  NULLIFY (int_hh, int_ss)
920  rho0_atom => rho0_atom_set(iatom)
921  CALL get_rho_atom(rho_atom=rho0_atom, ga_vlocal_gb_h=int_hh, ga_vlocal_gb_s=int_ss)
922  ALLOCATE (fint_ss(nspins), fint_hh(nspins))
923  DO ns = 1, nspins
924  nf = SIZE(int_ss(ns)%r_coef, 1)
925  ALLOCATE (fint_ss(ns)%r_coef(nf, nf))
926  nf = SIZE(int_hh(ns)%r_coef, 1)
927  ALLOCATE (fint_hh(ns)%r_coef(nf, nf))
928  END DO
929 
930  ! RHO0
931  rho0_h = 0.0_dp
932  rho0_s = 0.0_dp
933  rho0_atom => rho0_atom_set(iatom)
934  IF (gradient_f) THEN
935  NULLIFY (r_h, r_s, dr_h, dr_s, r_h_d, r_s_d)
936  CALL get_rho_atom(rho_atom=rho0_atom, rho_rad_h=r_h, rho_rad_s=r_s, drho_rad_h=dr_h, &
937  drho_rad_s=dr_s, rho_rad_h_d=r_h_d, rho_rad_s_d=r_s_d)
938  drho0_h = 0.0_dp
939  drho0_s = 0.0_dp
940  ELSE
941  NULLIFY (r_h, r_s)
942  CALL get_rho_atom(rho_atom=rho0_atom, rho_rad_h=r_h, rho_rad_s=r_s)
943  rho_d = 0.0_dp
944  END IF
945  DO ir = 1, nr
946  CALL calc_rho_angular(grid_atom, harmonics, nspins, gradient_f, &
947  ir, r_h, r_s, rho0_h, rho0_s, dr_h, dr_s, &
948  r_h_d, r_s_d, drho0_h, drho0_s)
949  IF (donlcc) THEN
950  CALL calc_rho_nlcc(grid_atom, nspins, gradient_f, &
951  ir, rho_nlcc(:, 1), rho0_h, rho0_s, rho_nlcc(:, 2), drho0_h, drho0_s)
952  END IF
953  END DO
954  IF (tau_f) THEN
955  !compute tau on the grid all at once
956  CALL calc_tau_atom(tau0_h, tau0_s, rho0_atom, kind_set(ikind), nspins)
957  ELSE
958  tau_d = 0.0_dp
959  END IF
960  ! RHO1
961  rho1_h = 0.0_dp
962  rho1_s = 0.0_dp
963  rho1_atom => rho1_atom_set(iatom)
964  IF (gradient_f) THEN
965  NULLIFY (r_h, r_s, dr_h, dr_s, r_h_d, r_s_d)
966  CALL get_rho_atom(rho_atom=rho1_atom, rho_rad_h=r_h, rho_rad_s=r_s, drho_rad_h=dr_h, &
967  drho_rad_s=dr_s, rho_rad_h_d=r_h_d, rho_rad_s_d=r_s_d)
968  drho1_h = 0.0_dp
969  drho1_s = 0.0_dp
970  ELSE
971  NULLIFY (r_h, r_s)
972  CALL get_rho_atom(rho_atom=rho1_atom, rho_rad_h=r_h, rho_rad_s=r_s)
973  END IF
974  DO ir = 1, nr
975  CALL calc_rho_angular(grid_atom, harmonics, nspins, gradient_f, &
976  ir, r_h, r_s, rho1_h, rho1_s, dr_h, dr_s, &
977  r_h_d, r_s_d, drho1_h, drho1_s)
978  END DO
979  IF (tau_f) THEN
980  !compute tau on the grid all at once
981  CALL calc_tau_atom(tau1_h, tau1_s, rho1_atom, kind_set(ikind), nspins)
982  END IF
983  ! RHO2
984  rho2_atom => rho2_atom_set(iatom)
985 
986  DO istep = -nstep, nstep
987 
988  beta = real(istep, kind=dp)*epsrho
989 
990  IF (is_triplet) THEN
991  rho_h(:, :, 1) = rho0_h(:, :, 1) + beta*rho1_h(:, :, 1)
992  rho_h(:, :, 2) = rho0_h(:, :, 1)
993  rho_h = 0.5_dp*rho_h
994  rho_s(:, :, 1) = rho0_s(:, :, 1) + beta*rho1_s(:, :, 1)
995  rho_s(:, :, 2) = rho0_s(:, :, 1)
996  rho_s = 0.5_dp*rho_s
997  IF (gradient_f) THEN
998  drho_h(:, :, :, 1) = drho0_h(:, :, :, 1) + beta*drho1_h(:, :, :, 1)
999  drho_h(:, :, :, 2) = drho0_h(:, :, :, 1)
1000  drho_h = 0.5_dp*drho_h
1001  drho_s(:, :, :, 1) = drho0_s(:, :, :, 1) + beta*drho1_s(:, :, :, 1)
1002  drho_s(:, :, :, 2) = drho0_s(:, :, :, 1)
1003  drho_s = 0.5_dp*drho_s
1004  END IF
1005  IF (tau_f) THEN
1006  tau_h(:, :, 1) = tau0_h(:, :, 1) + beta*tau1_h(:, :, 1)
1007  tau_h(:, :, 2) = tau0_h(:, :, 1)
1008  tau_h = 0.5_dp*tau0_h
1009  tau_s(:, :, 1) = tau0_s(:, :, 1) + beta*tau1_s(:, :, 1)
1010  tau_s(:, :, 2) = tau0_s(:, :, 1)
1011  tau_s = 0.5_dp*tau0_s
1012  END IF
1013  ELSE
1014  rho_h = rho0_h + beta*rho1_h
1015  rho_s = rho0_s + beta*rho1_s
1016  IF (gradient_f) THEN
1017  drho_h = drho0_h + beta*drho1_h
1018  drho_s = drho0_s + beta*drho1_s
1019  END IF
1020  IF (tau_f) THEN
1021  tau_h = tau0_h + beta*tau1_h
1022  tau_s = tau0_s + beta*tau1_s
1023  END IF
1024  END IF
1025  !
1026  IF (gradient_f) THEN
1027  drho_h(4, :, :, :) = sqrt( &
1028  drho_h(1, :, :, :)*drho_h(1, :, :, :) + &
1029  drho_h(2, :, :, :)*drho_h(2, :, :, :) + &
1030  drho_h(3, :, :, :)*drho_h(3, :, :, :))
1031 
1032  drho_s(4, :, :, :) = sqrt( &
1033  drho_s(1, :, :, :)*drho_s(1, :, :, :) + &
1034  drho_s(2, :, :, :)*drho_s(2, :, :, :) + &
1035  drho_s(3, :, :, :)*drho_s(3, :, :, :))
1036  END IF
1037 
1038  DO ir = 1, nr
1039  IF (tau_f) THEN
1040  CALL fill_rho_set(rho_set_h, lsd, mspins, needs, rho_h, drho_h, tau_h, na, ir)
1041  CALL fill_rho_set(rho_set_s, lsd, mspins, needs, rho_s, drho_s, tau_s, na, ir)
1042  ELSE IF (gradient_f) THEN
1043  CALL fill_rho_set(rho_set_h, lsd, mspins, needs, rho_h, drho_h, tau_d, na, ir)
1044  CALL fill_rho_set(rho_set_s, lsd, mspins, needs, rho_s, drho_s, tau_d, na, ir)
1045  ELSE
1046  CALL fill_rho_set(rho_set_h, lsd, mspins, needs, rho_h, rho_d, tau_d, na, ir)
1047  CALL fill_rho_set(rho_set_s, lsd, mspins, needs, rho_s, rho_d, tau_d, na, ir)
1048  END IF
1049  END DO
1050 
1051  ! hard atom density !
1052  CALL xc_dset_zero_all(deriv_set)
1053  CALL vxc_of_r_new(xc_fun_section, rho_set_h, deriv_set, 1, needs, weight, &
1054  lsd, na, nr, exc_h, vxc_h, vxg_h, vtau_h)
1055  IF (is_triplet) THEN
1056  vxc_h(:, :, 1) = vxc_h(:, :, 1) - vxc_h(:, :, 2)
1057  IF (gradient_f) THEN
1058  vxg_h(:, :, :, 1) = vxg_h(:, :, :, 1) - vxg_h(:, :, :, 2)
1059  END IF
1060  IF (tau_f) THEN
1061  vtau_h(:, :, 1) = vtau_h(:, :, 1) - vtau_h(:, :, 2)
1062  END IF
1063  END IF
1064  ! soft atom density !
1065  CALL xc_dset_zero_all(deriv_set)
1066  CALL vxc_of_r_new(xc_fun_section, rho_set_s, deriv_set, 1, needs, weight, &
1067  lsd, na, nr, exc_s, vxc_s, vxg_s, vtau_s)
1068  IF (is_triplet) THEN
1069  vxc_s(:, :, 1) = vxc_s(:, :, 1) - vxc_s(:, :, 2)
1070  IF (gradient_f) THEN
1071  vxg_s(:, :, :, 1) = vxg_s(:, :, :, 1) - vxg_s(:, :, :, 2)
1072  END IF
1073  IF (tau_f) THEN
1074  vtau_s(:, :, 1) = vtau_s(:, :, 1) - vtau_s(:, :, 2)
1075  END IF
1076  END IF
1077  ! potentials
1078  DO ns = 1, nspins
1079  fint_hh(ns)%r_coef(:, :) = 0.0_dp
1080  fint_ss(ns)%r_coef(:, :) = 0.0_dp
1081  END DO
1082  IF (gradient_f) THEN
1083  CALL gavxcgb_gc(vxc_h, vxc_s, vxg_h, vxg_s, fint_hh, fint_ss, &
1084  grid_atom, basis_1c, harmonics, nspins)
1085  ELSE
1086  CALL gavxcgb_nogc(vxc_h, vxc_s, fint_hh, fint_ss, &
1087  grid_atom, basis_1c, harmonics, nspins)
1088  END IF
1089  IF (tau_f) THEN
1090  CALL dgavtaudgb(vtau_h, vtau_s, fint_hh, fint_ss, &
1091  grid_atom, basis_1c, harmonics, nspins)
1092  END IF
1093  ! first derivative fxc
1094  NULLIFY (int_hh, int_ss)
1095  CALL get_rho_atom(rho_atom=rho1_atom, ga_vlocal_gb_h=int_hh, ga_vlocal_gb_s=int_ss)
1096  DO ns = 1, nspins
1097  int_ss(ns)%r_coef(:, :) = int_ss(ns)%r_coef(:, :) + oeps1*ak(istep)*fint_ss(ns)%r_coef(:, :)
1098  int_hh(ns)%r_coef(:, :) = int_hh(ns)%r_coef(:, :) + oeps1*ak(istep)*fint_hh(ns)%r_coef(:, :)
1099  END DO
1100  ! second derivative gxc
1101  NULLIFY (int_hh, int_ss)
1102  CALL get_rho_atom(rho_atom=rho2_atom, ga_vlocal_gb_h=int_hh, ga_vlocal_gb_s=int_ss)
1103  DO ns = 1, nspins
1104  int_ss(ns)%r_coef(:, :) = int_ss(ns)%r_coef(:, :) + oeps2*bl(istep)*fint_ss(ns)%r_coef(:, :)
1105  int_hh(ns)%r_coef(:, :) = int_hh(ns)%r_coef(:, :) + oeps2*bl(istep)*fint_hh(ns)%r_coef(:, :)
1106  END DO
1107  END DO
1108  !
1109  DO ns = 1, nspins
1110  DEALLOCATE (fint_ss(ns)%r_coef)
1111  DEALLOCATE (fint_hh(ns)%r_coef)
1112  END DO
1113  DEALLOCATE (fint_ss, fint_hh)
1114 
1115  END DO ! iat
1116 
1117  ! Release the xc structure used to store the xc derivatives
1118  CALL xc_dset_release(deriv_set)
1119  CALL xc_rho_set_release(rho_set_h)
1120  CALL xc_rho_set_release(rho_set_s)
1121 
1122  DEALLOCATE (rho_h, rho_s, rho0_h, rho0_s, rho1_h, rho1_s)
1123  DEALLOCATE (vxc_h, vxc_s)
1124  IF (gradient_f) THEN
1125  DEALLOCATE (drho_h, drho_s, drho0_h, drho0_s, drho1_h, drho1_s)
1126  DEALLOCATE (vxg_h, vxg_s)
1127  END IF
1128  IF (tau_f) THEN
1129  DEALLOCATE (tau_h, tau_s, tau0_h, tau0_s, tau1_h, tau1_s)
1130  DEALLOCATE (vtau_h, vtau_s)
1131  END IF
1132 
1133  END DO ! ikind
1134 
1135  END IF !xc_none
1136 
1137  CALL timestop(handle)
1138 
1139  END SUBROUTINE calculate_gfxc_atom
1140 
1141 ! **************************************************************************************************
1142 !> \brief ...
1143 !> \param qs_env ...
1144 !> \param rho0_atom_set ...
1145 !> \param rho1_atom_set ...
1146 !> \param rho2_atom_set ...
1147 !> \param kind_set ...
1148 !> \param xc_section ...
1149 !> \param is_triplet ...
1150 !> \param accuracy ...
1151 !> \param epsrho ...
1152 ! **************************************************************************************************
1153  SUBROUTINE gfxc_atom_diff(qs_env, rho0_atom_set, rho1_atom_set, rho2_atom_set, &
1154  kind_set, xc_section, is_triplet, accuracy, epsrho)
1155 
1156  TYPE(qs_environment_type), POINTER :: qs_env
1157  TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho0_atom_set, rho1_atom_set, &
1158  rho2_atom_set
1159  TYPE(qs_kind_type), DIMENSION(:), POINTER :: kind_set
1160  TYPE(section_vals_type), OPTIONAL, POINTER :: xc_section
1161  LOGICAL, INTENT(IN) :: is_triplet
1162  INTEGER, INTENT(IN) :: accuracy
1163  REAL(kind=dp), INTENT(IN) :: epsrho
1164 
1165  CHARACTER(LEN=*), PARAMETER :: routinen = 'gfxc_atom_diff'
1166 
1167  INTEGER :: bo(2), handle, iat, iatom, ikind, ir, &
1168  istep, mspins, myfun, na, natom, nf, &
1169  nr, ns, nspins, nstep, num_pe
1170  INTEGER, DIMENSION(2, 3) :: bounds
1171  INTEGER, DIMENSION(:), POINTER :: atom_list
1172  LOGICAL :: donlcc, gradient_f, lsd, nlcc, paw_atom, &
1173  tau_f
1174  REAL(dp) :: beta, density_cut, gradient_cut, oeps1, &
1175  tau_cut
1176  REAL(dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: vxc_h, vxc_s
1177  REAL(dp), DIMENSION(1, 1, 1) :: tau_d
1178  REAL(dp), DIMENSION(1, 1, 1, 1) :: rho_d
1179  REAL(dp), DIMENSION(:, :), POINTER :: rho_nlcc, weight
1180  REAL(dp), DIMENSION(:, :, :), POINTER :: rho0_h, rho0_s, rho1_h, rho1_s, rho_h, &
1181  rho_s, tau0_h, tau0_s, tau1_h, tau1_s, &
1182  tau_h, tau_s, vtau_h, vtau_s
1183  REAL(dp), DIMENSION(:, :, :, :), POINTER :: drho0_h, drho0_s, drho1_h, drho1_s, &
1184  drho_h, drho_s, vxg_h, vxg_s
1185  REAL(kind=dp), DIMENSION(-4:4) :: ak
1186  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1187  TYPE(dft_control_type), POINTER :: dft_control
1188  TYPE(grid_atom_type), POINTER :: grid_atom
1189  TYPE(gto_basis_set_type), POINTER :: basis_1c
1190  TYPE(harmonics_atom_type), POINTER :: harmonics
1191  TYPE(mp_para_env_type), POINTER :: para_env
1192  TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: dr_h, dr_s, fint_hh, fint_ss, int_hh, &
1193  int_ss, r_h, r_s
1194  TYPE(rho_atom_coeff), DIMENSION(:, :), POINTER :: r_h_d, r_s_d
1195  TYPE(rho_atom_type), POINTER :: rho0_atom, rho1_atom, rho2_atom
1196  TYPE(section_vals_type), POINTER :: xc_fun_section
1197  TYPE(xc_derivative_set_type) :: deriv_set
1198  TYPE(xc_rho_cflags_type) :: needs
1199  TYPE(xc_rho_set_type) :: rho1_set_h, rho1_set_s, rho_set_h, &
1200  rho_set_s
1201 
1202  CALL timeset(routinen, handle)
1203 
1204  ak = 0.0_dp
1205  SELECT CASE (accuracy)
1206  CASE (:4)
1207  nstep = 2
1208  ak(-2:2) = (/1.0_dp, -8.0_dp, 0.0_dp, 8.0_dp, -1.0_dp/)/12.0_dp
1209  CASE (5:7)
1210  nstep = 3
1211  ak(-3:3) = (/-1.0_dp, 9.0_dp, -45.0_dp, 0.0_dp, 45.0_dp, -9.0_dp, 1.0_dp/)/60.0_dp
1212  CASE (8:)
1213  nstep = 4
1214  ak(-4:4) = (/1.0_dp, -32.0_dp/3.0_dp, 56.0_dp, -224.0_dp, 0.0_dp, &
1215  224.0_dp, -56.0_dp, 32.0_dp/3.0_dp, -1.0_dp/)/280.0_dp
1216  END SELECT
1217  oeps1 = 1.0_dp/epsrho
1218 
1219  CALL get_qs_env(qs_env=qs_env, &
1220  dft_control=dft_control, &
1221  para_env=para_env, &
1222  atomic_kind_set=atomic_kind_set)
1223 
1224  xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
1225  CALL section_vals_val_get(xc_fun_section, "_SECTION_PARAMETERS_", i_val=myfun)
1226 
1227  IF (myfun == xc_none) THEN
1228  ! no action needed?
1229  ELSE
1230  ! calculate fxc
1231  CALL calculate_xc_2nd_deriv_atom(rho0_atom_set, rho1_atom_set, qs_env, xc_section, para_env, &
1232  do_triplet=is_triplet, kind_set_external=kind_set)
1233 
1234  CALL section_vals_val_get(xc_section, "DENSITY_CUTOFF", r_val=density_cut)
1235  CALL section_vals_val_get(xc_section, "GRADIENT_CUTOFF", r_val=gradient_cut)
1236  CALL section_vals_val_get(xc_section, "TAU_CUTOFF", r_val=tau_cut)
1237 
1238  nlcc = has_nlcc(kind_set)
1239  lsd = dft_control%lsd
1240  nspins = dft_control%nspins
1241  mspins = nspins
1242  IF (is_triplet) THEN
1243  cpassert(nspins == 1)
1244  lsd = .true.
1245  mspins = 2
1246  END IF
1247  needs = xc_functionals_get_needs(xc_fun_section, lsd=lsd, calc_potential=.true.)
1248  gradient_f = (needs%drho .OR. needs%drho_spin)
1249  tau_f = (needs%tau .OR. needs%tau_spin)
1250 
1251  ! Here starts the loop over all the atoms
1252  DO ikind = 1, SIZE(atomic_kind_set)
1253  CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
1254  CALL get_qs_kind(kind_set(ikind), paw_atom=paw_atom, &
1255  harmonics=harmonics, grid_atom=grid_atom)
1256  CALL get_qs_kind(kind_set(ikind), basis_set=basis_1c, basis_type="GAPW_1C")
1257 
1258  IF (.NOT. paw_atom) cycle
1259 
1260  nr = grid_atom%nr
1261  na = grid_atom%ng_sphere
1262 
1263  ! Prepare the structures needed to calculate and store the xc derivatives
1264 
1265  ! Array dimension: here anly one dimensional arrays are used,
1266  ! i.e. only the first column of deriv_data is read.
1267  ! The other to dimensions are set to size equal 1
1268  bounds(1:2, 1:3) = 1
1269  bounds(2, 1) = na
1270  bounds(2, 2) = nr
1271 
1272  ! create a place where to put the derivatives
1273  CALL xc_dset_create(deriv_set, local_bounds=bounds)
1274  ! create the place where to store the argument for the functionals
1275  CALL xc_rho_set_create(rho_set_h, bounds, rho_cutoff=density_cut, &
1276  drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
1277  CALL xc_rho_set_create(rho_set_s, bounds, rho_cutoff=density_cut, &
1278  drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
1279  CALL xc_rho_set_create(rho1_set_h, bounds, rho_cutoff=density_cut, &
1280  drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
1281  CALL xc_rho_set_create(rho1_set_s, bounds, rho_cutoff=density_cut, &
1282  drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
1283 
1284  ! allocate the required 3d arrays where to store rho and drho
1285  CALL xc_rho_set_atom_update(rho_set_h, needs, mspins, bounds)
1286  CALL xc_rho_set_atom_update(rho_set_s, needs, mspins, bounds)
1287  CALL xc_rho_set_atom_update(rho1_set_h, needs, mspins, bounds)
1288  CALL xc_rho_set_atom_update(rho1_set_s, needs, mspins, bounds)
1289 
1290  weight => grid_atom%weight
1291 
1292  ALLOCATE (rho_h(na, nr, nspins), rho_s(na, nr, nspins), &
1293  rho0_h(na, nr, nspins), rho0_s(na, nr, nspins), &
1294  rho1_h(na, nr, nspins), rho1_s(na, nr, nspins))
1295  ALLOCATE (vxc_h(na, nr, nspins), vxc_s(na, nr, nspins))
1296  IF (gradient_f) THEN
1297  ALLOCATE (drho_h(4, na, nr, nspins), drho_s(4, na, nr, nspins), &
1298  drho0_h(4, na, nr, nspins), drho0_s(4, na, nr, nspins), &
1299  drho1_h(4, na, nr, nspins), drho1_s(4, na, nr, nspins))
1300  ALLOCATE (vxg_h(3, na, nr, nspins), vxg_s(3, na, nr, nspins))
1301  END IF
1302  IF (tau_f) THEN
1303  ALLOCATE (tau_h(na, nr, nspins), tau_s(na, nr, nspins), &
1304  tau0_h(na, nr, nspins), tau0_s(na, nr, nspins), &
1305  tau1_h(na, nr, nspins), tau1_s(na, nr, nspins))
1306  ALLOCATE (vtau_h(na, nr, nspins), vtau_s(na, nr, nspins))
1307  END IF
1308  !
1309  ! NLCC: prepare rho and drho of the core charge for this KIND
1310  donlcc = .false.
1311  IF (nlcc) THEN
1312  NULLIFY (rho_nlcc)
1313  rho_nlcc => kind_set(ikind)%nlcc_pot
1314  IF (ASSOCIATED(rho_nlcc)) donlcc = .true.
1315  END IF
1316 
1317  ! Distribute the atoms of this kind
1318  num_pe = para_env%num_pe
1319  bo = get_limit(natom, num_pe, para_env%mepos)
1320 
1321  DO iat = bo(1), bo(2)
1322  iatom = atom_list(iat)
1323  !
1324  NULLIFY (int_hh, int_ss)
1325  rho0_atom => rho0_atom_set(iatom)
1326  CALL get_rho_atom(rho_atom=rho0_atom, ga_vlocal_gb_h=int_hh, ga_vlocal_gb_s=int_ss)
1327  ALLOCATE (fint_ss(nspins), fint_hh(nspins))
1328  DO ns = 1, nspins
1329  nf = SIZE(int_ss(ns)%r_coef, 1)
1330  ALLOCATE (fint_ss(ns)%r_coef(nf, nf))
1331  nf = SIZE(int_hh(ns)%r_coef, 1)
1332  ALLOCATE (fint_hh(ns)%r_coef(nf, nf))
1333  END DO
1334 
1335  ! RHO0
1336  rho0_h = 0.0_dp
1337  rho0_s = 0.0_dp
1338  rho0_atom => rho0_atom_set(iatom)
1339  IF (gradient_f) THEN
1340  NULLIFY (r_h, r_s, dr_h, dr_s, r_h_d, r_s_d)
1341  CALL get_rho_atom(rho_atom=rho0_atom, rho_rad_h=r_h, rho_rad_s=r_s, drho_rad_h=dr_h, &
1342  drho_rad_s=dr_s, rho_rad_h_d=r_h_d, rho_rad_s_d=r_s_d)
1343  drho0_h = 0.0_dp
1344  drho0_s = 0.0_dp
1345  ELSE
1346  NULLIFY (r_h, r_s)
1347  CALL get_rho_atom(rho_atom=rho0_atom, rho_rad_h=r_h, rho_rad_s=r_s)
1348  rho_d = 0.0_dp
1349  END IF
1350  DO ir = 1, nr
1351  CALL calc_rho_angular(grid_atom, harmonics, nspins, gradient_f, &
1352  ir, r_h, r_s, rho0_h, rho0_s, dr_h, dr_s, &
1353  r_h_d, r_s_d, drho0_h, drho0_s)
1354  IF (donlcc) THEN
1355  CALL calc_rho_nlcc(grid_atom, nspins, gradient_f, &
1356  ir, rho_nlcc(:, 1), rho0_h, rho0_s, rho_nlcc(:, 2), drho0_h, drho0_s)
1357  END IF
1358  END DO
1359  IF (tau_f) THEN
1360  !compute tau on the grid all at once
1361  CALL calc_tau_atom(tau0_h, tau0_s, rho0_atom, kind_set(ikind), nspins)
1362  ELSE
1363  tau_d = 0.0_dp
1364  END IF
1365  ! RHO1
1366  rho1_h = 0.0_dp
1367  rho1_s = 0.0_dp
1368  rho1_atom => rho1_atom_set(iatom)
1369  IF (gradient_f) THEN
1370  NULLIFY (r_h, r_s, dr_h, dr_s, r_h_d, r_s_d)
1371  CALL get_rho_atom(rho_atom=rho1_atom, rho_rad_h=r_h, rho_rad_s=r_s, drho_rad_h=dr_h, &
1372  drho_rad_s=dr_s, rho_rad_h_d=r_h_d, rho_rad_s_d=r_s_d)
1373  drho1_h = 0.0_dp
1374  drho1_s = 0.0_dp
1375  ELSE
1376  NULLIFY (r_h, r_s)
1377  CALL get_rho_atom(rho_atom=rho1_atom, rho_rad_h=r_h, rho_rad_s=r_s)
1378  END IF
1379  DO ir = 1, nr
1380  CALL calc_rho_angular(grid_atom, harmonics, nspins, gradient_f, &
1381  ir, r_h, r_s, rho1_h, rho1_s, dr_h, dr_s, &
1382  r_h_d, r_s_d, drho1_h, drho1_s)
1383  END DO
1384  IF (tau_f) THEN
1385  !compute tau on the grid all at once
1386  CALL calc_tau_atom(tau1_h, tau1_s, rho1_atom, kind_set(ikind), nspins)
1387  END IF
1388 
1389  DO ir = 1, nr
1390  IF (tau_f) THEN
1391  CALL fill_rho_set(rho1_set_h, lsd, nspins, needs, rho1_h, drho1_h, tau1_h, na, ir)
1392  CALL fill_rho_set(rho1_set_s, lsd, nspins, needs, rho1_s, drho1_s, tau1_s, na, ir)
1393  ELSE IF (gradient_f) THEN
1394  CALL fill_rho_set(rho1_set_h, lsd, nspins, needs, rho1_h, drho1_h, tau_d, na, ir)
1395  CALL fill_rho_set(rho1_set_s, lsd, nspins, needs, rho1_s, drho1_s, tau_d, na, ir)
1396  ELSE
1397  CALL fill_rho_set(rho1_set_h, lsd, nspins, needs, rho1_h, rho_d, tau_d, na, ir)
1398  CALL fill_rho_set(rho1_set_s, lsd, nspins, needs, rho1_s, rho_d, tau_d, na, ir)
1399  END IF
1400  END DO
1401 
1402  ! RHO2
1403  rho2_atom => rho2_atom_set(iatom)
1404 
1405  DO istep = -nstep, nstep
1406 
1407  beta = real(istep, kind=dp)*epsrho
1408 
1409  rho_h = rho0_h + beta*rho1_h
1410  rho_s = rho0_s + beta*rho1_s
1411  IF (gradient_f) THEN
1412  drho_h = drho0_h + beta*drho1_h
1413  drho_s = drho0_s + beta*drho1_s
1414  END IF
1415  IF (tau_f) THEN
1416  tau_h = tau0_h + beta*tau1_h
1417  tau_s = tau0_s + beta*tau1_s
1418  END IF
1419  !
1420  IF (gradient_f) THEN
1421  drho_h(4, :, :, :) = sqrt( &
1422  drho_h(1, :, :, :)*drho_h(1, :, :, :) + &
1423  drho_h(2, :, :, :)*drho_h(2, :, :, :) + &
1424  drho_h(3, :, :, :)*drho_h(3, :, :, :))
1425 
1426  drho_s(4, :, :, :) = sqrt( &
1427  drho_s(1, :, :, :)*drho_s(1, :, :, :) + &
1428  drho_s(2, :, :, :)*drho_s(2, :, :, :) + &
1429  drho_s(3, :, :, :)*drho_s(3, :, :, :))
1430  END IF
1431 
1432  DO ir = 1, nr
1433  IF (tau_f) THEN
1434  CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, drho_h, tau_h, na, ir)
1435  CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, drho_s, tau_s, na, ir)
1436  ELSE IF (gradient_f) THEN
1437  CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, drho_h, tau_d, na, ir)
1438  CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, drho_s, tau_d, na, ir)
1439  ELSE
1440  CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, rho_d, tau_d, na, ir)
1441  CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, rho_d, tau_d, na, ir)
1442  END IF
1443  END DO
1444 
1445  ! hard atom density !
1446  CALL xc_dset_zero_all(deriv_set)
1447  CALL xc_2nd_deriv_of_r(xc_section=xc_section, &
1448  rho_set=rho_set_h, rho1_set=rho1_set_h, &
1449  deriv_set=deriv_set, &
1450  w=weight, vxc=vxc_h, vxg=vxg_h, do_triplet=is_triplet)
1451  ! soft atom density !
1452  CALL xc_dset_zero_all(deriv_set)
1453  CALL xc_2nd_deriv_of_r(xc_section=xc_section, &
1454  rho_set=rho_set_s, rho1_set=rho1_set_s, &
1455  deriv_set=deriv_set, &
1456  w=weight, vxc=vxc_s, vxg=vxg_s, do_triplet=is_triplet)
1457  ! potentials
1458  DO ns = 1, nspins
1459  fint_hh(ns)%r_coef(:, :) = 0.0_dp
1460  fint_ss(ns)%r_coef(:, :) = 0.0_dp
1461  END DO
1462  IF (gradient_f) THEN
1463  CALL gavxcgb_gc(vxc_h, vxc_s, vxg_h, vxg_s, fint_hh, fint_ss, &
1464  grid_atom, basis_1c, harmonics, nspins)
1465  ELSE
1466  CALL gavxcgb_nogc(vxc_h, vxc_s, fint_hh, fint_ss, &
1467  grid_atom, basis_1c, harmonics, nspins)
1468  END IF
1469  IF (tau_f) THEN
1470  CALL dgavtaudgb(vtau_h, vtau_s, fint_hh, fint_ss, &
1471  grid_atom, basis_1c, harmonics, nspins)
1472  END IF
1473  ! second derivative gxc
1474  NULLIFY (int_hh, int_ss)
1475  CALL get_rho_atom(rho_atom=rho2_atom, ga_vlocal_gb_h=int_hh, ga_vlocal_gb_s=int_ss)
1476  DO ns = 1, nspins
1477  int_ss(ns)%r_coef(:, :) = int_ss(ns)%r_coef(:, :) + oeps1*ak(istep)*fint_ss(ns)%r_coef(:, :)
1478  int_hh(ns)%r_coef(:, :) = int_hh(ns)%r_coef(:, :) + oeps1*ak(istep)*fint_hh(ns)%r_coef(:, :)
1479  END DO
1480  END DO
1481  !
1482  DO ns = 1, nspins
1483  DEALLOCATE (fint_ss(ns)%r_coef)
1484  DEALLOCATE (fint_hh(ns)%r_coef)
1485  END DO
1486  DEALLOCATE (fint_ss, fint_hh)
1487 
1488  END DO ! iat
1489 
1490  ! Release the xc structure used to store the xc derivatives
1491  CALL xc_dset_release(deriv_set)
1492  CALL xc_rho_set_release(rho_set_h)
1493  CALL xc_rho_set_release(rho_set_s)
1494  CALL xc_rho_set_release(rho1_set_h)
1495  CALL xc_rho_set_release(rho1_set_s)
1496 
1497  DEALLOCATE (rho_h, rho_s, rho0_h, rho0_s, rho1_h, rho1_s)
1498  DEALLOCATE (vxc_h, vxc_s)
1499  IF (gradient_f) THEN
1500  DEALLOCATE (drho_h, drho_s, drho0_h, drho0_s, drho1_h, drho1_s)
1501  DEALLOCATE (vxg_h, vxg_s)
1502  END IF
1503  IF (tau_f) THEN
1504  DEALLOCATE (tau_h, tau_s, tau0_h, tau0_s, tau1_h, tau1_s)
1505  DEALLOCATE (vtau_h, vtau_s)
1506  END IF
1507 
1508  END DO ! ikind
1509 
1510  END IF !xc_none
1511 
1512  CALL timestop(handle)
1513 
1514  END SUBROUTINE gfxc_atom_diff
1515 
1516 ! **************************************************************************************************
1517 !> \brief ...
1518 !> \param grid_atom ...
1519 !> \param harmonics ...
1520 !> \param nspins ...
1521 !> \param grad_func ...
1522 !> \param ir ...
1523 !> \param r_h ...
1524 !> \param r_s ...
1525 !> \param rho_h ...
1526 !> \param rho_s ...
1527 !> \param dr_h ...
1528 !> \param dr_s ...
1529 !> \param r_h_d ...
1530 !> \param r_s_d ...
1531 !> \param drho_h ...
1532 !> \param drho_s ...
1533 ! **************************************************************************************************
1534  SUBROUTINE calc_rho_angular(grid_atom, harmonics, nspins, grad_func, &
1535  ir, r_h, r_s, rho_h, rho_s, &
1536  dr_h, dr_s, r_h_d, r_s_d, drho_h, drho_s)
1537 
1538  TYPE(grid_atom_type), POINTER :: grid_atom
1539  TYPE(harmonics_atom_type), POINTER :: harmonics
1540  INTEGER, INTENT(IN) :: nspins
1541  LOGICAL, INTENT(IN) :: grad_func
1542  INTEGER, INTENT(IN) :: ir
1543  TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: r_h, r_s
1544  REAL(kind=dp), DIMENSION(:, :, :), POINTER :: rho_h, rho_s
1545  TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: dr_h, dr_s
1546  TYPE(rho_atom_coeff), DIMENSION(:, :), POINTER :: r_h_d, r_s_d
1547  REAL(kind=dp), DIMENSION(:, :, :, :), POINTER :: drho_h, drho_s
1548 
1549  INTEGER :: ia, iso, ispin, na
1550  REAL(kind=dp) :: rad, urad
1551 
1552  cpassert(ASSOCIATED(r_h))
1553  cpassert(ASSOCIATED(r_s))
1554  cpassert(ASSOCIATED(rho_h))
1555  cpassert(ASSOCIATED(rho_s))
1556  IF (grad_func) THEN
1557  cpassert(ASSOCIATED(dr_h))
1558  cpassert(ASSOCIATED(dr_s))
1559  cpassert(ASSOCIATED(r_h_d))
1560  cpassert(ASSOCIATED(r_s_d))
1561  cpassert(ASSOCIATED(drho_h))
1562  cpassert(ASSOCIATED(drho_s))
1563  END IF
1564 
1565  na = grid_atom%ng_sphere
1566  rad = grid_atom%rad(ir)
1567  urad = grid_atom%oorad2l(ir, 1)
1568  DO ispin = 1, nspins
1569  DO iso = 1, harmonics%max_iso_not0
1570  DO ia = 1, na
1571  rho_h(ia, ir, ispin) = rho_h(ia, ir, ispin) + &
1572  r_h(ispin)%r_coef(ir, iso)*harmonics%slm(ia, iso)
1573  rho_s(ia, ir, ispin) = rho_s(ia, ir, ispin) + &
1574  r_s(ispin)%r_coef(ir, iso)*harmonics%slm(ia, iso)
1575  END DO ! ia
1576  END DO ! iso
1577  END DO ! ispin
1578 
1579  IF (grad_func) THEN
1580  DO ispin = 1, nspins
1581  DO iso = 1, harmonics%max_iso_not0
1582  DO ia = 1, na
1583 
1584  ! components of the gradient of rho1 hard
1585  drho_h(1, ia, ir, ispin) = drho_h(1, ia, ir, ispin) + &
1586  dr_h(ispin)%r_coef(ir, iso)* &
1587  harmonics%a(1, ia)*harmonics%slm(ia, iso) + &
1588  r_h_d(1, ispin)%r_coef(ir, iso)* &
1589  harmonics%slm(ia, iso)
1590 
1591  drho_h(2, ia, ir, ispin) = drho_h(2, ia, ir, ispin) + &
1592  dr_h(ispin)%r_coef(ir, iso)* &
1593  harmonics%a(2, ia)*harmonics%slm(ia, iso) + &
1594  r_h_d(2, ispin)%r_coef(ir, iso)* &
1595  harmonics%slm(ia, iso)
1596 
1597  drho_h(3, ia, ir, ispin) = drho_h(3, ia, ir, ispin) + &
1598  dr_h(ispin)%r_coef(ir, iso)* &
1599  harmonics%a(3, ia)*harmonics%slm(ia, iso) + &
1600  r_h_d(3, ispin)%r_coef(ir, iso)* &
1601  harmonics%slm(ia, iso)
1602 
1603  ! components of the gradient of rho1 soft
1604  drho_s(1, ia, ir, ispin) = drho_s(1, ia, ir, ispin) + &
1605  dr_s(ispin)%r_coef(ir, iso)* &
1606  harmonics%a(1, ia)*harmonics%slm(ia, iso) + &
1607  r_s_d(1, ispin)%r_coef(ir, iso)* &
1608  harmonics%slm(ia, iso)
1609 
1610  drho_s(2, ia, ir, ispin) = drho_s(2, ia, ir, ispin) + &
1611  dr_s(ispin)%r_coef(ir, iso)* &
1612  harmonics%a(2, ia)*harmonics%slm(ia, iso) + &
1613  r_s_d(2, ispin)%r_coef(ir, iso)* &
1614  harmonics%slm(ia, iso)
1615 
1616  drho_s(3, ia, ir, ispin) = drho_s(3, ia, ir, ispin) + &
1617  dr_s(ispin)%r_coef(ir, iso)* &
1618  harmonics%a(3, ia)*harmonics%slm(ia, iso) + &
1619  r_s_d(3, ispin)%r_coef(ir, iso)* &
1620  harmonics%slm(ia, iso)
1621 
1622  drho_h(4, ia, ir, ispin) = sqrt( &
1623  drho_h(1, ia, ir, ispin)*drho_h(1, ia, ir, ispin) + &
1624  drho_h(2, ia, ir, ispin)*drho_h(2, ia, ir, ispin) + &
1625  drho_h(3, ia, ir, ispin)*drho_h(3, ia, ir, ispin))
1626 
1627  drho_s(4, ia, ir, ispin) = sqrt( &
1628  drho_s(1, ia, ir, ispin)*drho_s(1, ia, ir, ispin) + &
1629  drho_s(2, ia, ir, ispin)*drho_s(2, ia, ir, ispin) + &
1630  drho_s(3, ia, ir, ispin)*drho_s(3, ia, ir, ispin))
1631 
1632  END DO ! ia
1633  END DO ! iso
1634  END DO ! ispin
1635  END IF
1636 
1637  END SUBROUTINE calc_rho_angular
1638 
1639 ! **************************************************************************************************
1640 !> \brief Computes tau hard and soft on the atomic grids for meta-GGA calculations
1641 !> \param tau_h the hard part of tau
1642 !> \param tau_s the soft part of tau
1643 !> \param rho_atom ...
1644 !> \param qs_kind ...
1645 !> \param nspins ...
1646 !> \note This is a rewrite to correct a meta-GGA GAPW bug. This is more brute force than the original,
1647 !> which was done along in qs_rho_atom_methods.F, but makes sure that no corner is cut in
1648 !> terms of accuracy (A. Bussy)
1649 ! **************************************************************************************************
1650  SUBROUTINE calc_tau_atom(tau_h, tau_s, rho_atom, qs_kind, nspins)
1651 
1652  REAL(dp), DIMENSION(:, :, :), INTENT(INOUT) :: tau_h, tau_s
1653  TYPE(rho_atom_type), POINTER :: rho_atom
1654  TYPE(qs_kind_type), INTENT(IN) :: qs_kind
1655  INTEGER, INTENT(IN) :: nspins
1656 
1657  CHARACTER(len=*), PARAMETER :: routinen = 'calc_tau_atom'
1658 
1659  INTEGER :: dir, handle, ia, ip, ipgf, ir, iset, &
1660  iso, ispin, jp, jpgf, jset, jso, l, &
1661  maxso, na, nr, nset, starti, startj
1662  INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf, o2nindex
1663  REAL(dp) :: cpc_h, cpc_s
1664  REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: fga, fgr, r1, r2
1665  REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: a1, a2
1666  REAL(dp), DIMENSION(:, :), POINTER :: slm, zet
1667  REAL(dp), DIMENSION(:, :, :), POINTER :: dslm_dxyz
1668  TYPE(grid_atom_type), POINTER :: grid_atom
1669  TYPE(gto_basis_set_type), POINTER :: basis_1c
1670  TYPE(harmonics_atom_type), POINTER :: harmonics
1671 
1672  NULLIFY (harmonics, grid_atom, basis_1c, lmax, lmin, npgf, zet, slm, dslm_dxyz, o2nindex)
1673 
1674  CALL timeset(routinen, handle)
1675 
1676  !We need to put 0.5* grad_g1 dot grad_gw on the grid
1677  !For this we need grid info, basis info, and projector info
1678  CALL get_qs_kind(qs_kind, grid_atom=grid_atom, harmonics=harmonics)
1679  CALL get_qs_kind(qs_kind, basis_set=basis_1c, basis_type="GAPW_1C")
1680 
1681  nr = grid_atom%nr
1682  na = grid_atom%ng_sphere
1683 
1684  slm => harmonics%slm
1685  dslm_dxyz => harmonics%dslm_dxyz
1686 
1687  CALL get_paw_basis_info(basis_1c, o2nindex=o2nindex)
1688 
1689  !zeroing tau, assuming it is already allocated
1690  tau_h = 0.0_dp
1691  tau_s = 0.0_dp
1692 
1693  CALL get_gto_basis_set(gto_basis_set=basis_1c, lmax=lmax, lmin=lmin, npgf=npgf, &
1694  nset=nset, zet=zet, maxso=maxso)
1695 
1696  !Separate the functions into purely r and purely angular parts, precompute them all
1697  ALLOCATE (a1(na, nset*maxso, 3), a2(na, nset*maxso, 3))
1698  ALLOCATE (r1(nr, nset*maxso), r2(nr, nset*maxso))
1699  a1 = 0.0_dp; a2 = 0.0_dp
1700  r1 = 0.0_dp; r2 = 0.0_dp
1701 
1702  DO iset = 1, nset
1703  DO ipgf = 1, npgf(iset)
1704  starti = (iset - 1)*maxso + (ipgf - 1)*nsoset(lmax(iset))
1705  DO iso = nsoset(lmin(iset) - 1) + 1, nsoset(lmax(iset))
1706  l = indso(1, iso)
1707 
1708  ! The x derivative of the spherical orbital, divided in angular and radial parts
1709  ! Two of each are needed because d/dx(r^l Y_lm) * exp(-al*r^2) + r^l Y_lm * ! d/dx(exp-al*r^2)
1710 
1711  ! the purely radial part of d/dx(r^l Y_lm) * exp(-al*r^2) (same for y, z)
1712  r1(1:nr, starti + iso) = grid_atom%rad(1:nr)**(l - 1)*exp(-zet(ipgf, iset)*grid_atom%rad2(1:nr))
1713 
1714  ! the purely radial part of r^l Y_lm * d/dx(exp-al*r^2) (same for y, z)
1715  r2(1:nr, starti + iso) = -2.0_dp*zet(ipgf, iset)*grid_atom%rad(1:nr)**(l + 1) &
1716  *exp(-zet(ipgf, iset)*grid_atom%rad2(1:nr))
1717 
1718  DO dir = 1, 3
1719  ! the purely angular part of d/dx(r^l Y_lm) * exp(-al*r^2)
1720  a1(1:na, starti + iso, dir) = dslm_dxyz(dir, 1:na, iso)
1721 
1722  ! the purely angular part of r^l Y_lm * d/dx(exp-al*r^2)
1723  a2(1:na, starti + iso, dir) = harmonics%a(dir, 1:na)*slm(1:na, iso)
1724  END DO
1725 
1726  END DO !iso
1727  END DO !ipgf
1728  END DO !iset
1729 
1730  !Compute the matrix products
1731  ALLOCATE (fga(na, 1))
1732  ALLOCATE (fgr(nr, 1))
1733  fga = 0.0_dp; fgr = 0.0_dp
1734 
1735  DO iset = 1, nset
1736  DO jset = 1, nset
1737  DO ipgf = 1, npgf(iset)
1738  starti = (iset - 1)*maxso + (ipgf - 1)*nsoset(lmax(iset))
1739  DO jpgf = 1, npgf(jset)
1740  startj = (jset - 1)*maxso + (jpgf - 1)*nsoset(lmax(jset))
1741 
1742  DO iso = nsoset(lmin(iset) - 1) + 1, nsoset(lmax(iset))
1743  DO jso = nsoset(lmin(jset) - 1) + 1, nsoset(lmax(jset))
1744 
1745  ip = o2nindex(starti + iso)
1746  jp = o2nindex(startj + jso)
1747 
1748  !Two component per function => 4 terms in total
1749 
1750  ! take r1*a1(dir) * r1*a1(dir)
1751  fgr(1:nr, 1) = r1(1:nr, starti + iso)*r1(1:nr, startj + jso)
1752  DO dir = 1, 3
1753  fga(1:na, 1) = a1(1:na, starti + iso, dir)*a1(1:na, startj + jso, dir)
1754 
1755  DO ispin = 1, nspins
1756  !get the projectors
1757  cpc_h = rho_atom%cpc_h(ispin)%r_coef(ip, jp)
1758  cpc_s = rho_atom%cpc_s(ispin)%r_coef(ip, jp)
1759 
1760  !compute contribution to tau
1761  DO ir = 1, nr
1762  DO ia = 1, na
1763  tau_h(ia, ir, ispin) = tau_h(ia, ir, ispin) + 0.5_dp*cpc_h* &
1764  fgr(ir, 1)*fga(ia, 1)
1765 
1766  tau_s(ia, ir, ispin) = tau_s(ia, ir, ispin) + 0.5_dp*cpc_s* &
1767  fgr(ir, 1)*fga(ia, 1)
1768  END DO
1769  END DO
1770 
1771  END DO !ispin
1772  END DO !dir
1773 
1774  ! add r1*a1(dir) * r2*a2(dir)
1775  fgr(1:nr, 1) = r1(1:nr, starti + iso)*r2(1:nr, startj + jso)
1776  DO dir = 1, 3
1777  fga(1:na, 1) = a1(1:na, starti + iso, dir)*a2(1:na, startj + jso, dir)
1778 
1779  DO ispin = 1, nspins
1780  !get the projectors
1781  cpc_h = rho_atom%cpc_h(ispin)%r_coef(ip, jp)
1782  cpc_s = rho_atom%cpc_s(ispin)%r_coef(ip, jp)
1783 
1784  !compute contribution to tau
1785  DO ir = 1, nr
1786  DO ia = 1, na
1787  tau_h(ia, ir, ispin) = tau_h(ia, ir, ispin) + 0.5_dp*cpc_h* &
1788  fgr(ir, 1)*fga(ia, 1)
1789 
1790  tau_s(ia, ir, ispin) = tau_s(ia, ir, ispin) + 0.5_dp*cpc_s* &
1791  fgr(ir, 1)*fga(ia, 1)
1792  END DO
1793  END DO
1794 
1795  END DO !ispin
1796  END DO !dir
1797 
1798  ! add r2*a2(dir) * V * r1*a1(dir)
1799  fgr(1:nr, 1) = r2(1:nr, starti + iso)*r1(1:nr, startj + jso)
1800  DO dir = 1, 3
1801  fga(1:na, 1) = a2(1:na, starti + iso, dir)*a1(1:na, startj + jso, dir)
1802 
1803  DO ispin = 1, nspins
1804  !get the projectors
1805  cpc_h = rho_atom%cpc_h(ispin)%r_coef(ip, jp)
1806  cpc_s = rho_atom%cpc_s(ispin)%r_coef(ip, jp)
1807 
1808  !compute contribution to tau
1809  DO ir = 1, nr
1810  DO ia = 1, na
1811  tau_h(ia, ir, ispin) = tau_h(ia, ir, ispin) + 0.5_dp*cpc_h* &
1812  fgr(ir, 1)*fga(ia, 1)
1813 
1814  tau_s(ia, ir, ispin) = tau_s(ia, ir, ispin) + 0.5_dp*cpc_s* &
1815  fgr(ir, 1)*fga(ia, 1)
1816  END DO
1817  END DO
1818 
1819  END DO !ispin
1820  END DO !dir
1821 
1822  ! add the last term: r2*a2(dir) * r2*a2(dir)
1823  fgr(1:nr, 1) = r2(1:nr, starti + iso)*r2(1:nr, startj + jso)
1824  DO dir = 1, 3
1825  fga(1:na, 1) = a2(1:na, starti + iso, dir)*a2(1:na, startj + jso, dir)
1826 
1827  DO ispin = 1, nspins
1828  !get the projectors
1829  cpc_h = rho_atom%cpc_h(ispin)%r_coef(ip, jp)
1830  cpc_s = rho_atom%cpc_s(ispin)%r_coef(ip, jp)
1831 
1832  !compute contribution to tau
1833  DO ir = 1, nr
1834  DO ia = 1, na
1835  tau_h(ia, ir, ispin) = tau_h(ia, ir, ispin) + 0.5_dp*cpc_h* &
1836  fgr(ir, 1)*fga(ia, 1)
1837 
1838  tau_s(ia, ir, ispin) = tau_s(ia, ir, ispin) + 0.5_dp*cpc_s* &
1839  fgr(ir, 1)*fga(ia, 1)
1840  END DO
1841  END DO
1842 
1843  END DO !ispin
1844  END DO !dir
1845 
1846  END DO !jso
1847  END DO !iso
1848 
1849  END DO !jpgf
1850  END DO !ipgf
1851  END DO !jset
1852  END DO !iset
1853 
1854  DEALLOCATE (o2nindex)
1855 
1856  CALL timestop(handle)
1857 
1858  END SUBROUTINE calc_tau_atom
1859 
1860 ! **************************************************************************************************
1861 !> \brief ...
1862 !> \param grid_atom ...
1863 !> \param nspins ...
1864 !> \param grad_func ...
1865 !> \param ir ...
1866 !> \param rho_nlcc ...
1867 !> \param rho_h ...
1868 !> \param rho_s ...
1869 !> \param drho_nlcc ...
1870 !> \param drho_h ...
1871 !> \param drho_s ...
1872 ! **************************************************************************************************
1873  SUBROUTINE calc_rho_nlcc(grid_atom, nspins, grad_func, &
1874  ir, rho_nlcc, rho_h, rho_s, drho_nlcc, drho_h, drho_s)
1875 
1876  TYPE(grid_atom_type), POINTER :: grid_atom
1877  INTEGER, INTENT(IN) :: nspins
1878  LOGICAL, INTENT(IN) :: grad_func
1879  INTEGER, INTENT(IN) :: ir
1880  REAL(kind=dp), DIMENSION(:) :: rho_nlcc
1881  REAL(kind=dp), DIMENSION(:, :, :), POINTER :: rho_h, rho_s
1882  REAL(kind=dp), DIMENSION(:) :: drho_nlcc
1883  REAL(kind=dp), DIMENSION(:, :, :, :), POINTER :: drho_h, drho_s
1884 
1885  INTEGER :: ia, ispin, na
1886  REAL(kind=dp) :: drho, dx, dy, dz, rad, rho, urad, xsp
1887 
1888  cpassert(ASSOCIATED(rho_h))
1889  cpassert(ASSOCIATED(rho_s))
1890  IF (grad_func) THEN
1891  cpassert(ASSOCIATED(drho_h))
1892  cpassert(ASSOCIATED(drho_s))
1893  END IF
1894 
1895  na = grid_atom%ng_sphere
1896  rad = grid_atom%rad(ir)
1897  urad = grid_atom%oorad2l(ir, 1)
1898 
1899  xsp = real(nspins, kind=dp)
1900  rho = rho_nlcc(ir)/xsp
1901  DO ispin = 1, nspins
1902  rho_h(1:na, ir, ispin) = rho_h(1:na, ir, ispin) + rho
1903  rho_s(1:na, ir, ispin) = rho_s(1:na, ir, ispin) + rho
1904  END DO ! ispin
1905 
1906  IF (grad_func) THEN
1907  drho = drho_nlcc(ir)/xsp
1908  DO ispin = 1, nspins
1909  DO ia = 1, na
1910  IF (grid_atom%azi(ia) == 0.0_dp) THEN
1911  dx = 0.0_dp
1912  dy = 0.0_dp
1913  ELSE
1914  dx = grid_atom%sin_pol(ia)*grid_atom%sin_azi(ia)
1915  dy = grid_atom%sin_pol(ia)*grid_atom%cos_azi(ia)
1916  END IF
1917  dz = grid_atom%cos_pol(ia)
1918  ! components of the gradient of rho1 hard
1919  drho_h(1, ia, ir, ispin) = drho_h(1, ia, ir, ispin) + drho*dx
1920  drho_h(2, ia, ir, ispin) = drho_h(2, ia, ir, ispin) + drho*dy
1921  drho_h(3, ia, ir, ispin) = drho_h(3, ia, ir, ispin) + drho*dz
1922  ! components of the gradient of rho1 soft
1923  drho_s(1, ia, ir, ispin) = drho_s(1, ia, ir, ispin) + drho*dx
1924  drho_s(2, ia, ir, ispin) = drho_s(2, ia, ir, ispin) + drho*dy
1925  drho_s(3, ia, ir, ispin) = drho_s(3, ia, ir, ispin) + drho*dz
1926  ! norm of gradient
1927  drho_h(4, ia, ir, ispin) = sqrt( &
1928  drho_h(1, ia, ir, ispin)*drho_h(1, ia, ir, ispin) + &
1929  drho_h(2, ia, ir, ispin)*drho_h(2, ia, ir, ispin) + &
1930  drho_h(3, ia, ir, ispin)*drho_h(3, ia, ir, ispin))
1931 
1932  drho_s(4, ia, ir, ispin) = sqrt( &
1933  drho_s(1, ia, ir, ispin)*drho_s(1, ia, ir, ispin) + &
1934  drho_s(2, ia, ir, ispin)*drho_s(2, ia, ir, ispin) + &
1935  drho_s(3, ia, ir, ispin)*drho_s(3, ia, ir, ispin))
1936  END DO ! ia
1937  END DO ! ispin
1938  END IF
1939 
1940  END SUBROUTINE calc_rho_nlcc
1941 
1942 ! **************************************************************************************************
1943 !> \brief ...
1944 !> \param vxc_h ...
1945 !> \param vxc_s ...
1946 !> \param int_hh ...
1947 !> \param int_ss ...
1948 !> \param grid_atom ...
1949 !> \param basis_1c ...
1950 !> \param harmonics ...
1951 !> \param nspins ...
1952 ! **************************************************************************************************
1953  SUBROUTINE gavxcgb_nogc(vxc_h, vxc_s, int_hh, int_ss, grid_atom, basis_1c, harmonics, nspins)
1954 
1955  REAL(dp), DIMENSION(:, :, :), POINTER :: vxc_h, vxc_s
1956  TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: int_hh, int_ss
1957  TYPE(grid_atom_type), POINTER :: grid_atom
1958  TYPE(gto_basis_set_type), POINTER :: basis_1c
1959  TYPE(harmonics_atom_type), POINTER :: harmonics
1960  INTEGER, INTENT(IN) :: nspins
1961 
1962  CHARACTER(len=*), PARAMETER :: routinen = 'gaVxcgb_noGC'
1963 
1964  INTEGER :: handle, ia, ic, icg, ipgf1, ipgf2, ir, iset1, iset2, iso, iso1, iso2, ispin, l, &
1965  ld, lmax12, lmax_expansion, lmin12, m1, m2, max_iso_not0, max_iso_not0_local, max_s_harm, &
1966  maxl, maxso, n1, n2, na, ngau1, ngau2, nngau1, nr, nset, size1
1967  INTEGER, ALLOCATABLE, DIMENSION(:) :: cg_n_list
1968  INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cg_list
1969  INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf
1970  REAL(dp), ALLOCATABLE, DIMENSION(:) :: g1, g2
1971  REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: gg, gvg_h, gvg_s, matso_h, matso_s, vx
1972  REAL(dp), DIMENSION(:, :), POINTER :: zet
1973  REAL(dp), DIMENSION(:, :, :), POINTER :: my_cg
1974 
1975  CALL timeset(routinen, handle)
1976 
1977  NULLIFY (lmin, lmax, npgf, zet, my_cg)
1978 
1979  CALL get_gto_basis_set(gto_basis_set=basis_1c, lmax=lmax, lmin=lmin, &
1980  maxso=maxso, maxl=maxl, npgf=npgf, &
1981  nset=nset, zet=zet)
1982 
1983  nr = grid_atom%nr
1984  na = grid_atom%ng_sphere
1985  my_cg => harmonics%my_CG
1986  max_iso_not0 = harmonics%max_iso_not0
1987  lmax_expansion = indso(1, max_iso_not0)
1988  max_s_harm = harmonics%max_s_harm
1989 
1990  ALLOCATE (g1(nr), g2(nr), gg(nr, 0:2*maxl))
1991  ALLOCATE (gvg_h(na, 0:2*maxl), gvg_s(na, 0:2*maxl))
1992  ALLOCATE (matso_h(nsoset(maxl), nsoset(maxl)), &
1993  matso_s(nsoset(maxl), nsoset(maxl)))
1994  ALLOCATE (vx(na, nr))
1995  ALLOCATE (cg_list(2, nsoset(maxl)**2, max_s_harm), cg_n_list(max_s_harm))
1996 
1997  g1 = 0.0_dp
1998  g2 = 0.0_dp
1999  m1 = 0
2000  DO iset1 = 1, nset
2001  n1 = nsoset(lmax(iset1))
2002  m2 = 0
2003  DO iset2 = 1, nset
2004  CALL get_none0_cg_list(my_cg, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
2005  max_s_harm, lmax_expansion, cg_list, cg_n_list, max_iso_not0_local)
2006  cpassert(max_iso_not0_local .LE. max_iso_not0)
2007 
2008  n2 = nsoset(lmax(iset2))
2009  DO ipgf1 = 1, npgf(iset1)
2010  ngau1 = n1*(ipgf1 - 1) + m1
2011  size1 = nsoset(lmax(iset1)) - nsoset(lmin(iset1) - 1)
2012  nngau1 = nsoset(lmin(iset1) - 1) + ngau1
2013 
2014  g1(1:nr) = exp(-zet(ipgf1, iset1)*grid_atom%rad2(1:nr))
2015  DO ipgf2 = 1, npgf(iset2)
2016  ngau2 = n2*(ipgf2 - 1) + m2
2017 
2018  g2(1:nr) = exp(-zet(ipgf2, iset2)*grid_atom%rad2(1:nr))
2019  lmin12 = lmin(iset1) + lmin(iset2)
2020  lmax12 = lmax(iset1) + lmax(iset2)
2021 
2022  ! reduce expansion local densities
2023  IF (lmin12 .LE. lmax_expansion) THEN
2024 
2025  gg = 0.0_dp
2026  IF (lmin12 == 0) THEN
2027  gg(1:nr, lmin12) = g1(1:nr)*g2(1:nr)
2028  ELSE
2029  gg(1:nr, lmin12) = grid_atom%rad2l(1:nr, lmin12)*g1(1:nr)*g2(1:nr)
2030  END IF
2031 
2032  ! limit the expansion of the local densities to a max L
2033  IF (lmax12 .GT. lmax_expansion) lmax12 = lmax_expansion
2034 
2035  DO l = lmin12 + 1, lmax12
2036  gg(1:nr, l) = grid_atom%rad(1:nr)*gg(:, l - 1)
2037  END DO
2038 
2039  DO ispin = 1, nspins
2040  ld = lmax12 + 1
2041  DO ir = 1, nr
2042  vx(1:na, ir) = vxc_h(1:na, ir, ispin)
2043  END DO
2044  CALL dgemm('N', 'N', na, ld, nr, 1.0_dp, vx(1:na, 1:nr), na, &
2045  gg(1:nr, 0:lmax12), nr, 0.0_dp, gvg_h(1:na, 0:lmax12), na)
2046  DO ir = 1, nr
2047  vx(1:na, ir) = vxc_s(1:na, ir, ispin)
2048  END DO
2049  CALL dgemm('N', 'N', na, ld, nr, 1.0_dp, vx(1:na, 1:nr), na, &
2050  gg(1:nr, 0:lmax12), nr, 0.0_dp, gvg_s(1:na, 0:lmax12), na)
2051 
2052  matso_h = 0.0_dp
2053  matso_s = 0.0_dp
2054  DO iso = 1, max_iso_not0_local
2055  DO icg = 1, cg_n_list(iso)
2056  iso1 = cg_list(1, icg, iso)
2057  iso2 = cg_list(2, icg, iso)
2058  l = indso(1, iso1) + indso(1, iso2)
2059 
2060  cpassert(l <= lmax_expansion)
2061  DO ia = 1, na
2062  matso_h(iso1, iso2) = matso_h(iso1, iso2) + &
2063  gvg_h(ia, l)* &
2064  my_cg(iso1, iso2, iso)* &
2065  harmonics%slm(ia, iso)
2066  matso_s(iso1, iso2) = matso_s(iso1, iso2) + &
2067  gvg_s(ia, l)* &
2068  my_cg(iso1, iso2, iso)* &
2069  harmonics%slm(ia, iso)
2070  END DO
2071  END DO
2072  END DO
2073 
2074  ! Write in the global matrix
2075  DO ic = nsoset(lmin(iset2) - 1) + 1, nsoset(lmax(iset2))
2076  iso1 = nsoset(lmin(iset1) - 1) + 1
2077  iso2 = ngau2 + ic
2078  CALL daxpy(size1, 1.0_dp, matso_h(iso1, ic), 1, &
2079  int_hh(ispin)%r_coef(nngau1 + 1, iso2), 1)
2080  CALL daxpy(size1, 1.0_dp, matso_s(iso1, ic), 1, &
2081  int_ss(ispin)%r_coef(nngau1 + 1, iso2), 1)
2082  END DO
2083 
2084  END DO ! ispin
2085 
2086  END IF ! lmax_expansion
2087 
2088  END DO ! ipfg2
2089  END DO ! ipfg1
2090  m2 = m2 + maxso
2091  END DO ! iset2
2092  m1 = m1 + maxso
2093  END DO ! iset1
2094 
2095  DEALLOCATE (g1, g2, gg, matso_h, matso_s, gvg_s, gvg_h, vx)
2096 
2097  DEALLOCATE (cg_list, cg_n_list)
2098 
2099  CALL timestop(handle)
2100 
2101  END SUBROUTINE gavxcgb_nogc
2102 
2103 ! **************************************************************************************************
2104 !> \brief ...
2105 !> \param vxc_h ...
2106 !> \param vxc_s ...
2107 !> \param vxg_h ...
2108 !> \param vxg_s ...
2109 !> \param int_hh ...
2110 !> \param int_ss ...
2111 !> \param grid_atom ...
2112 !> \param basis_1c ...
2113 !> \param harmonics ...
2114 !> \param nspins ...
2115 ! **************************************************************************************************
2116  SUBROUTINE gavxcgb_gc(vxc_h, vxc_s, vxg_h, vxg_s, int_hh, int_ss, &
2117  grid_atom, basis_1c, harmonics, nspins)
2118 
2119  REAL(dp), DIMENSION(:, :, :), POINTER :: vxc_h, vxc_s
2120  REAL(dp), DIMENSION(:, :, :, :), POINTER :: vxg_h, vxg_s
2121  TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: int_hh, int_ss
2122  TYPE(grid_atom_type), POINTER :: grid_atom
2123  TYPE(gto_basis_set_type), POINTER :: basis_1c
2124  TYPE(harmonics_atom_type), POINTER :: harmonics
2125  INTEGER, INTENT(IN) :: nspins
2126 
2127  CHARACTER(len=*), PARAMETER :: routinen = 'gaVxcgb_GC'
2128 
2129  INTEGER :: dmax_iso_not0_local, handle, ia, ic, icg, ipgf1, ipgf2, ir, iset1, iset2, iso, &
2130  iso1, iso2, ispin, l, lmax12, lmax_expansion, lmin12, m1, m2, max_iso_not0, &
2131  max_iso_not0_local, max_s_harm, maxl, maxso, n1, n2, na, ngau1, ngau2, nngau1, nr, nset, &
2132  size1
2133  INTEGER, ALLOCATABLE, DIMENSION(:) :: cg_n_list, dcg_n_list
2134  INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cg_list, dcg_list
2135  INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf
2136  REAL(dp) :: urad
2137  REAL(dp), ALLOCATABLE, DIMENSION(:) :: g1, g2
2138  REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: dgg, gg, gvxcg_h, gvxcg_s, matso_h, &
2139  matso_s
2140  REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: gvxgg_h, gvxgg_s
2141  REAL(dp), DIMENSION(:, :), POINTER :: zet
2142  REAL(dp), DIMENSION(:, :, :), POINTER :: my_cg
2143  REAL(dp), DIMENSION(:, :, :, :), POINTER :: my_cg_dxyz
2144 
2145  CALL timeset(routinen, handle)
2146 
2147  NULLIFY (lmin, lmax, npgf, zet, my_cg, my_cg_dxyz)
2148 
2149  CALL get_gto_basis_set(gto_basis_set=basis_1c, lmax=lmax, lmin=lmin, &
2150  maxso=maxso, maxl=maxl, npgf=npgf, &
2151  nset=nset, zet=zet)
2152 
2153  nr = grid_atom%nr
2154  na = grid_atom%ng_sphere
2155  my_cg => harmonics%my_CG
2156  my_cg_dxyz => harmonics%my_CG_dxyz
2157  max_iso_not0 = harmonics%max_iso_not0
2158  lmax_expansion = indso(1, max_iso_not0)
2159  max_s_harm = harmonics%max_s_harm
2160 
2161  ALLOCATE (g1(nr), g2(nr), gg(nr, 0:2*maxl), dgg(nr, 0:2*maxl))
2162  ALLOCATE (gvxcg_h(na, 0:2*maxl), gvxcg_s(na, 0:2*maxl))
2163  ALLOCATE (gvxgg_h(3, na, 0:2*maxl), gvxgg_s(3, na, 0:2*maxl))
2164  ALLOCATE (cg_list(2, nsoset(maxl)**2, max_s_harm), cg_n_list(max_s_harm), &
2165  dcg_list(2, nsoset(maxl)**2, max_s_harm), dcg_n_list(max_s_harm))
2166 
2167  ALLOCATE (matso_h(nsoset(maxl), nsoset(maxl)), &
2168  matso_s(nsoset(maxl), nsoset(maxl)))
2169 
2170  DO ispin = 1, nspins
2171 
2172  g1 = 0.0_dp
2173  g2 = 0.0_dp
2174  m1 = 0
2175  DO iset1 = 1, nset
2176  n1 = nsoset(lmax(iset1))
2177  m2 = 0
2178  DO iset2 = 1, nset
2179  CALL get_none0_cg_list(my_cg, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
2180  max_s_harm, lmax_expansion, cg_list, cg_n_list, max_iso_not0_local)
2181  cpassert(max_iso_not0_local .LE. max_iso_not0)
2182  CALL get_none0_cg_list(my_cg_dxyz, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
2183  max_s_harm, lmax_expansion, dcg_list, dcg_n_list, dmax_iso_not0_local)
2184 
2185  n2 = nsoset(lmax(iset2))
2186  DO ipgf1 = 1, npgf(iset1)
2187  ngau1 = n1*(ipgf1 - 1) + m1
2188  size1 = nsoset(lmax(iset1)) - nsoset(lmin(iset1) - 1)
2189  nngau1 = nsoset(lmin(iset1) - 1) + ngau1
2190 
2191  g1(1:nr) = exp(-zet(ipgf1, iset1)*grid_atom%rad2(1:nr))
2192  DO ipgf2 = 1, npgf(iset2)
2193  ngau2 = n2*(ipgf2 - 1) + m2
2194 
2195  g2(1:nr) = exp(-zet(ipgf2, iset2)*grid_atom%rad2(1:nr))
2196  lmin12 = lmin(iset1) + lmin(iset2)
2197  lmax12 = lmax(iset1) + lmax(iset2)
2198 
2199  !test reduce expansion local densities
2200  IF (lmin12 .LE. lmax_expansion) THEN
2201 
2202  gg = 0.0_dp
2203  dgg = 0.0_dp
2204 
2205  IF (lmin12 == 0) THEN
2206  gg(1:nr, lmin12) = g1(1:nr)*g2(1:nr)
2207  ELSE
2208  gg(1:nr, lmin12) = grid_atom%rad2l(1:nr, lmin12)*g1(1:nr)*g2(1:nr)
2209  END IF
2210 
2211  !test reduce expansion local densities
2212  IF (lmax12 .GT. lmax_expansion) lmax12 = lmax_expansion
2213 
2214  DO l = lmin12 + 1, lmax12
2215  gg(1:nr, l) = grid_atom%rad(1:nr)*gg(:, l - 1)
2216  dgg(1:nr, l - 1) = dgg(1:nr, l - 1) - 2.0_dp*(zet(ipgf1, iset1) + &
2217  zet(ipgf2, iset2))*gg(1:nr, l)
2218  END DO
2219  dgg(1:nr, lmax12) = dgg(1:nr, lmax12) - 2.0_dp*(zet(ipgf1, iset1) + &
2220  zet(ipgf2, iset2))*grid_atom%rad(1:nr)* &
2221  gg(1:nr, lmax12)
2222 
2223  gvxcg_h = 0.0_dp
2224  gvxcg_s = 0.0_dp
2225  gvxgg_h = 0.0_dp
2226  gvxgg_s = 0.0_dp
2227 
2228  ! Cross Term
2229  DO l = lmin12, lmax12
2230  DO ia = 1, na
2231  DO ir = 1, nr
2232  gvxcg_h(ia, l) = gvxcg_h(ia, l) + &
2233  gg(ir, l)*vxc_h(ia, ir, ispin) + &
2234  dgg(ir, l)* &
2235  (vxg_h(1, ia, ir, ispin)*harmonics%a(1, ia) + &
2236  vxg_h(2, ia, ir, ispin)*harmonics%a(2, ia) + &
2237  vxg_h(3, ia, ir, ispin)*harmonics%a(3, ia))
2238 
2239  gvxcg_s(ia, l) = gvxcg_s(ia, l) + &
2240  gg(ir, l)*vxc_s(ia, ir, ispin) + &
2241  dgg(ir, l)* &
2242  (vxg_s(1, ia, ir, ispin)*harmonics%a(1, ia) + &
2243  vxg_s(2, ia, ir, ispin)*harmonics%a(2, ia) + &
2244  vxg_s(3, ia, ir, ispin)*harmonics%a(3, ia))
2245 
2246  urad = grid_atom%oorad2l(ir, 1)
2247 
2248  gvxgg_h(1, ia, l) = gvxgg_h(1, ia, l) + &
2249  vxg_h(1, ia, ir, ispin)* &
2250  gg(ir, l)*urad
2251 
2252  gvxgg_h(2, ia, l) = gvxgg_h(2, ia, l) + &
2253  vxg_h(2, ia, ir, ispin)* &
2254  gg(ir, l)*urad
2255 
2256  gvxgg_h(3, ia, l) = gvxgg_h(3, ia, l) + &
2257  vxg_h(3, ia, ir, ispin)* &
2258  gg(ir, l)*urad
2259 
2260  gvxgg_s(1, ia, l) = gvxgg_s(1, ia, l) + &
2261  vxg_s(1, ia, ir, ispin)* &
2262  gg(ir, l)*urad
2263 
2264  gvxgg_s(2, ia, l) = gvxgg_s(2, ia, l) + &
2265  vxg_s(2, ia, ir, ispin)* &
2266  gg(ir, l)*urad
2267 
2268  gvxgg_s(3, ia, l) = gvxgg_s(3, ia, l) + &
2269  vxg_s(3, ia, ir, ispin)* &
2270  gg(ir, l)*urad
2271 
2272  END DO ! ir
2273  END DO ! ia
2274  END DO ! l
2275 
2276  matso_h = 0.0_dp
2277  matso_s = 0.0_dp
2278  DO iso = 1, max_iso_not0_local
2279  DO icg = 1, cg_n_list(iso)
2280  iso1 = cg_list(1, icg, iso)
2281  iso2 = cg_list(2, icg, iso)
2282 
2283  l = indso(1, iso1) + indso(1, iso2)
2284 
2285  !test reduce expansion local densities
2286  cpassert(l <= lmax_expansion)
2287  DO ia = 1, na
2288  matso_h(iso1, iso2) = matso_h(iso1, iso2) + &
2289  gvxcg_h(ia, l)* &
2290  harmonics%slm(ia, iso)* &
2291  my_cg(iso1, iso2, iso)
2292  matso_s(iso1, iso2) = matso_s(iso1, iso2) + &
2293  gvxcg_s(ia, l)* &
2294  harmonics%slm(ia, iso)* &
2295  my_cg(iso1, iso2, iso)
2296  END DO ! ia
2297 
2298  !test reduce expansion local densities
2299 
2300  END DO
2301 
2302  END DO ! iso
2303 
2304  DO iso = 1, dmax_iso_not0_local
2305  DO icg = 1, dcg_n_list(iso)
2306  iso1 = dcg_list(1, icg, iso)
2307  iso2 = dcg_list(2, icg, iso)
2308 
2309  l = indso(1, iso1) + indso(1, iso2)
2310  !test reduce expansion local densities
2311  cpassert(l <= lmax_expansion)
2312  DO ia = 1, na
2313  matso_h(iso1, iso2) = matso_h(iso1, iso2) + &
2314  (gvxgg_h(1, ia, l)*my_cg_dxyz(1, iso1, iso2, iso) + &
2315  gvxgg_h(2, ia, l)*my_cg_dxyz(2, iso1, iso2, iso) + &
2316  gvxgg_h(3, ia, l)*my_cg_dxyz(3, iso1, iso2, iso))* &
2317  harmonics%slm(ia, iso)
2318 
2319  matso_s(iso1, iso2) = matso_s(iso1, iso2) + &
2320  (gvxgg_s(1, ia, l)*my_cg_dxyz(1, iso1, iso2, iso) + &
2321  gvxgg_s(2, ia, l)*my_cg_dxyz(2, iso1, iso2, iso) + &
2322  gvxgg_s(3, ia, l)*my_cg_dxyz(3, iso1, iso2, iso))* &
2323  harmonics%slm(ia, iso)
2324 
2325  END DO ! ia
2326 
2327  !test reduce expansion local densities
2328 
2329  END DO ! icg
2330  END DO ! iso
2331  !test reduce expansion local densities
2332  END IF ! lmax_expansion
2333 
2334  ! Write in the global matrix
2335  DO ic = nsoset(lmin(iset2) - 1) + 1, nsoset(lmax(iset2))
2336  iso1 = nsoset(lmin(iset1) - 1) + 1
2337  iso2 = ngau2 + ic
2338  CALL daxpy(size1, 1.0_dp, matso_h(iso1, ic), 1, &
2339  int_hh(ispin)%r_coef(nngau1 + 1, iso2), 1)
2340  CALL daxpy(size1, 1.0_dp, matso_s(iso1, ic), 1, &
2341  int_ss(ispin)%r_coef(nngau1 + 1, iso2), 1)
2342  END DO
2343 
2344  END DO ! ipfg2
2345  END DO ! ipfg1
2346  m2 = m2 + maxso
2347  END DO ! iset2
2348  m1 = m1 + maxso
2349  END DO ! iset1
2350  END DO ! ispin
2351 
2352  DEALLOCATE (g1, g2, gg, dgg, matso_h, matso_s, gvxcg_h, gvxcg_s, gvxgg_h, gvxgg_s)
2353  DEALLOCATE (cg_list, cg_n_list, dcg_list, dcg_n_list)
2354 
2355  CALL timestop(handle)
2356 
2357  END SUBROUTINE gavxcgb_gc
2358 
2359 ! **************************************************************************************************
2360 !> \brief Integrates 0.5 * grad_ga .dot. (V_tau * grad_gb) on the atomic grid for meta-GGA
2361 !> \param vtau_h the har tau potential
2362 !> \param vtau_s the soft tau potential
2363 !> \param int_hh ...
2364 !> \param int_ss ...
2365 !> \param grid_atom ...
2366 !> \param basis_1c ...
2367 !> \param harmonics ...
2368 !> \param nspins ...
2369 !> \note This is a rewrite to correct meta-GGA GAPW bug. This is more brute force than the original
2370 !> but makes sure that no corner is cut in terms of accuracy (A. Bussy)
2371 ! **************************************************************************************************
2372  SUBROUTINE dgavtaudgb(vtau_h, vtau_s, int_hh, int_ss, &
2373  grid_atom, basis_1c, harmonics, nspins)
2374 
2375  REAL(dp), DIMENSION(:, :, :), POINTER :: vtau_h, vtau_s
2376  TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: int_hh, int_ss
2377  TYPE(grid_atom_type), POINTER :: grid_atom
2378  TYPE(gto_basis_set_type), POINTER :: basis_1c
2379  TYPE(harmonics_atom_type), POINTER :: harmonics
2380  INTEGER, INTENT(IN) :: nspins
2381 
2382  CHARACTER(len=*), PARAMETER :: routinen = 'dgaVtaudgb'
2383 
2384  INTEGER :: dir, handle, ipgf, iset, iso, ispin, &
2385  jpgf, jset, jso, l, maxso, na, nr, &
2386  nset, starti, startj
2387  INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf
2388  REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: fga, fgr, r1, r2, work
2389  REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: a1, a2, intso_h, intso_s
2390  REAL(dp), DIMENSION(:, :), POINTER :: slm, zet
2391  REAL(dp), DIMENSION(:, :, :), POINTER :: dslm_dxyz
2392 
2393  CALL timeset(routinen, handle)
2394 
2395  NULLIFY (zet, slm, dslm_dxyz, lmax, lmin, npgf)
2396 
2397  CALL get_gto_basis_set(gto_basis_set=basis_1c, lmax=lmax, lmin=lmin, &
2398  maxso=maxso, npgf=npgf, nset=nset, zet=zet)
2399 
2400  na = grid_atom%ng_sphere
2401  nr = grid_atom%nr
2402 
2403  slm => harmonics%slm
2404  dslm_dxyz => harmonics%dslm_dxyz
2405 
2406  ! Separate the functions into purely r and purely angular parts and precompute them all
2407  ! Not memory intensive since 1D arrays
2408  ALLOCATE (a1(na, nset*maxso, 3), a2(na, nset*maxso, 3))
2409  ALLOCATE (r1(nr, nset*maxso), r2(nr, nset*maxso))
2410  a1 = 0.0_dp; a2 = 0.0_dp
2411  r1 = 0.0_dp; r2 = 0.0_dp
2412 
2413  DO iset = 1, nset
2414  DO ipgf = 1, npgf(iset)
2415  starti = (iset - 1)*maxso + (ipgf - 1)*nsoset(lmax(iset))
2416  DO iso = nsoset(lmin(iset) - 1) + 1, nsoset(lmax(iset))
2417  l = indso(1, iso)
2418 
2419  ! The x derivative of the spherical orbital, divided in angular and radial parts
2420  ! Two of each are needed because d/dx(r^l Y_lm) * exp(-al*r^2) + r^l Y_lm * d/dx(exp-al*r^2)
2421 
2422  ! the purely radial part of d/dx(r^l Y_lm) * exp(-al*r^2) (same for y,z)
2423  r1(1:nr, starti + iso) = grid_atom%rad(1:nr)**(l - 1)*exp(-zet(ipgf, iset)*grid_atom%rad2(1:nr))
2424 
2425  ! the purely radial part of r^l Y_lm * d/dx(exp-al*r^2) (same for y,z)
2426  r2(1:nr, starti + iso) = -2.0_dp*zet(ipgf, iset)*grid_atom%rad(1:nr)**(l + 1) &
2427  *exp(-zet(ipgf, iset)*grid_atom%rad2(1:nr))
2428 
2429  DO dir = 1, 3
2430  ! the purely angular part of d/dx(r^l Y_lm) * exp(-al*r^2)
2431  a1(1:na, starti + iso, dir) = dslm_dxyz(dir, 1:na, iso)
2432 
2433  ! the purely angular part of r^l Y_lm * d/dx(exp-al*r^2)
2434  a2(1:na, starti + iso, dir) = harmonics%a(dir, 1:na)*slm(1:na, iso)
2435  END DO
2436 
2437  END DO !iso
2438  END DO !ipgf
2439  END DO !iset
2440 
2441  !Do the integration in terms of matrix-matrix multiplications
2442  ALLOCATE (intso_h(nset*maxso, nset*maxso, nspins))
2443  ALLOCATE (intso_s(nset*maxso, nset*maxso, nspins))
2444  intso_h = 0.0_dp; intso_s = 0.0_dp
2445 
2446  ALLOCATE (fga(na, 1))
2447  ALLOCATE (fgr(nr, 1))
2448  ALLOCATE (work(na, 1))
2449  fga = 0.0_dp; fgr = 0.0_dp; work = 0.0_dp
2450 
2451  DO iset = 1, nset
2452  DO jset = 1, nset
2453  DO ipgf = 1, npgf(iset)
2454  starti = (iset - 1)*maxso + (ipgf - 1)*nsoset(lmax(iset))
2455  DO jpgf = 1, npgf(jset)
2456  startj = (jset - 1)*maxso + (jpgf - 1)*nsoset(lmax(jset))
2457 
2458  DO iso = nsoset(lmin(iset) - 1) + 1, nsoset(lmax(iset))
2459  DO jso = nsoset(lmin(jset) - 1) + 1, nsoset(lmax(jset))
2460 
2461  !Two component per function => 4 terms in total
2462 
2463  ! take 0.5*r1*a1(dir) * V * r1*a1(dir)
2464  fgr(1:nr, 1) = r1(1:nr, starti + iso)*r1(1:nr, startj + jso)
2465  DO dir = 1, 3
2466  fga(1:na, 1) = a1(1:na, starti + iso, dir)*a1(1:na, startj + jso, dir)
2467 
2468  DO ispin = 1, nspins
2469  CALL dgemm('N', 'N', na, 1, nr, 0.5_dp, vtau_h(:, :, ispin), na, fgr, &
2470  nr, 0.0_dp, work, na)
2471  CALL dgemm('T', 'N', 1, 1, na, 1.0_dp, work, na, fga, na, 1.0_dp, &
2472  intso_h(starti + iso:, startj + jso, ispin), 1)
2473 
2474  CALL dgemm('N', 'N', na, 1, nr, 0.5_dp, vtau_s(:, :, ispin), na, fgr, &
2475  nr, 0.0_dp, work, na)
2476  CALL dgemm('T', 'N', 1, 1, na, 1.0_dp, work, na, fga, na, 1.0_dp, &
2477  intso_s(starti + iso:, startj + jso, ispin), 1)
2478  END DO
2479  END DO !dir
2480 
2481  ! add 0.5*r1*a1(dir) * V * r2*a2(dir)
2482  fgr(1:nr, 1) = r1(1:nr, starti + iso)*r2(1:nr, startj + jso)
2483  DO dir = 1, 3
2484  fga(1:na, 1) = a1(1:na, starti + iso, dir)*a2(1:na, startj + jso, dir)
2485 
2486  DO ispin = 1, nspins
2487  CALL dgemm('N', 'N', na, 1, nr, 0.5_dp, vtau_h(:, :, ispin), na, fgr, &
2488  nr, 0.0_dp, work, na)
2489  CALL dgemm('T', 'N', 1, 1, na, 1.0_dp, work, na, fga, na, 1.0_dp, &
2490  intso_h(starti + iso:, startj + jso, ispin), 1)
2491 
2492  CALL dgemm('N', 'N', na, 1, nr, 0.5_dp, vtau_s(:, :, ispin), na, fgr, &
2493  nr, 0.0_dp, work, na)
2494  CALL dgemm('T', 'N', 1, 1, na, 1.0_dp, work, na, fga, na, 1.0_dp, &
2495  intso_s(starti + iso:, startj + jso, ispin), 1)
2496  END DO
2497  END DO !dir
2498 
2499  ! add 0.5*r2*a2(dir) * V * r1*a1(dir)
2500  fgr(1:nr, 1) = r2(1:nr, starti + iso)*r1(1:nr, startj + jso)
2501  DO dir = 1, 3
2502  fga(1:na, 1) = a2(1:na, starti + iso, dir)*a1(1:na, startj + jso, dir)
2503 
2504  DO ispin = 1, nspins
2505  CALL dgemm('N', 'N', na, 1, nr, 0.5_dp, vtau_h(:, :, ispin), na, fgr, &
2506  nr, 0.0_dp, work, na)
2507  CALL dgemm('T', 'N', 1, 1, na, 1.0_dp, work, na, fga, na, 1.0_dp, &
2508  intso_h(starti + iso:, startj + jso, ispin), 1)
2509 
2510  CALL dgemm('N', 'N', na, 1, nr, 0.5_dp, vtau_s(:, :, ispin), na, fgr, &
2511  nr, 0.0_dp, work, na)
2512  CALL dgemm('T', 'N', 1, 1, na, 1.0_dp, work, na, fga, na, 1.0_dp, &
2513  intso_s(starti + iso:, startj + jso, ispin), 1)
2514  END DO
2515  END DO !dir
2516 
2517  ! add the last term: 0.5*r2*a2(dir) * V * r2*a2(dir)
2518  fgr(1:nr, 1) = r2(1:nr, starti + iso)*r2(1:nr, startj + jso)
2519  DO dir = 1, 3
2520  fga(1:na, 1) = a2(1:na, starti + iso, dir)*a2(1:na, startj + jso, dir)
2521 
2522  DO ispin = 1, nspins
2523  CALL dgemm('N', 'N', na, 1, nr, 0.5_dp, vtau_h(:, :, ispin), na, fgr, &
2524  nr, 0.0_dp, work, na)
2525  CALL dgemm('T', 'N', 1, 1, na, 1.0_dp, work, na, fga, na, 1.0_dp, &
2526  intso_h(starti + iso:, startj + jso, ispin), 1)
2527 
2528  CALL dgemm('N', 'N', na, 1, nr, 0.5_dp, vtau_s(:, :, ispin), na, fgr, &
2529  nr, 0.0_dp, work, na)
2530  CALL dgemm('T', 'N', 1, 1, na, 1.0_dp, work, na, fga, na, 1.0_dp, &
2531  intso_s(starti + iso:, startj + jso, ispin), 1)
2532  END DO
2533  END DO !dir
2534 
2535  END DO !jso
2536  END DO !iso
2537 
2538  END DO !jpgf
2539  END DO !ipgf
2540  END DO !jset
2541  END DO !iset
2542 
2543  ! Put the integrals in the rho_atom data structure
2544  DO ispin = 1, nspins
2545  int_hh(ispin)%r_coef(:, :) = int_hh(ispin)%r_coef(:, :) + intso_h(:, :, ispin)
2546  int_ss(ispin)%r_coef(:, :) = int_ss(ispin)%r_coef(:, :) + intso_s(:, :, ispin)
2547  END DO
2548 
2549  CALL timestop(handle)
2550 
2551  END SUBROUTINE dgavtaudgb
2552 
2553 END MODULE qs_vxc_atom
dgemm
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
Definition: grid_cpu_task_list.c:195
atom
Definition: atom.F:9
atomic_kind_types
Define the atomic kind types and their sub types.
Definition: atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition: atomic_kind_types.F:138
basis_set_types
Definition: basis_set_types.F:15
basis_set_types::get_gto_basis_set
subroutine, public get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, lmin, lx, ly, lz, m, ncgf_set, npgf, nsgf_set, nshell, cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, last_cgf, last_sgf, n, gcc, maxco, maxl, maxpgf, maxsgf_set, maxshell, maxso, nco_sum, npgf_sum, nshell_sum, maxder, short_kind_radius)
...
Definition: basis_set_types.F:615
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition: cp_control_types.F:12
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition: input_constants.F:17
input_constants::tddfpt_excitations
integer, parameter, public tddfpt_excitations
Definition: input_constants.F:590
input_constants::tddfpt_triplet
integer, parameter, public tddfpt_triplet
Definition: input_constants.F:575
input_constants::xc_none
integer, parameter, public xc_none
Definition: input_constants.F:551
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition: input_section_types.F:15
input_section_types::section_vals_get_subs_vals
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
Definition: input_section_types.F:724
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition: input_section_types.F:1040
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
memory_utilities
Utility routines for the memory handling.
Definition: memory_utilities.F:14
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition: orbital_pointers.F:15
orbital_pointers::nsoset
integer, dimension(:), allocatable, public nsoset
Definition: orbital_pointers.F:42
orbital_pointers::indso
integer, dimension(:, :), allocatable, public indso
Definition: orbital_pointers.F:43
paw_basis_types
Definition: paw_basis_types.F:13
paw_basis_types::get_paw_basis_info
subroutine, public get_paw_basis_info(basis_1c, o2nindex, n2oindex, nsatbas)
Return some info on the PAW basis derived from a GTO basis set.
Definition: paw_basis_types.F:40
qs_environment_types
Definition: qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Definition: qs_environment_types.F:502
qs_grid_atom
Definition: qs_grid_atom.F:8
qs_harmonics_atom
Definition: qs_harmonics_atom.F:8
qs_kind_types
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
qs_kind_types::has_nlcc
logical function, public has_nlcc(qs_kind_set)
finds if a given qs run needs to use nlcc
Definition: qs_kind_types.F:3499
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_r3d_rs_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, U_of_dft_plus_u, J_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, J0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition: qs_kind_types.F:451
qs_linres_types
Type definitiona for linear response calculations.
Definition: qs_linres_types.F:12
qs_rho_atom_types
Definition: qs_rho_atom_types.F:8
qs_rho_atom_types::get_rho_atom
subroutine, public get_rho_atom(rho_atom, cpc_h, cpc_s, rho_rad_h, rho_rad_s, drho_rad_h, drho_rad_s, vrho_rad_h, vrho_rad_s, rho_rad_h_d, rho_rad_s_d, ga_Vlocal_gb_h, ga_Vlocal_gb_s, int_scr_h, int_scr_s)
...
Definition: qs_rho_atom_types.F:195
qs_vxc_atom
routines that build the integrals of the Vxc potential calculated for the atomic density in the basis...
Definition: qs_vxc_atom.F:12
qs_vxc_atom::calculate_gfxc_atom
subroutine, public calculate_gfxc_atom(qs_env, rho0_atom_set, rho1_atom_set, rho2_atom_set, kind_set, xc_section, is_triplet, accuracy)
...
Definition: qs_vxc_atom.F:755
qs_vxc_atom::calc_rho_angular
subroutine, public calc_rho_angular(grid_atom, harmonics, nspins, grad_func, ir, r_h, r_s, rho_h, rho_s, dr_h, dr_s, r_h_d, r_s_d, drho_h, drho_s)
...
Definition: qs_vxc_atom.F:1537
qs_vxc_atom::calculate_xc_2nd_deriv_atom
subroutine, public calculate_xc_2nd_deriv_atom(rho_atom_set, rho1_atom_set, qs_env, xc_section, para_env, do_tddft, do_tddfpt2, do_triplet, kind_set_external)
...
Definition: qs_vxc_atom.F:446
qs_vxc_atom::gavxcgb_nogc
subroutine, public gavxcgb_nogc(vxc_h, vxc_s, int_hh, int_ss, grid_atom, basis_1c, harmonics, nspins)
...
Definition: qs_vxc_atom.F:1954
qs_vxc_atom::calculate_vxc_atom
subroutine, public calculate_vxc_atom(qs_env, energy_only, exc1, gradient_atom_set, adiabatic_rescale_factor, kind_set_external, rho_atom_set_external, xc_section_external)
...
Definition: qs_vxc_atom.F:87
qs_vxc_atom::gfxc_atom_diff
subroutine, public gfxc_atom_diff(qs_env, rho0_atom_set, rho1_atom_set, rho2_atom_set, kind_set, xc_section, is_triplet, accuracy, epsrho)
...
Definition: qs_vxc_atom.F:1155
util
All kind of helpful little routines.
Definition: util.F:14
util::get_limit
pure integer function, dimension(2), public get_limit(m, n, me)
divide m entries into n parts, return size of part me
Definition: util.F:333
xc_atom
Definition: xc_atom.F:9
xc_atom::xc_rho_set_atom_update
subroutine, public xc_rho_set_atom_update(rho_set, needs, nspins, bo)
...
Definition: xc_atom.F:432
xc_atom::fill_rho_set
subroutine, public fill_rho_set(rho_set, lsd, nspins, needs, rho, drho, tau, na, ir)
...
Definition: xc_atom.F:592
xc_atom::xc_2nd_deriv_of_r
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
xc_atom::vxc_of_r_new
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
xc_derivative_set_types
represent a group ofunctional derivatives
Definition: xc_derivative_set_types.F:14
xc_derivative_set_types::xc_dset_zero_all
subroutine, public xc_dset_zero_all(deriv_set)
...
Definition: xc_derivative_set_types.F:167
xc_derivative_set_types::xc_dset_release
subroutine, public xc_dset_release(derivative_set)
releases a derivative set
Definition: xc_derivative_set_types.F:144
xc_derivative_set_types::xc_dset_create
subroutine, public xc_dset_create(derivative_set, pw_pool, local_bounds)
creates a derivative set object
Definition: xc_derivative_set_types.F:112
xc_derivatives
Definition: xc_derivatives.F:9
xc_derivatives::xc_functionals_get_needs
type(xc_rho_cflags_type) function, public xc_functionals_get_needs(functionals, lsd, calc_potential)
...
Definition: xc_derivatives.F:600
xc_rho_cflags_types
contains the structure
Definition: xc_rho_cflags_types.F:14
xc_rho_set_types
contains the structure
Definition: xc_rho_set_types.F:14
xc_rho_set_types::xc_rho_set_create
subroutine, public xc_rho_set_create(rho_set, local_bounds, rho_cutoff, drho_cutoff, tau_cutoff)
allocates and does (minimal) initialization of a rho_set
Definition: xc_rho_set_types.F:111
xc_rho_set_types::xc_rho_set_release
subroutine, public xc_rho_set_release(rho_set, pw_pool)
releases the given rho_set
Definition: xc_rho_set_types.F:129