(git:e7e05ae)
xas_tdp_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 This module deals with all the integrals done on local atomic grids in xas_tdp. This is
10 !> mostly used to compute the xc kernel matrix elements wrt two RI basis elements (centered
11 !> on the same excited atom) <P|fxc(r)|Q>, where the kernel fxc is purely a function of the
12 !> ground state density and r. This is also used to compute the SOC matrix elements in the
13 !> orbital basis
14 ! **************************************************************************************************
15 MODULE xas_tdp_atom
16  USE ai_contraction_sphi, ONLY: ab_contract
17  USE atom_operators, ONLY: calculate_model_potential
18  USE basis_set_types, ONLY: get_gto_basis_set,&
19  gto_basis_set_p_type,&
20  gto_basis_set_type
21  USE cell_types, ONLY: cell_type,&
22  pbc
23  USE cp_array_utils, ONLY: cp_1d_i_p_type,&
24  cp_1d_r_p_type,&
25  cp_2d_r_p_type,&
26  cp_3d_r_p_type
27  USE cp_blacs_env, ONLY: cp_blacs_env_type
28  USE cp_control_types, ONLY: dft_control_type,&
29  qs_control_type
30  USE cp_dbcsr_cholesky, ONLY: cp_dbcsr_cholesky_decompose,&
31  cp_dbcsr_cholesky_invert
32  USE cp_dbcsr_operations, ONLY: dbcsr_deallocate_matrix_set
33  USE cp_log_handling, ONLY: cp_logger_get_default_io_unit
34  USE dbcsr_api, ONLY: &
35  dbcsr_copy, dbcsr_create, dbcsr_distribution_get, dbcsr_distribution_new, &
36  dbcsr_distribution_release, dbcsr_distribution_type, dbcsr_filter, dbcsr_finalize, &
37  dbcsr_get_block_p, dbcsr_get_stored_coordinates, dbcsr_iterator_blocks_left, &
38  dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
39  dbcsr_p_type, dbcsr_put_block, dbcsr_release, dbcsr_replicate_all, dbcsr_set, dbcsr_type, &
40  dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, dbcsr_type_symmetric
41  USE dbt_api, ONLY: dbt_destroy,&
42  dbt_get_block,&
43  dbt_iterator_blocks_left,&
44  dbt_iterator_next_block,&
45  dbt_iterator_start,&
46  dbt_iterator_stop,&
47  dbt_iterator_type,&
48  dbt_type
49  USE input_constants, ONLY: do_potential_id
50  USE input_section_types, ONLY: section_vals_get_subs_vals,&
51  section_vals_type
52  USE kinds, ONLY: default_string_length,&
53  dp
54  USE lebedev, ONLY: deallocate_lebedev_grids,&
55  get_number_of_lebedev_grid,&
56  init_lebedev_grids,&
57  lebedev_grid
58  USE libint_2c_3c, ONLY: libint_potential_type
59  USE mathlib, ONLY: get_diag,&
60  invmat_symm
61  USE memory_utilities, ONLY: reallocate
62  USE message_passing, ONLY: mp_comm_type,&
63  mp_para_env_type,&
64  mp_request_type,&
65  mp_waitall
66  USE orbital_pointers, ONLY: indco,&
67  indso,&
68  nco,&
69  ncoset,&
70  nso,&
71  nsoset
72  USE orbital_transformation_matrices, ONLY: orbtramat
73  USE particle_methods, ONLY: get_particle_set
74  USE particle_types, ONLY: particle_type
75  USE physcon, ONLY: c_light_au
76  USE qs_environment_types, ONLY: get_qs_env,&
77  qs_environment_type
78  USE qs_grid_atom, ONLY: allocate_grid_atom,&
79  create_grid_atom,&
80  grid_atom_type
81  USE qs_harmonics_atom, ONLY: allocate_harmonics_atom,&
82  create_harmonics_atom,&
83  get_maxl_cg,&
84  get_none0_cg_list,&
85  harmonics_atom_type
86  USE qs_integral_utils, ONLY: basis_set_list_setup
87  USE qs_kind_types, ONLY: get_qs_kind,&
88  get_qs_kind_set,&
89  qs_kind_type
90  USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type,&
91  release_neighbor_list_sets
92  USE qs_overlap, ONLY: build_overlap_matrix_simple
93  USE qs_rho_types, ONLY: qs_rho_get,&
94  qs_rho_type
95  USE qs_tddfpt2_soc_types, ONLY: soc_atom_env_type
96  USE spherical_harmonics, ONLY: clebsch_gordon,&
97  clebsch_gordon_deallocate,&
98  clebsch_gordon_init
99  USE util, ONLY: get_limit,&
100  sort_unique
101  USE xas_tdp_integrals, ONLY: build_xas_tdp_3c_nl,&
102  build_xas_tdp_ovlp_nl,&
103  create_pqx_tensor,&
104  fill_pqx_tensor
105  USE xas_tdp_types, ONLY: batch_info_type,&
106  get_proc_batch_sizes,&
107  release_batch_info,&
108  xas_atom_env_type,&
109  xas_tdp_control_type,&
110  xas_tdp_env_type
111  USE xc, ONLY: divide_by_norm_drho
112  USE xc_atom, ONLY: xc_rho_set_atom_update
113  USE xc_derivative_desc, ONLY: deriv_norm_drho,&
114  deriv_norm_drhoa,&
115  deriv_norm_drhob,&
116  deriv_rhoa,&
117  deriv_rhob
118  USE xc_derivative_set_types, ONLY: xc_derivative_set_type,&
119  xc_dset_create,&
120  xc_dset_get_derivative,&
121  xc_dset_release
122  USE xc_derivative_types, ONLY: xc_derivative_get,&
123  xc_derivative_type
124  USE xc_derivatives, ONLY: xc_functionals_eval,&
125  xc_functionals_get_needs
126  USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
127  USE xc_rho_set_types, ONLY: xc_rho_set_create,&
128  xc_rho_set_release,&
129  xc_rho_set_type
130 
131 !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num
132 #include "./base/base_uses.f90"
133 
134  IMPLICIT NONE
135 
136  PRIVATE
137 
138  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xas_tdp_atom'
139 
140  PUBLIC :: init_xas_atom_env, &
141  integrate_fxc_atoms, &
142  integrate_soc_atoms, &
143  calculate_density_coeffs, &
144  compute_sphi_so, &
145  truncate_radial_grid, &
146  init_xas_atom_grid_harmo
147 
148 CONTAINS
149 
150 ! **************************************************************************************************
151 !> \brief Initializes a xas_atom_env type given the qs_enxas_atom_env, qs_envv
152 !> \param xas_atom_env the xas_atom_env to initialize
153 !> \param xas_tdp_env ...
154 !> \param xas_tdp_control ...
155 !> \param qs_env ...
156 !> \param ltddfpt ...
157 ! **************************************************************************************************
158  SUBROUTINE init_xas_atom_env(xas_atom_env, xas_tdp_env, xas_tdp_control, qs_env, ltddfpt)
159 
160  TYPE(xas_atom_env_type), POINTER :: xas_atom_env
161  TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
162  TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
163  TYPE(qs_environment_type), POINTER :: qs_env
164  LOGICAL, OPTIONAL :: ltddfpt
165 
166  CHARACTER(len=*), PARAMETER :: routinen = 'init_xas_atom_env'
167 
168  INTEGER :: handle, ikind, natom, nex_atoms, &
169  nex_kinds, nkind, nspins
170  LOGICAL :: do_xc
171  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
172  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
173 
174  NULLIFY (qs_kind_set, particle_set)
175 
176  CALL timeset(routinen, handle)
177 
178 ! Initializing the type
179  CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, natom=natom, particle_set=particle_set)
180 
181  nkind = SIZE(qs_kind_set)
182  nex_kinds = xas_tdp_env%nex_kinds
183  nex_atoms = xas_tdp_env%nex_atoms
184  do_xc = xas_tdp_control%do_xc
185  IF (PRESENT(ltddfpt) .AND. ltddfpt) do_xc = .false.
186  nspins = 1; IF (xas_tdp_control%do_uks .OR. xas_tdp_control%do_roks) nspins = 2
187 
188  xas_atom_env%nspins = nspins
189  xas_atom_env%ri_radius = xas_tdp_control%ri_radius
190  ALLOCATE (xas_atom_env%grid_atom_set(nkind))
191  ALLOCATE (xas_atom_env%harmonics_atom_set(nkind))
192  ALLOCATE (xas_atom_env%ri_sphi_so(nkind))
193  ALLOCATE (xas_atom_env%orb_sphi_so(nkind))
194  IF (do_xc) THEN
195  ALLOCATE (xas_atom_env%gr(nkind))
196  ALLOCATE (xas_atom_env%ga(nkind))
197  ALLOCATE (xas_atom_env%dgr1(nkind))
198  ALLOCATE (xas_atom_env%dgr2(nkind))
199  ALLOCATE (xas_atom_env%dga1(nkind))
200  ALLOCATE (xas_atom_env%dga2(nkind))
201  END IF
202 
203  xas_atom_env%excited_atoms => xas_tdp_env%ex_atom_indices
204  xas_atom_env%excited_kinds => xas_tdp_env%ex_kind_indices
205 
206 ! Allocate and initialize the atomic grids and harmonics
207  CALL init_xas_atom_grid_harmo(xas_atom_env, xas_tdp_control%grid_info, do_xc, qs_env)
208 
209 ! Compute the contraction coefficients for spherical orbitals
210  DO ikind = 1, nkind
211  NULLIFY (xas_atom_env%orb_sphi_so(ikind)%array, xas_atom_env%ri_sphi_so(ikind)%array)
212  CALL compute_sphi_so(ikind, "ORB", xas_atom_env%orb_sphi_so(ikind)%array, qs_env)
213  IF (any(xas_atom_env%excited_kinds == ikind)) THEN
214  CALL compute_sphi_so(ikind, "RI_XAS", xas_atom_env%ri_sphi_so(ikind)%array, qs_env)
215  END IF
216  END DO !ikind
217 
218 ! Compute the coefficients to expand the density in the RI_XAS basis, if requested
219  IF (do_xc .AND. (.NOT. xas_tdp_control%xps_only)) THEN
220  CALL calculate_density_coeffs(xas_atom_env=xas_atom_env, qs_env=qs_env)
221  END IF
222 
223  CALL timestop(handle)
224 
225  END SUBROUTINE init_xas_atom_env
226 
227 ! **************************************************************************************************
228 !> \brief Initializes the atomic grids and harmonics for the RI atomic calculations
229 !> \param xas_atom_env ...
230 !> \param grid_info ...
231 !> \param do_xc Whether the xc kernel will ne computed on the atomic grids. If not, the harmonics
232 !> are built for the orbital basis for all kinds.
233 !> \param qs_env ...
234 !> \note Largely inspired by init_rho_atom subroutine
235 ! **************************************************************************************************
236  SUBROUTINE init_xas_atom_grid_harmo(xas_atom_env, grid_info, do_xc, qs_env)
237 
238  TYPE(xas_atom_env_type), POINTER :: xas_atom_env
239  CHARACTER(len=default_string_length), &
240  DIMENSION(:, :), POINTER :: grid_info
241  LOGICAL, INTENT(IN) :: do_xc
242  TYPE(qs_environment_type), POINTER :: qs_env
243 
244  CHARACTER(LEN=default_string_length) :: kind_name
245  INTEGER :: igrid, ikind, il, iso, iso1, iso2, l1, l1l2, l2, la, lc1, lc2, lcleb, ll, llmax, &
246  lp, m1, m2, max_s_harm, max_s_set, maxl, maxlgto, maxs, mm, mp, na, nr, quadrature, stat
247  REAL(dp) :: kind_radius
248  REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: rga
249  REAL(dp), DIMENSION(:, :, :), POINTER :: my_cg
250  TYPE(dft_control_type), POINTER :: dft_control
251  TYPE(grid_atom_type), POINTER :: grid_atom
252  TYPE(gto_basis_set_type), POINTER :: tmp_basis
253  TYPE(harmonics_atom_type), POINTER :: harmonics
254  TYPE(qs_control_type), POINTER :: qs_control
255  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
256 
257  NULLIFY (my_cg, qs_kind_set, tmp_basis, grid_atom, harmonics, qs_control, dft_control)
258 
259 ! Initialization of some integer for the CG coeff generation
260  CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, dft_control=dft_control)
261  IF (do_xc) THEN
262  CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto, basis_type="RI_XAS")
263  ELSE
264  CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto, basis_type="ORB")
265  END IF
266  qs_control => dft_control%qs_control
267 
268  !maximum expansion
269  llmax = 2*maxlgto
270  max_s_harm = nsoset(llmax)
271  max_s_set = nsoset(maxlgto)
272  lcleb = llmax
273 
274 ! Allocate and compute the CG coeffs (copied from init_rho_atom)
275  CALL clebsch_gordon_init(lcleb)
276  CALL reallocate(my_cg, 1, max_s_set, 1, max_s_set, 1, max_s_harm)
277 
278  ALLOCATE (rga(lcleb, 2))
279  DO lc1 = 0, maxlgto
280  DO iso1 = nsoset(lc1 - 1) + 1, nsoset(lc1)
281  l1 = indso(1, iso1)
282  m1 = indso(2, iso1)
283  DO lc2 = 0, maxlgto
284  DO iso2 = nsoset(lc2 - 1) + 1, nsoset(lc2)
285  l2 = indso(1, iso2)
286  m2 = indso(2, iso2)
287  CALL clebsch_gordon(l1, m1, l2, m2, rga)
288  IF (l1 + l2 > llmax) THEN
289  l1l2 = llmax
290  ELSE
291  l1l2 = l1 + l2
292  END IF
293  mp = m1 + m2
294  mm = m1 - m2
295  IF (m1*m2 < 0 .OR. (m1*m2 == 0 .AND. (m1 < 0 .OR. m2 < 0))) THEN
296  mp = -abs(mp)
297  mm = -abs(mm)
298  ELSE
299  mp = abs(mp)
300  mm = abs(mm)
301  END IF
302  DO lp = mod(l1 + l2, 2), l1l2, 2
303  il = lp/2 + 1
304  IF (abs(mp) <= lp) THEN
305  IF (mp >= 0) THEN
306  iso = nsoset(lp - 1) + lp + 1 + mp
307  ELSE
308  iso = nsoset(lp - 1) + lp + 1 - abs(mp)
309  END IF
310  my_cg(iso1, iso2, iso) = rga(il, 1)
311  END IF
312  IF (mp /= mm .AND. abs(mm) <= lp) THEN
313  IF (mm >= 0) THEN
314  iso = nsoset(lp - 1) + lp + 1 + mm
315  ELSE
316  iso = nsoset(lp - 1) + lp + 1 - abs(mm)
317  END IF
318  my_cg(iso1, iso2, iso) = rga(il, 2)
319  END IF
320  END DO
321  END DO ! iso2
322  END DO ! lc2
323  END DO ! iso1
324  END DO ! lc1
325  DEALLOCATE (rga)
326  CALL clebsch_gordon_deallocate()
327 
328 ! Create the Lebedev grids and compute the spherical harmonics
329  CALL init_lebedev_grids()
330  quadrature = qs_control%gapw_control%quadrature
331 
332  DO ikind = 1, SIZE(xas_atom_env%grid_atom_set)
333 
334 ! Allocate the grid and the harmonics for this kind
335  NULLIFY (xas_atom_env%grid_atom_set(ikind)%grid_atom)
336  NULLIFY (xas_atom_env%harmonics_atom_set(ikind)%harmonics_atom)
337  CALL allocate_grid_atom(xas_atom_env%grid_atom_set(ikind)%grid_atom)
338  CALL allocate_harmonics_atom(xas_atom_env%harmonics_atom_set(ikind)%harmonics_atom)
339 
340  NULLIFY (grid_atom, harmonics)
341  grid_atom => xas_atom_env%grid_atom_set(ikind)%grid_atom
342  harmonics => xas_atom_env%harmonics_atom_set(ikind)%harmonics_atom
343 
344 ! Initialize some integers
345  CALL get_qs_kind(qs_kind_set(ikind), ngrid_rad=nr, ngrid_ang=na, name=kind_name)
346 
347  !take the grid dimension given as input, if none, take the GAPW ones above
348  IF (any(grid_info == kind_name)) THEN
349  DO igrid = 1, SIZE(grid_info, 1)
350  IF (grid_info(igrid, 1) == kind_name) THEN
351 
352  !hack to convert string into integer
353  READ (grid_info(igrid, 2), *, iostat=stat) na
354  IF (stat .NE. 0) cpabort("The 'na' value for the GRID keyword must be an integer")
355  READ (grid_info(igrid, 3), *, iostat=stat) nr
356  IF (stat .NE. 0) cpabort("The 'nr' value for the GRID keyword must be an integer")
357  EXIT
358  END IF
359  END DO
360  ELSE IF (do_xc .AND. any(xas_atom_env%excited_kinds == ikind)) THEN
361  !need good integration grids for the xc kernel, but taking the default GAPW grid
362  cpwarn("The default (and possibly small) GAPW grid is being used for one excited KIND")
363  END IF
364 
365  ll = get_number_of_lebedev_grid(n=na)
366  na = lebedev_grid(ll)%n
367  la = lebedev_grid(ll)%l
368  grid_atom%ng_sphere = na
369  grid_atom%nr = nr
370 
371 ! If this is an excited kind, create the harmonics with the RI_XAS basis, otherwise the ORB
372  IF (any(xas_atom_env%excited_kinds == ikind) .AND. do_xc) THEN
373  CALL get_qs_kind(qs_kind_set(ikind), basis_set=tmp_basis, basis_type="RI_XAS")
374  ELSE
375  CALL get_qs_kind(qs_kind_set(ikind), basis_set=tmp_basis, basis_type="ORB")
376  END IF
377  CALL get_gto_basis_set(gto_basis_set=tmp_basis, maxl=maxl, kind_radius=kind_radius)
378 
379  CALL create_grid_atom(grid_atom, nr, na, llmax, ll, quadrature)
380  CALL truncate_radial_grid(grid_atom, kind_radius)
381 
382  maxs = nsoset(maxl)
383  CALL create_harmonics_atom(harmonics, &
384  my_cg, na, llmax, maxs, max_s_harm, ll, grid_atom%wa, &
385  grid_atom%azi, grid_atom%pol)
386  CALL get_maxl_cg(harmonics, tmp_basis, llmax, max_s_harm)
387  END DO
388 
389  CALL deallocate_lebedev_grids()
390  DEALLOCATE (my_cg)
391 
392  END SUBROUTINE init_xas_atom_grid_harmo
393 
394 ! **************************************************************************************************
395 !> \brief Reduces the radial extension of an atomic grid such that it only covers a given radius
396 !> \param grid_atom the atomic grid from which we truncate the radial part
397 !> \param max_radius the maximal radial extension of the resulting grid
398 !> \note Since the RI density used for <P|fxc|Q> is only of quality close to the atom, and the
399 !> integrand only non-zero within the radius of the gaussian P,Q. One can reduce the grid
400 !> extansion to the largest radius of the RI basis set
401 ! **************************************************************************************************
402  SUBROUTINE truncate_radial_grid(grid_atom, max_radius)
403 
404  TYPE(grid_atom_type), POINTER :: grid_atom
405  REAL(dp), INTENT(IN) :: max_radius
406 
407  INTEGER :: first_ir, ir, llmax_p1, na, new_nr, nr
408 
409  nr = grid_atom%nr
410  na = grid_atom%ng_sphere
411  llmax_p1 = SIZE(grid_atom%rad2l, 2) - 1
412 
413 ! Find the index corresponding to the limiting radius (small ir => large radius)
414  DO ir = 1, nr
415  IF (grid_atom%rad(ir) < max_radius) THEN
416  first_ir = ir
417  EXIT
418  END IF
419  END DO
420  new_nr = nr - first_ir + 1
421 
422 ! Reallcoate everything that depends on r
423  grid_atom%nr = new_nr
424 
425  grid_atom%rad(1:new_nr) = grid_atom%rad(first_ir:nr)
426  grid_atom%rad2(1:new_nr) = grid_atom%rad2(first_ir:nr)
427  grid_atom%wr(1:new_nr) = grid_atom%wr(first_ir:nr)
428  grid_atom%weight(:, 1:new_nr) = grid_atom%weight(:, first_ir:nr)
429  grid_atom%rad2l(1:new_nr, :) = grid_atom%rad2l(first_ir:nr, :)
430  grid_atom%oorad2l(1:new_nr, :) = grid_atom%oorad2l(first_ir:nr, :)
431 
432  CALL reallocate(grid_atom%rad, 1, new_nr)
433  CALL reallocate(grid_atom%rad2, 1, new_nr)
434  CALL reallocate(grid_atom%wr, 1, new_nr)
435  CALL reallocate(grid_atom%weight, 1, na, 1, new_nr)
436  CALL reallocate(grid_atom%rad2l, 1, new_nr, 0, llmax_p1)
437  CALL reallocate(grid_atom%oorad2l, 1, new_nr, 0, llmax_p1)
438 
439  END SUBROUTINE truncate_radial_grid
440 
441 ! **************************************************************************************************
442 !> \brief Computes the contraction coefficients to go from spherical orbitals to sgf for a given
443 !> atomic kind and a given basis type.
444 !> \param ikind the kind for which we compute the coefficients
445 !> \param basis_type the type of basis for which we compute
446 !> \param sphi_so where the new contraction coefficients are stored (not yet allocated)
447 !> \param qs_env ...
448 ! **************************************************************************************************
449  SUBROUTINE compute_sphi_so(ikind, basis_type, sphi_so, qs_env)
450 
451  INTEGER, INTENT(IN) :: ikind
452  CHARACTER(len=*), INTENT(IN) :: basis_type
453  REAL(dp), DIMENSION(:, :), POINTER :: sphi_so
454  TYPE(qs_environment_type), POINTER :: qs_env
455 
456  INTEGER :: ico, ipgf, iset, iso, l, lx, ly, lz, &
457  maxso, nset, sgfi, start_c, start_s
458  INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf, nsgf_set
459  INTEGER, DIMENSION(:, :), POINTER :: first_sgf
460  REAL(dp) :: factor
461  REAL(dp), DIMENSION(:, :), POINTER :: sphi
462  TYPE(gto_basis_set_type), POINTER :: basis
463  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
464 
465  NULLIFY (basis, lmax, lmin, npgf, nsgf_set, qs_kind_set, first_sgf, sphi)
466 
467  CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set)
468  CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis, basis_type=basis_type)
469  CALL get_gto_basis_set(basis, lmax=lmax, nset=nset, npgf=npgf, maxso=maxso, lmin=lmin, &
470  nsgf_set=nsgf_set, sphi=sphi, first_sgf=first_sgf)
471 
472  ALLOCATE (sphi_so(maxso, sum(nsgf_set)))
473  sphi_so = 0.0_dp
474 
475  DO iset = 1, nset
476  sgfi = first_sgf(1, iset)
477  DO ipgf = 1, npgf(iset)
478  start_s = (ipgf - 1)*nsoset(lmax(iset))
479  start_c = (ipgf - 1)*ncoset(lmax(iset))
480  DO l = lmin(iset), lmax(iset)
481  DO iso = 1, nso(l)
482  DO ico = 1, nco(l)
483  lx = indco(1, ico + ncoset(l - 1))
484  ly = indco(2, ico + ncoset(l - 1))
485  lz = indco(3, ico + ncoset(l - 1))
486 !MK factor = orbtramat(l)%s2c(iso, ico) &
487 !MK *SQRT(4.0_dp*pi/dfac(2*l + 1)*dfac(2*lx - 1)*dfac(2*ly - 1)*dfac(2*lz - 1))
488  factor = orbtramat(l)%slm_inv(iso, ico)
489  sphi_so(start_s + nsoset(l - 1) + iso, sgfi:sgfi + nsgf_set(iset) - 1) = &
490  sphi_so(start_s + nsoset(l - 1) + iso, sgfi:sgfi + nsgf_set(iset) - 1) + &
491  factor*sphi(start_c + ncoset(l - 1) + ico, sgfi:sgfi + nsgf_set(iset) - 1)
492  END DO ! ico
493  END DO ! iso
494  END DO ! l
495  END DO ! ipgf
496  END DO ! iset
497 
498  END SUBROUTINE compute_sphi_so
499 
500 ! **************************************************************************************************
501 !> \brief Find the neighbors of a given set of atoms based on the non-zero blocks of a provided
502 !> overlap matrix. Optionally returns an array containing the indices of all involved atoms
503 !> (the given subset plus all their neighbors, without repetition) AND/OR an array of arrays
504 !> providing the indices of the neighbors of each input atom.
505 !> \param base_atoms the set of atoms for which we search neighbors
506 !> \param mat_s the overlap matrix used to find neighbors
507 !> \param radius the cutoff radius after which atoms are not considered neighbors
508 !> \param qs_env ...
509 !> \param all_neighbors the array uniquely contatining all indices of all atoms involved
510 !> \param neighbor_set array of arrays containing the neighbors of all given atoms
511 ! **************************************************************************************************
512  SUBROUTINE find_neighbors(base_atoms, mat_s, radius, qs_env, all_neighbors, neighbor_set)
513 
514  INTEGER, DIMENSION(:), INTENT(INOUT) :: base_atoms
515  TYPE(dbcsr_type), INTENT(IN) :: mat_s
516  REAL(dp) :: radius
517  TYPE(qs_environment_type), POINTER :: qs_env
518  INTEGER, DIMENSION(:), OPTIONAL, POINTER :: all_neighbors
519  TYPE(cp_1d_i_p_type), DIMENSION(:), OPTIONAL, &
520  POINTER :: neighbor_set
521 
522  INTEGER :: blk, i, iat, ibase, iblk, jblk, mepos, &
523  natom, nb, nbase
524  INTEGER, ALLOCATABLE, DIMENSION(:) :: blk_to_base, inb, who_is_there
525  INTEGER, ALLOCATABLE, DIMENSION(:, :) :: n_neighbors
526  LOGICAL, ALLOCATABLE, DIMENSION(:) :: is_base_atom
527  REAL(dp) :: dist2, rad2, ri(3), rij(3), rj(3)
528  TYPE(cell_type), POINTER :: cell
529  TYPE(cp_1d_i_p_type), DIMENSION(:), POINTER :: my_neighbor_set
530  TYPE(dbcsr_iterator_type) :: iter
531  TYPE(mp_para_env_type), POINTER :: para_env
532  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
533 
534  NULLIFY (particle_set, para_env, my_neighbor_set, cell)
535 
536  ! Initialization
537  CALL get_qs_env(qs_env, para_env=para_env, natom=natom, particle_set=particle_set, cell=cell)
538  mepos = para_env%mepos
539  nbase = SIZE(base_atoms)
540  !work with the neighbor_set structure, see at the end if we keep it
541  ALLOCATE (my_neighbor_set(nbase))
542  rad2 = radius**2
543 
544  ALLOCATE (blk_to_base(natom), is_base_atom(natom))
545  blk_to_base = 0; is_base_atom = .false.
546  DO ibase = 1, nbase
547  blk_to_base(base_atoms(ibase)) = ibase
548  is_base_atom(base_atoms(ibase)) = .true.
549  END DO
550 
551  ! First loop over S => count the number of neighbors
552  ALLOCATE (n_neighbors(nbase, 0:para_env%num_pe - 1))
553  n_neighbors = 0
554 
555  CALL dbcsr_iterator_start(iter, mat_s)
556  DO WHILE (dbcsr_iterator_blocks_left(iter))
557 
558  CALL dbcsr_iterator_next_block(iter, row=iblk, column=jblk, blk=blk)
559 
560  !avoid self-neighbors
561  IF (iblk == jblk) cycle
562 
563  !test distance
564  ri = pbc(particle_set(iblk)%r, cell)
565  rj = pbc(particle_set(jblk)%r, cell)
566  rij = pbc(ri, rj, cell)
567  dist2 = sum(rij**2)
568  IF (dist2 > rad2) cycle
569 
570  IF (is_base_atom(iblk)) THEN
571  ibase = blk_to_base(iblk)
572  n_neighbors(ibase, mepos) = n_neighbors(ibase, mepos) + 1
573  END IF
574  IF (is_base_atom(jblk)) THEN
575  ibase = blk_to_base(jblk)
576  n_neighbors(ibase, mepos) = n_neighbors(ibase, mepos) + 1
577  END IF
578 
579  END DO !iter
580  CALL dbcsr_iterator_stop(iter)
581  CALL para_env%sum(n_neighbors)
582 
583  ! Allocate the neighbor_set arrays at the correct length
584  DO ibase = 1, nbase
585  ALLOCATE (my_neighbor_set(ibase)%array(sum(n_neighbors(ibase, :))))
586  my_neighbor_set(ibase)%array = 0
587  END DO
588 
589  ! Loop a second time over S, this time fill the neighbors details
590  CALL dbcsr_iterator_start(iter, mat_s)
591  ALLOCATE (inb(nbase))
592  inb = 1
593  DO WHILE (dbcsr_iterator_blocks_left(iter))
594 
595  CALL dbcsr_iterator_next_block(iter, row=iblk, column=jblk, blk=blk)
596  IF (iblk == jblk) cycle
597 
598  !test distance
599  ri = pbc(particle_set(iblk)%r, cell)
600  rj = pbc(particle_set(jblk)%r, cell)
601  rij = pbc(ri, rj, cell)
602  dist2 = sum(rij**2)
603  IF (dist2 > rad2) cycle
604 
605  IF (is_base_atom(iblk)) THEN
606  ibase = blk_to_base(iblk)
607  my_neighbor_set(ibase)%array(sum(n_neighbors(ibase, 0:mepos - 1)) + inb(ibase)) = jblk
608  inb(ibase) = inb(ibase) + 1
609  END IF
610  IF (is_base_atom(jblk)) THEN
611  ibase = blk_to_base(jblk)
612  my_neighbor_set(ibase)%array(sum(n_neighbors(ibase, 0:mepos - 1)) + inb(ibase)) = iblk
613  inb(ibase) = inb(ibase) + 1
614  END IF
615 
616  END DO !iter
617  CALL dbcsr_iterator_stop(iter)
618 
619  ! Make sure the info is shared among the procs
620  DO ibase = 1, nbase
621  CALL para_env%sum(my_neighbor_set(ibase)%array)
622  END DO
623 
624  ! Gather all indices if asked for it
625  IF (PRESENT(all_neighbors)) THEN
626  ALLOCATE (who_is_there(natom))
627  who_is_there = 0
628 
629  DO ibase = 1, nbase
630  who_is_there(base_atoms(ibase)) = 1
631  DO nb = 1, SIZE(my_neighbor_set(ibase)%array)
632  who_is_there(my_neighbor_set(ibase)%array(nb)) = 1
633  END DO
634  END DO
635 
636  ALLOCATE (all_neighbors(natom))
637  i = 0
638  DO iat = 1, natom
639  IF (who_is_there(iat) == 1) THEN
640  i = i + 1
641  all_neighbors(i) = iat
642  END IF
643  END DO
644  CALL reallocate(all_neighbors, 1, i)
645  END IF
646 
647  ! If not asked for the neighbor set, deallocate it
648  IF (PRESENT(neighbor_set)) THEN
649  neighbor_set => my_neighbor_set
650  ELSE
651  DO ibase = 1, nbase
652  DEALLOCATE (my_neighbor_set(ibase)%array)
653  END DO
654  DEALLOCATE (my_neighbor_set)
655  END IF
656 
657  END SUBROUTINE find_neighbors
658 
659 ! **************************************************************************************************
660 !> \brief Returns the RI inverse overlap for a subset of the RI_XAS matrix spaning a given
661 !> excited atom and its neighbors.
662 !> \param ri_sinv the inverse overlap as a dbcsr matrix
663 !> \param whole_s the whole RI overlap matrix
664 !> \param neighbors the indeces of the excited atom and their neighbors
665 !> \param idx_to_nb array telling where any atom can be found in neighbors (if there at all)
666 !> \param basis_set_ri the RI basis set list for all kinds
667 !> \param qs_env ...
668 !> \note It is assumed that the neighbors are sorted, the output matrix is assumed to be small and
669 !> is replicated on all processors
670 ! **************************************************************************************************
671  SUBROUTINE get_exat_ri_sinv(ri_sinv, whole_s, neighbors, idx_to_nb, basis_set_ri, qs_env)
672 
673  TYPE(dbcsr_type) :: ri_sinv, whole_s
674  INTEGER, DIMENSION(:), INTENT(IN) :: neighbors, idx_to_nb
675  TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_ri
676  TYPE(qs_environment_type), POINTER :: qs_env
677 
678  CHARACTER(len=*), PARAMETER :: routinen = 'get_exat_ri_sinv'
679 
680  INTEGER :: blk, dest, group_handle, handle, iat, &
681  ikind, inb, ir, is, jat, jnb, natom, &
682  nnb, npcols, nprows, source, tag
683  INTEGER, DIMENSION(:), POINTER :: col_dist, nsgf, row_dist
684  INTEGER, DIMENSION(:, :), POINTER :: pgrid
685  LOGICAL :: found_risinv, found_whole
686  LOGICAL, ALLOCATABLE, DIMENSION(:) :: is_neighbor
687  REAL(dp) :: ri(3), rij(3), rj(3)
688  REAL(dp), ALLOCATABLE, DIMENSION(:) :: radius
689  REAL(dp), DIMENSION(:, :), POINTER :: block_risinv, block_whole
690  TYPE(cell_type), POINTER :: cell
691  TYPE(cp_2d_r_p_type), ALLOCATABLE, DIMENSION(:) :: recv_buff, send_buff
692  TYPE(cp_blacs_env_type), POINTER :: blacs_env
693  TYPE(dbcsr_distribution_type) :: sinv_dist
694  TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
695  TYPE(dbcsr_iterator_type) :: iter
696  TYPE(mp_comm_type) :: group
697  TYPE(mp_para_env_type), POINTER :: para_env
698  TYPE(mp_request_type), ALLOCATABLE, DIMENSION(:) :: recv_req, send_req
699  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
700 
701  NULLIFY (pgrid, dbcsr_dist, row_dist, col_dist, nsgf, particle_set, block_whole, block_risinv)
702  NULLIFY (cell, para_env, blacs_env)
703 
704  CALL timeset(routinen, handle)
705 
706  CALL get_qs_env(qs_env, dbcsr_dist=dbcsr_dist, particle_set=particle_set, natom=natom, &
707  para_env=para_env, blacs_env=blacs_env, cell=cell)
708  nnb = SIZE(neighbors)
709  ALLOCATE (nsgf(nnb), is_neighbor(natom), radius(nnb))
710  is_neighbor = .false.
711  DO inb = 1, nnb
712  iat = neighbors(inb)
713  ikind = particle_set(iat)%atomic_kind%kind_number
714  CALL get_gto_basis_set(basis_set_ri(ikind)%gto_basis_set, nsgf=nsgf(inb), kind_radius=radius(inb))
715  is_neighbor(iat) = .true.
716  END DO
717 
718  !Create the ri_sinv matrix based on some arbitrary dbcsr_dist
719  CALL dbcsr_distribution_get(dbcsr_dist, group=group_handle, pgrid=pgrid, nprows=nprows, npcols=npcols)
720  CALL group%set_handle(group_handle)
721 
722  ALLOCATE (row_dist(nnb), col_dist(nnb))
723  DO inb = 1, nnb
724  row_dist(inb) = modulo(nprows - inb, nprows)
725  col_dist(inb) = modulo(npcols - inb, npcols)
726  END DO
727 
728  CALL dbcsr_distribution_new(sinv_dist, group=group_handle, pgrid=pgrid, row_dist=row_dist, &
729  col_dist=col_dist)
730 
731  CALL dbcsr_create(matrix=ri_sinv, name="RI_SINV", matrix_type=dbcsr_type_symmetric, &
732  dist=sinv_dist, row_blk_size=nsgf, col_blk_size=nsgf)
733  !reserving the blocks in the correct pattern
734  DO inb = 1, nnb
735  ri = pbc(particle_set(neighbors(inb))%r, cell)
736  DO jnb = inb, nnb
737 
738  !do the atom overlap ?
739  rj = pbc(particle_set(neighbors(jnb))%r, cell)
740  rij = pbc(ri, rj, cell)
741  IF (sum(rij**2) > (radius(inb) + radius(jnb))**2) cycle
742 
743  CALL dbcsr_get_stored_coordinates(ri_sinv, inb, jnb, blk)
744  IF (para_env%mepos == blk) THEN
745  ALLOCATE (block_risinv(nsgf(inb), nsgf(jnb)))
746  block_risinv = 0.0_dp
747  CALL dbcsr_put_block(ri_sinv, inb, jnb, block_risinv)
748  DEALLOCATE (block_risinv)
749  END IF
750  END DO
751  END DO
752  CALL dbcsr_finalize(ri_sinv)
753 
754  CALL dbcsr_distribution_release(sinv_dist)
755  DEALLOCATE (row_dist, col_dist)
756 
757  !prepare the send and recv buffers we will need for change of dist between the two matrices
758  !worst case scenario: all neighbors are on same procs => need to send nnb**2 messages
759  ALLOCATE (send_buff(nnb**2), recv_buff(nnb**2))
760  ALLOCATE (send_req(nnb**2), recv_req(nnb**2))
761  is = 0; ir = 0
762 
763  !Loop over the whole RI overlap matrix and pick the blocks we need
764  CALL dbcsr_iterator_start(iter, whole_s)
765  DO WHILE (dbcsr_iterator_blocks_left(iter))
766 
767  CALL dbcsr_iterator_next_block(iter, row=iat, column=jat, blk=blk)
768  CALL dbcsr_get_block_p(whole_s, iat, jat, block_whole, found_whole)
769 
770  !only interested in neighbors
771  IF (.NOT. found_whole) cycle
772  IF (.NOT. is_neighbor(iat)) cycle
773  IF (.NOT. is_neighbor(jat)) cycle
774 
775  inb = idx_to_nb(iat)
776  jnb = idx_to_nb(jat)
777 
778  !If blocks are on the same proc for both matrices, simply copy
779  CALL dbcsr_get_block_p(ri_sinv, inb, jnb, block_risinv, found_risinv)
780  IF (found_risinv) THEN
781  CALL dcopy(nsgf(inb)*nsgf(jnb), block_whole, 1, block_risinv, 1)
782  ELSE
783 
784  !send the block with unique tag to the proc where inb,jnb is in ri_sinv
785  CALL dbcsr_get_stored_coordinates(ri_sinv, inb, jnb, dest)
786  is = is + 1
787  send_buff(is)%array => block_whole
788  tag = natom*inb + jnb
789  CALL group%isend(msgin=send_buff(is)%array, dest=dest, request=send_req(is), tag=tag)
790 
791  END IF
792 
793  END DO !dbcsr iter
794  CALL dbcsr_iterator_stop(iter)
795 
796  !Loop over ri_sinv and receive all those blocks
797  CALL dbcsr_iterator_start(iter, ri_sinv)
798  DO WHILE (dbcsr_iterator_blocks_left(iter))
799 
800  CALL dbcsr_iterator_next_block(iter, row=inb, column=jnb, blk=blk)
801  CALL dbcsr_get_block_p(ri_sinv, inb, jnb, block_risinv, found_risinv)
802 
803  IF (.NOT. found_risinv) cycle
804 
805  iat = neighbors(inb)
806  jat = neighbors(jnb)
807 
808  !If blocks are on the same proc on both matrices do nothing
809  CALL dbcsr_get_stored_coordinates(whole_s, iat, jat, source)
810  IF (para_env%mepos == source) cycle
811 
812  tag = natom*inb + jnb
813  ir = ir + 1
814  recv_buff(ir)%array => block_risinv
815  CALL group%irecv(msgout=recv_buff(ir)%array, source=source, request=recv_req(ir), &
816  tag=tag)
817 
818  END DO
819  CALL dbcsr_iterator_stop(iter)
820 
821  !make sure that all comm is over before proceeding
822  CALL mp_waitall(send_req(1:is))
823  CALL mp_waitall(recv_req(1:ir))
824 
825  !Invert. 2 cases: with or without neighbors. If no neighbors, easier to invert on one proc and
826  !avoid the whole fm to dbcsr to fm that is quite expensive
827  IF (nnb == 1) THEN
828 
829  CALL dbcsr_get_block_p(ri_sinv, 1, 1, block_risinv, found_risinv)
830  IF (found_risinv) THEN
831  CALL invmat_symm(block_risinv)
832  END IF
833 
834  ELSE
835  CALL cp_dbcsr_cholesky_decompose(ri_sinv, para_env=para_env, blacs_env=blacs_env)
836  CALL cp_dbcsr_cholesky_invert(ri_sinv, para_env=para_env, blacs_env=blacs_env, upper_to_full=.true.)
837  CALL dbcsr_filter(ri_sinv, 1.e-10_dp) !make sure ri_sinv is sparse coming out of fm routines
838  END IF
839  CALL dbcsr_replicate_all(ri_sinv)
840 
841  !clean-up
842  DEALLOCATE (nsgf)
843 
844  CALL timestop(handle)
845 
846  END SUBROUTINE get_exat_ri_sinv
847 
848 ! **************************************************************************************************
849 !> \brief Compute the coefficients to project the density on a partial RI_XAS basis
850 !> \param xas_atom_env ...
851 !> \param qs_env ...
852 !> \note The density is n = sum_ab P_ab*phi_a*phi_b, the RI basis covers the products of orbital sgfs
853 !> => n = sum_ab sum_cd P_ab (phi_a phi_b xi_c) S_cd^-1 xi_d
854 !> = sum_d coeff_d xi_d , where xi are the RI basis func.
855 !> In this case, with the partial RI projection, the RI basis is restricted to an excited atom
856 !> and its neighbors at a time. Leads to smaller overlap matrix to invert and less 3-center
857 !> overlap to compute. The procedure is repeated for each excited atom
858 ! **************************************************************************************************
859  SUBROUTINE calculate_density_coeffs(xas_atom_env, qs_env)
860 
861  TYPE(xas_atom_env_type), POINTER :: xas_atom_env
862  TYPE(qs_environment_type), POINTER :: qs_env
863 
864  CHARACTER(len=*), PARAMETER :: routinen = 'calculate_density_coeffs'
865 
866  INTEGER :: exat, handle, i, iat, iatom, iex, inb, &
867  ind(3), ispin, jatom, jnb, katom, &
868  natom, nex, nkind, nnb, nspins, &
869  output_unit, ri_at
870  INTEGER, ALLOCATABLE, DIMENSION(:) :: blk_size_orb, blk_size_ri, idx_to_nb, &
871  neighbors
872  INTEGER, DIMENSION(:), POINTER :: all_ri_atoms
873  LOGICAL :: pmat_found, pmat_foundt, sinv_found, &
874  sinv_foundt, tensor_found, unique
875  REAL(dp) :: factor, prefac
876  REAL(dp), ALLOCATABLE, DIMENSION(:) :: work2
877  REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: work1
878  REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: t_block
879  REAL(dp), DIMENSION(:, :), POINTER :: pmat_block, pmat_blockt, sinv_block, &
880  sinv_blockt
881  TYPE(cp_blacs_env_type), POINTER :: blacs_env
882  TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
883  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: overlap, rho_ao
884  TYPE(dbcsr_type) :: ri_sinv
885  TYPE(dbcsr_type), POINTER :: ri_mats
886  TYPE(dbt_iterator_type) :: iter
887  TYPE(dbt_type) :: pqx
888  TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_orb, basis_set_ri
889  TYPE(libint_potential_type) :: pot
890  TYPE(mp_para_env_type), POINTER :: para_env
891  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
892  POINTER :: ab_list, ac_list, sab_ri
893  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
894  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
895  TYPE(qs_rho_type), POINTER :: rho
896 
897  NULLIFY (qs_kind_set, basis_set_ri, basis_set_orb, ac_list, rho, rho_ao, sab_ri, ri_mats)
898  NULLIFY (particle_set, para_env, all_ri_atoms, overlap, pmat_blockt)
899  NULLIFY (blacs_env, pmat_block, ab_list, dbcsr_dist, sinv_block, sinv_blockt)
900 
901  !Idea: We don't do a full RI here as it would be too expensive in many ways (inversion of a
902  ! large matrix, many 3-center overlaps, density coefficients for all atoms, etc...)
903  ! Instead, we go excited atom by excited atom and only do a RI expansion on its basis
904  ! and that of its closest neighbors (defined through RI_RADIUS), such that we only have
905  ! very small matrices to invert and only a few (abc) overlp integrals with c on the
906  ! excited atom its neighbors. This is physically sound since we only need the density
907  ! well defined on the excited atom grid as we do (aI|P)*(P|Q)^-1*(Q|fxc|R)*(R|S)^-1*(S|Jb)
908 
909  CALL timeset(routinen, handle)
910 
911  CALL get_qs_env(qs_env, nkind=nkind, qs_kind_set=qs_kind_set, rho=rho, &
912  natom=natom, particle_set=particle_set, dbcsr_dist=dbcsr_dist, &
913  para_env=para_env, blacs_env=blacs_env)
914  nspins = xas_atom_env%nspins
915  nex = SIZE(xas_atom_env%excited_atoms)
916  CALL qs_rho_get(rho, rho_ao=rho_ao)
917 
918 ! Create the needed neighbor list and basis set lists.
919  ALLOCATE (basis_set_ri(nkind))
920  ALLOCATE (basis_set_orb(nkind))
921  CALL basis_set_list_setup(basis_set_ri, "RI_XAS", qs_kind_set)
922  CALL basis_set_list_setup(basis_set_orb, "ORB", qs_kind_set)
923 
924 ! Compute the RI overlap matrix on the whole system
925  CALL build_xas_tdp_ovlp_nl(sab_ri, basis_set_ri, basis_set_ri, qs_env)
926  CALL build_overlap_matrix_simple(qs_env%ks_env, overlap, basis_set_ri, basis_set_ri, sab_ri)
927  ri_mats => overlap(1)%matrix
928  CALL release_neighbor_list_sets(sab_ri)
929 
930 ! Get the neighbors of the excited atoms (= all the atoms where density coeffs are needed)
931  CALL find_neighbors(xas_atom_env%excited_atoms, ri_mats, xas_atom_env%ri_radius, &
932  qs_env, all_neighbors=all_ri_atoms, neighbor_set=xas_atom_env%exat_neighbors)
933 
934  !keep in mind that double occupation is included in rho_ao in case of closed-shell
935  factor = 0.5_dp; IF (nspins == 2) factor = 1.0_dp
936 
937 ! Allocate space for the projected density coefficients. On all ri_atoms.
938 ! Note: the sub-region where we project the density changes from excited atom to excited atom
939 ! => need different sets of RI coeffs
940  ALLOCATE (blk_size_ri(natom))
941  CALL get_particle_set(particle_set, qs_kind_set, nsgf=blk_size_ri, basis=basis_set_ri)
942  ALLOCATE (xas_atom_env%ri_dcoeff(natom, nspins, nex))
943  DO iex = 1, nex
944  DO ispin = 1, nspins
945  DO iat = 1, natom
946  NULLIFY (xas_atom_env%ri_dcoeff(iat, ispin, iex)%array)
947  IF ((.NOT. any(xas_atom_env%exat_neighbors(iex)%array == iat)) &
948  .AND. (.NOT. xas_atom_env%excited_atoms(iex) == iat)) cycle
949  ALLOCATE (xas_atom_env%ri_dcoeff(iat, ispin, iex)%array(blk_size_ri(iat)))
950  xas_atom_env%ri_dcoeff(iat, ispin, iex)%array = 0.0_dp
951  END DO
952  END DO
953  END DO
954 
955  output_unit = cp_logger_get_default_io_unit()
956  IF (output_unit > 0) THEN
957  WRITE (output_unit, fmt="(/,T7,A,/,T7,A)") &
958  "Excited atom, natoms in RI_REGION:", &
959  "---------------------------------"
960  END IF
961 
962  !We go atom by atom, first computing the integrals themselves that we put into a tensor, then we do
963  !the contraction with the density. We do that in the original dist, which is optimized for overlap
964 
965  ALLOCATE (blk_size_orb(natom))
966  CALL get_particle_set(particle_set, qs_kind_set, nsgf=blk_size_orb, basis=basis_set_orb)
967 
968  DO iex = 1, nex
969 
970  !get neighbors of current atom
971  exat = xas_atom_env%excited_atoms(iex)
972  nnb = 1 + SIZE(xas_atom_env%exat_neighbors(iex)%array)
973  ALLOCATE (neighbors(nnb))
974  neighbors(1) = exat
975  neighbors(2:nnb) = xas_atom_env%exat_neighbors(iex)%array(:)
976  CALL sort_unique(neighbors, unique)
977 
978  !link the atoms to their position in neighbors
979  ALLOCATE (idx_to_nb(natom))
980  idx_to_nb = 0
981  DO inb = 1, nnb
982  idx_to_nb(neighbors(inb)) = inb
983  END DO
984 
985  IF (output_unit > 0) THEN
986  WRITE (output_unit, fmt="(T7,I12,I21)") &
987  exat, nnb
988  END IF
989 
990  !Get the neighbor lists for the overlap integrals (abc), centers c on the current
991  !excited atom and its neighbors defined by RI_RADIUS
992  CALL build_xas_tdp_ovlp_nl(ab_list, basis_set_orb, basis_set_orb, qs_env)
993  CALL build_xas_tdp_3c_nl(ac_list, basis_set_orb, basis_set_ri, do_potential_id, &
994  qs_env, excited_atoms=neighbors)
995 
996  !Compute the 3-center overlap integrals
997  pot%potential_type = do_potential_id
998 
999  CALL create_pqx_tensor(pqx, ab_list, ac_list, dbcsr_dist, blk_size_orb, blk_size_orb, &
1000  blk_size_ri)
1001  CALL fill_pqx_tensor(pqx, ab_list, ac_list, basis_set_orb, basis_set_orb, basis_set_ri, &
1002  pot, qs_env)
1003 
1004  !Compute the RI inverse overlap matrix on the reduced RI basis that spans the excited
1005  !atom and its neighbors, ri_sinv is replicated over all procs
1006  CALL get_exat_ri_sinv(ri_sinv, ri_mats, neighbors, idx_to_nb, basis_set_ri, qs_env)
1007 
1008  !Do the actual contraction: coeff_y = sum_pq sum_x P_pq (phi_p phi_q xi_x) S_xy^-1
1009 
1010 !$OMP PARALLEL DEFAULT(NONE) &
1011 !$OMP SHARED(pqX,rho_ao,ri_sinv,xas_atom_env) &
1012 !$OMP SHARED(blk_size_ri,idx_to_nb,nspins,nnb,neighbors,iex,factor) &
1013 !$OMP PRIVATE(iter,ind,t_block,tensor_found,iatom,jatom,katom,inb,prefac,ispin) &
1014 !$OMP PRIVATE(pmat_block,pmat_found,pmat_blockt,pmat_foundt,work1,work2,jnb,ri_at) &
1015 !$OMP PRIVATE(sinv_block,sinv_found,sinv_blockt,sinv_foundt)
1016  CALL dbt_iterator_start(iter, pqx)
1017  DO WHILE (dbt_iterator_blocks_left(iter))
1018  CALL dbt_iterator_next_block(iter, ind)
1019  CALL dbt_get_block(pqx, ind, t_block, tensor_found)
1020 
1021  iatom = ind(1)
1022  jatom = ind(2)
1023  katom = ind(3)
1024  inb = idx_to_nb(katom)
1025 
1026  !non-diagonal elements need to be counted twice
1027  prefac = 2.0_dp
1028  IF (iatom == jatom) prefac = 1.0_dp
1029 
1030  DO ispin = 1, nspins
1031 
1032  !rho_ao is symmetric, block can be in either location
1033  CALL dbcsr_get_block_p(rho_ao(ispin)%matrix, iatom, jatom, pmat_block, pmat_found)
1034  CALL dbcsr_get_block_p(rho_ao(ispin)%matrix, jatom, iatom, pmat_blockt, pmat_foundt)
1035  IF ((.NOT. pmat_found) .AND. (.NOT. pmat_foundt)) cycle
1036 
1037  ALLOCATE (work1(blk_size_ri(katom), 1))
1038  work1 = 0.0_dp
1039 
1040  !first contraction with the density matrix
1041  IF (pmat_found) THEN
1042  DO i = 1, blk_size_ri(katom)
1043  work1(i, 1) = prefac*sum(pmat_block(:, :)*t_block(:, :, i))
1044  END DO
1045  ELSE
1046  DO i = 1, blk_size_ri(katom)
1047  work1(i, 1) = prefac*sum(transpose(pmat_blockt(:, :))*t_block(:, :, i))
1048  END DO
1049  END IF
1050 
1051  !loop over neighbors
1052  DO jnb = 1, nnb
1053 
1054  ri_at = neighbors(jnb)
1055 
1056  !ri_sinv is a symmetric matrix => actual block is one of the two
1057  CALL dbcsr_get_block_p(ri_sinv, inb, jnb, sinv_block, sinv_found)
1058  CALL dbcsr_get_block_p(ri_sinv, jnb, inb, sinv_blockt, sinv_foundt)
1059  IF ((.NOT. sinv_found) .AND. (.NOT. sinv_foundt)) cycle
1060 
1061  !second contraction with the inverse RI overlap
1062  ALLOCATE (work2(SIZE(xas_atom_env%ri_dcoeff(ri_at, ispin, iex)%array)))
1063  work2 = 0.0_dp
1064 
1065  IF (sinv_found) THEN
1066  DO i = 1, blk_size_ri(katom)
1067  work2(:) = work2(:) + factor*work1(i, 1)*sinv_block(i, :)
1068  END DO
1069  ELSE
1070  DO i = 1, blk_size_ri(katom)
1071  work2(:) = work2(:) + factor*work1(i, 1)*sinv_blockt(:, i)
1072  END DO
1073  END IF
1074  DO i = 1, SIZE(work2)
1075 !$OMP ATOMIC
1076  xas_atom_env%ri_dcoeff(ri_at, ispin, iex)%array(i) = &
1077  xas_atom_env%ri_dcoeff(ri_at, ispin, iex)%array(i) + work2(i)
1078  END DO
1079 
1080  DEALLOCATE (work2)
1081  END DO !jnb
1082 
1083  DEALLOCATE (work1)
1084  END DO
1085 
1086  DEALLOCATE (t_block)
1087  END DO !iter
1088  CALL dbt_iterator_stop(iter)
1089 !$OMP END PARALLEL
1090 
1091  !clean-up
1092  CALL dbcsr_release(ri_sinv)
1093  CALL dbt_destroy(pqx)
1094  CALL release_neighbor_list_sets(ab_list)
1095  CALL release_neighbor_list_sets(ac_list)
1096  DEALLOCATE (neighbors, idx_to_nb)
1097 
1098  END DO !iex
1099 
1100  !making sure all procs have the same info
1101  DO iex = 1, nex
1102  DO ispin = 1, nspins
1103  DO iat = 1, natom
1104  IF ((.NOT. any(xas_atom_env%exat_neighbors(iex)%array == iat)) &
1105  .AND. (.NOT. xas_atom_env%excited_atoms(iex) == iat)) cycle
1106  CALL para_env%sum(xas_atom_env%ri_dcoeff(iat, ispin, iex)%array)
1107  END DO !iat
1108  END DO !ispin
1109  END DO !iex
1110 
1111 ! clean-up
1112  CALL dbcsr_deallocate_matrix_set(overlap)
1113  DEALLOCATE (basis_set_ri, basis_set_orb, all_ri_atoms)
1114 
1115  CALL timestop(handle)
1116 
1117  END SUBROUTINE calculate_density_coeffs
1118 
1119 ! **************************************************************************************************
1120 !> \brief Evaluates the density on a given atomic grid
1121 !> \param rho_set where the densities are stored
1122 !> \param ri_dcoeff the arrays containing the RI density coefficients of this atom, for each spin
1123 !> \param atom_kind the kind of the atom in question
1124 !> \param do_gga whether the gradient of the density should also be put on the grid
1125 !> \param batch_info how the so are distributed
1126 !> \param xas_atom_env ...
1127 !> \param qs_env ...
1128 !> \note The density is expressed as n = sum_d coeff_d*xi_d. Knowing the coordinate of each grid
1129 !> grid point, one can simply evaluate xi_d(r)
1130 ! **************************************************************************************************
1131  SUBROUTINE put_density_on_atomic_grid(rho_set, ri_dcoeff, atom_kind, do_gga, batch_info, &
1132  xas_atom_env, qs_env)
1133 
1134  TYPE(xc_rho_set_type), INTENT(INOUT) :: rho_set
1135  TYPE(cp_1d_r_p_type), DIMENSION(:) :: ri_dcoeff
1136  INTEGER, INTENT(IN) :: atom_kind
1137  LOGICAL, INTENT(IN) :: do_gga
1138  TYPE(batch_info_type) :: batch_info
1139  TYPE(xas_atom_env_type), POINTER :: xas_atom_env
1140  TYPE(qs_environment_type), POINTER :: qs_env
1141 
1142  CHARACTER(len=*), PARAMETER :: routinen = 'put_density_on_atomic_grid'
1143 
1144  INTEGER :: dir, handle, ipgf, iset, iso, iso_proc, &
1145  ispin, maxso, n, na, nr, nset, nsgfi, &
1146  nsoi, nspins, sgfi, starti
1147  INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf, nsgf_set
1148  INTEGER, DIMENSION(:, :), POINTER :: first_sgf
1149  REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: so
1150  REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: dso
1151  REAL(dp), DIMENSION(:, :), POINTER :: dgr1, dgr2, ga, gr, ri_sphi_so, zet
1152  REAL(dp), DIMENSION(:, :, :), POINTER :: dga1, dga2, rhoa, rhob
1153  TYPE(cp_1d_r_p_type), DIMENSION(:), POINTER :: ri_dcoeff_so
1154  TYPE(grid_atom_type), POINTER :: grid_atom
1155  TYPE(gto_basis_set_type), POINTER :: ri_basis
1156  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1157 
1158  NULLIFY (grid_atom, ri_basis, qs_kind_set, lmax, npgf, zet, nsgf_set, ri_sphi_so)
1159  NULLIFY (lmin, first_sgf, rhoa, rhob, ga, gr, dgr1, dgr2, dga1, dga2, ri_dcoeff_so)
1160 
1161  CALL timeset(routinen, handle)
1162 
1163 ! Strategy: it makes sense to evaluate the spherical orbital on the grid (because of symmetry)
1164 ! From there, one can directly contract into sgf using ri_sphi_so and then take the weight
1165 ! The spherical orbital were precomputed and split in a purely radial and a purely
1166 ! angular part. The full values on each grid point are obtain through gemm
1167 
1168 ! Generalities
1169  CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set)
1170  CALL get_qs_kind(qs_kind_set(atom_kind), basis_set=ri_basis, basis_type="RI_XAS")
1171  CALL get_gto_basis_set(ri_basis, lmax=lmax, npgf=npgf, zet=zet, nset=nset, nsgf_set=nsgf_set, &
1172  first_sgf=first_sgf, lmin=lmin, maxso=maxso)
1173 
1174 ! Get the grid and the info we need from it
1175  grid_atom => xas_atom_env%grid_atom_set(atom_kind)%grid_atom
1176  na = grid_atom%ng_sphere
1177  nr = grid_atom%nr
1178  n = na*nr
1179  nspins = xas_atom_env%nspins
1180  ri_sphi_so => xas_atom_env%ri_sphi_so(atom_kind)%array
1181 
1182 ! Point to the rho_set densities
1183  rhoa => rho_set%rhoa
1184  rhob => rho_set%rhob
1185  rhoa = 0.0_dp; rhob = 0.0_dp;
1186  IF (do_gga) THEN
1187  DO dir = 1, 3
1188  rho_set%drhoa(dir)%array = 0.0_dp
1189  rho_set%drhob(dir)%array = 0.0_dp
1190  END DO
1191  END IF
1192 
1193 ! Point to the precomputed SO
1194  ga => xas_atom_env%ga(atom_kind)%array
1195  gr => xas_atom_env%gr(atom_kind)%array
1196  IF (do_gga) THEN
1197  dga1 => xas_atom_env%dga1(atom_kind)%array
1198  dga2 => xas_atom_env%dga2(atom_kind)%array
1199  dgr1 => xas_atom_env%dgr1(atom_kind)%array
1200  dgr2 => xas_atom_env%dgr2(atom_kind)%array
1201  ELSE
1202  dga1 => xas_atom_env%dga1(atom_kind)%array
1203  dga2 => xas_atom_env%dga2(atom_kind)%array
1204  dgr1 => xas_atom_env%dgr1(atom_kind)%array
1205  dgr2 => xas_atom_env%dgr2(atom_kind)%array
1206  END IF
1207 
1208 ! Need to express the ri_dcoeffs in terms of so (and not sgf)
1209  ALLOCATE (ri_dcoeff_so(nspins))
1210  DO ispin = 1, nspins
1211  ALLOCATE (ri_dcoeff_so(ispin)%array(nset*maxso))
1212  ri_dcoeff_so(ispin)%array = 0.0_dp
1213 
1214  !for a given so, loop over sgf and sum
1215  DO iset = 1, nset
1216  sgfi = first_sgf(1, iset)
1217  nsoi = npgf(iset)*nsoset(lmax(iset))
1218  nsgfi = nsgf_set(iset)
1219 
1220  CALL dgemv('N', nsoi, nsgfi, 1.0_dp, ri_sphi_so(1:nsoi, sgfi:sgfi + nsgfi - 1), nsoi, &
1221  ri_dcoeff(ispin)%array(sgfi:sgfi + nsgfi - 1), 1, 0.0_dp, &
1222  ri_dcoeff_so(ispin)%array((iset - 1)*maxso + 1:(iset - 1)*maxso + nsoi), 1)
1223 
1224  END DO
1225  END DO
1226 
1227  !allocate space to store the spherical orbitals on the grid
1228  ALLOCATE (so(na, nr))
1229  IF (do_gga) ALLOCATE (dso(na, nr, 3))
1230 
1231 ! Loop over the spherical orbitals on this proc
1232  DO iso_proc = 1, batch_info%nso_proc(atom_kind)
1233  iset = batch_info%so_proc_info(atom_kind)%array(1, iso_proc)
1234  ipgf = batch_info%so_proc_info(atom_kind)%array(2, iso_proc)
1235  iso = batch_info%so_proc_info(atom_kind)%array(3, iso_proc)
1236  IF (iso < 0) cycle
1237 
1238  starti = (iset - 1)*maxso + (ipgf - 1)*nsoset(lmax(iset))
1239 
1240  !the spherical orbital on the grid
1241  CALL dgemm('N', 'T', na, nr, 1, 1.0_dp, ga(:, iso_proc:iso_proc), na, &
1242  gr(:, iso_proc:iso_proc), nr, 0.0_dp, so(:, :), na)
1243 
1244  !the gradient on the grid
1245  IF (do_gga) THEN
1246 
1247  DO dir = 1, 3
1248  CALL dgemm('N', 'T', na, nr, 1, 1.0_dp, dga1(:, iso_proc:iso_proc, dir), na, &
1249  dgr1(:, iso_proc:iso_proc), nr, 0.0_dp, dso(:, :, dir), na)
1250  CALL dgemm('N', 'T', na, nr, 1, 1.0_dp, dga2(:, iso_proc:iso_proc, dir), na, &
1251  dgr2(:, iso_proc:iso_proc), nr, 1.0_dp, dso(:, :, dir), na)
1252  END DO
1253  END IF
1254 
1255  !put the so on the grid with the approriate coefficients and sum
1256  CALL daxpy(n, ri_dcoeff_so(1)%array(starti + iso), so, 1, rhoa(:, :, 1), 1)
1257 
1258  IF (nspins == 2) THEN
1259  CALL daxpy(n, ri_dcoeff_so(2)%array(starti + iso), so, 1, rhob(:, :, 1), 1)
1260  END IF
1261 
1262  IF (do_gga) THEN
1263 
1264  !put the gradient of the so on the grid with correspond RI coeff
1265  DO dir = 1, 3
1266  CALL daxpy(n, ri_dcoeff_so(1)%array(starti + iso), dso(:, :, dir), &
1267  1, rho_set%drhoa(dir)%array(:, :, 1), 1)
1268 
1269  IF (nspins == 2) THEN
1270  CALL daxpy(n, ri_dcoeff_so(2)%array(starti + iso), dso(:, :, dir), &
1271  1, rho_set%drhob(dir)%array(:, :, 1), 1)
1272  END IF
1273  END DO !dir
1274  END IF !do_gga
1275 
1276  END DO
1277 
1278 ! Treat spin restricted case (=> copy alpha into beta)
1279  IF (nspins == 1) THEN
1280  CALL dcopy(n, rhoa(:, :, 1), 1, rhob(:, :, 1), 1)
1281 
1282  IF (do_gga) THEN
1283  DO dir = 1, 3
1284  CALL dcopy(n, rho_set%drhoa(dir)%array(:, :, 1), 1, rho_set%drhob(dir)%array(:, :, 1), 1)
1285  END DO
1286  END IF
1287  END IF
1288 
1289 ! Note: sum over procs is done outside
1290 
1291 ! clean-up
1292  DO ispin = 1, nspins
1293  DEALLOCATE (ri_dcoeff_so(ispin)%array)
1294  END DO
1295  DEALLOCATE (ri_dcoeff_so)
1296 
1297  CALL timestop(handle)
1298 
1299  END SUBROUTINE put_density_on_atomic_grid
1300 
1301 ! **************************************************************************************************
1302 !> \brief Adds the density of a given source atom with source kind (with ri_dcoeff) on the atomic
1303 !> grid belonging to another target atom of target kind. The evaluations of the basis
1304 !> function first requires the evaluation of the x,y,z coordinates on each grid point of
1305 !> target atom wrt to the position of source atom
1306 !> \param rho_set where the densities are stored
1307 !> \param ri_dcoeff the arrays containing the RI density coefficient of source_iat, for each spin
1308 !> \param source_iat the index of the source atom
1309 !> \param source_ikind the kind of the source atom
1310 !> \param target_iat the index of the target atom
1311 !> \param target_ikind the kind of the target atom
1312 !> \param sr starting r index for the local grid
1313 !> \param er ending r index for the local grid
1314 !> \param do_gga whether the gradient of the density is needed
1315 !> \param xas_atom_env ...
1316 !> \param qs_env ...
1317 ! **************************************************************************************************
1318  SUBROUTINE put_density_on_other_grid(rho_set, ri_dcoeff, source_iat, source_ikind, target_iat, &
1319  target_ikind, sr, er, do_gga, xas_atom_env, qs_env)
1320 
1321  TYPE(xc_rho_set_type), INTENT(INOUT) :: rho_set
1322  TYPE(cp_1d_r_p_type), DIMENSION(:) :: ri_dcoeff
1323  INTEGER, INTENT(IN) :: source_iat, source_ikind, target_iat, &
1324  target_ikind, sr, er
1325  LOGICAL, INTENT(IN) :: do_gga
1326  TYPE(xas_atom_env_type), POINTER :: xas_atom_env
1327  TYPE(qs_environment_type), POINTER :: qs_env
1328 
1329  CHARACTER(len=*), PARAMETER :: routinen = 'put_density_on_other_grid'
1330 
1331  INTEGER :: dir, handle, ia, ico, ipgf, ir, iset, &
1332  isgf, lx, ly, lz, n, na, nr, nset, &
1333  nspins, sgfi, start
1334  INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf, nsgf_set
1335  INTEGER, DIMENSION(:, :), POINTER :: first_sgf
1336  REAL(dp) :: rmom
1337  REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: sgf
1338  REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: co, dsgf, pos
1339  REAL(dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: dco
1340  REAL(dp), DIMENSION(3) :: rs, rst, rt
1341  REAL(dp), DIMENSION(:, :), POINTER :: ri_sphi, zet
1342  REAL(dp), DIMENSION(:, :, :), POINTER :: rhoa, rhob
1343  TYPE(cell_type), POINTER :: cell
1344  TYPE(grid_atom_type), POINTER :: grid_atom
1345  TYPE(gto_basis_set_type), POINTER :: ri_basis
1346  TYPE(harmonics_atom_type), POINTER :: harmonics
1347  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1348  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1349 
1350  NULLIFY (qs_kind_set, ri_basis, lmax, npgf, nsgf_set, lmin, zet, first_sgf, grid_atom)
1351  NULLIFY (harmonics, rhoa, rhob, particle_set, cell, ri_sphi)
1352 
1353  !Same logic as the put_density_on_own_grid routine. Loop over orbitals, put them on the grid,
1354  !contract into sgf and daxpy with coeff. Notable difference: use cartesian orbitals instead of
1355  !spherical, since the center of the grid is not the origin and thus, spherical symmetry can't
1356  !be exploited so well
1357 
1358  CALL timeset(routinen, handle)
1359 
1360  !Generalities
1361  CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, particle_set=particle_set, cell=cell)
1362  !want basis of the source atom
1363  CALL get_qs_kind(qs_kind_set(source_ikind), basis_set=ri_basis, basis_type="RI_XAS")
1364  CALL get_gto_basis_set(ri_basis, lmax=lmax, npgf=npgf, zet=zet, nset=nset, nsgf_set=nsgf_set, &
1365  first_sgf=first_sgf, lmin=lmin, sphi=ri_sphi)
1366 
1367  ! Want the grid and harmonics of the target atom
1368  grid_atom => xas_atom_env%grid_atom_set(target_ikind)%grid_atom
1369  harmonics => xas_atom_env%harmonics_atom_set(target_ikind)%harmonics_atom
1370  na = grid_atom%ng_sphere
1371  nr = er - sr + 1
1372  n = na*nr
1373  nspins = xas_atom_env%nspins
1374 
1375  ! Point to the rho_set densities
1376  rhoa => rho_set%rhoa
1377  rhob => rho_set%rhob
1378 
1379  ! Need the source-target position vector
1380  rs = pbc(particle_set(source_iat)%r, cell)
1381  rt = pbc(particle_set(target_iat)%r, cell)
1382  rst = pbc(rs, rt, cell)
1383 
1384  ! Precompute the positions on the target grid
1385  ALLOCATE (pos(na, sr:er, 4))
1386 !$OMP PARALLEL DO COLLAPSE(2) SCHEDULE(STATIC) DEFAULT(NONE), &
1387 !$OMP SHARED(na,sr,er,pos,harmonics,grid_atom,rst), &
1388 !$OMP PRIVATE(ia,ir)
1389  DO ir = sr, er
1390  DO ia = 1, na
1391  pos(ia, ir, 1:3) = harmonics%a(:, ia)*grid_atom%rad(ir) + rst
1392  pos(ia, ir, 4) = pos(ia, ir, 1)**2 + pos(ia, ir, 2)**2 + pos(ia, ir, 3)**2
1393  END DO
1394  END DO
1395 !$OMP END PARALLEL DO
1396 
1397  ! Loop over the cartesian gaussian functions and evaluate them
1398  DO iset = 1, nset
1399 
1400  !allocate space to store the cartesian orbtial on the grid
1401  ALLOCATE (co(na, sr:er, npgf(iset)*ncoset(lmax(iset))))
1402  IF (do_gga) ALLOCATE (dco(na, sr:er, 3, npgf(iset)*ncoset(lmax(iset))))
1403 
1404 !$OMP PARALLEL DEFAULT(NONE), &
1405 !$OMP SHARED(co,npgf,ncoset,lmax,lmin,indco,pos,zet,iset,na,sr,er,do_gga,dco), &
1406 !$OMP PRIVATE(ipgf,start,ico,lx,ly,lz,ia,ir,rmom)
1407 
1408 !$OMP DO COLLAPSE(2) SCHEDULE(STATIC)
1409  DO ir = sr, er
1410  DO ia = 1, na
1411  co(ia, ir, :) = 0.0_dp
1412  IF (do_gga) THEN
1413  dco(ia, ir, :, :) = 0.0_dp
1414  END IF
1415  END DO
1416  END DO
1417 !$OMP END DO NOWAIT
1418 
1419  DO ipgf = 1, npgf(iset)
1420  start = (ipgf - 1)*ncoset(lmax(iset))
1421 
1422  !loop over the cartesian orbitals
1423  DO ico = ncoset(lmin(iset) - 1) + 1, ncoset(lmax(iset))
1424  lx = indco(1, ico)
1425  ly = indco(2, ico)
1426  lz = indco(3, ico)
1427 
1428  ! compute g = x**lx * y**ly * z**lz * exp(-zet * r**2)
1429 !$OMP DO COLLAPSE(2) SCHEDULE(STATIC)
1430  DO ir = sr, er
1431  DO ia = 1, na
1432  rmom = exp(-zet(ipgf, iset)*pos(ia, ir, 4))
1433  IF (lx /= 0) rmom = rmom*pos(ia, ir, 1)**lx
1434  IF (ly /= 0) rmom = rmom*pos(ia, ir, 2)**ly
1435  IF (lz /= 0) rmom = rmom*pos(ia, ir, 3)**lz
1436  co(ia, ir, start + ico) = rmom
1437  !co(ia, ir, start + ico) = pos(ia, ir, 1)**lx*pos(ia, ir, 2)**ly*pos(ia, ir, 3)**lz &
1438  ! *EXP(-zet(ipgf, iset)*pos(ia, ir, 4))
1439  END DO
1440  END DO
1441 !$OMP END DO NOWAIT
1442 
1443  IF (do_gga) THEN
1444  !the gradient: dg_x = lx*x**(lx-1) * y**ly * z**lz * exp(-zet * r**2)
1445  ! -2*zet* x**(lx+1) * y**ly * z**lz * exp(-zet * r**2)
1446  ! = (lx*x**(lx-1) - 2*zet*x**(lx+1)) * y**ly * z**lz * exp(-zet * r**2)
1447 
1448  !x direction, special case if lx == 0
1449  IF (lx == 0) THEN
1450 !$OMP DO COLLAPSE(2) SCHEDULE(STATIC)
1451  DO ir = sr, er
1452  DO ia = 1, na
1453  rmom = -2.0_dp*pos(ia, ir, 1)*zet(ipgf, iset)*exp(-zet(ipgf, iset)*pos(ia, ir, 4))
1454  IF (ly /= 0) rmom = rmom*pos(ia, ir, 2)**ly
1455  IF (lz /= 0) rmom = rmom*pos(ia, ir, 3)**lz
1456  dco(ia, ir, 1, start + ico) = rmom
1457 ! dco(ia, ir, 1, start + ico) = -2.0_dp*pos(ia, ir, 1)*zet(ipgf, iset) &
1458 ! *pos(ia, ir, 2)**ly*pos(ia, ir, 3)**lz &
1459 ! *EXP(-zet(ipgf, iset)*pos(ia, ir, 4))
1460  END DO
1461  END DO
1462 !$OMP END DO NOWAIT
1463  ELSE
1464 !$OMP DO COLLAPSE(2) SCHEDULE(STATIC)
1465  DO ir = sr, er
1466  DO ia = 1, na
1467  IF (lx /= 1) THEN
1468  rmom = (lx*pos(ia, ir, 1)**(lx - 1) - 2.0_dp*pos(ia, ir, 1)**(lx + 1)* &
1469  zet(ipgf, iset))*exp(-zet(ipgf, iset)*pos(ia, ir, 4))
1470  ELSE
1471  rmom = (1.0_dp - 2.0_dp*pos(ia, ir, 1)**2*zet(ipgf, iset))* &
1472  exp(-zet(ipgf, iset)*pos(ia, ir, 4))
1473  END IF
1474  IF (ly /= 0) rmom = rmom*pos(ia, ir, 2)**ly
1475  IF (lz /= 0) rmom = rmom*pos(ia, ir, 3)**lz
1476  dco(ia, ir, 1, start + ico) = rmom
1477 ! dco(ia, ir, 1, start + ico) = (lx*pos(ia, ir, 1)**(lx - 1) &
1478 ! - 2.0_dp*pos(ia, ir, 1)**(lx + 1)*zet(ipgf, iset)) &
1479 ! *pos(ia, ir, 2)**ly*pos(ia, ir, 3)**lz &
1480 ! *EXP(-zet(ipgf, iset)*pos(ia, ir, 4))
1481  END DO
1482  END DO
1483 !$OMP END DO NOWAIT
1484  END IF !lx == 0
1485 
1486  !y direction, special case if ly == 0
1487  IF (ly == 0) THEN
1488 !$OMP DO COLLAPSE(2) SCHEDULE(STATIC)
1489  DO ir = sr, er
1490  DO ia = 1, na
1491  rmom = -2.0_dp*pos(ia, ir, 2)*zet(ipgf, iset)*exp(-zet(ipgf, iset)*pos(ia, ir, 4))
1492  IF (lx /= 0) rmom = rmom*pos(ia, ir, 1)**lx
1493  IF (lz /= 0) rmom = rmom*pos(ia, ir, 3)**lz
1494  dco(ia, ir, 2, start + ico) = rmom
1495 ! dco(ia, ir, 2, start + ico) = -2.0_dp*pos(ia, ir, 2)*zet(ipgf, iset) &
1496 ! *pos(ia, ir, 1)**lx*pos(ia, ir, 3)**lz &
1497 ! *EXP(-zet(ipgf, iset)*pos(ia, ir, 4))
1498  END DO
1499  END DO
1500 !$OMP END DO NOWAIT
1501  ELSE
1502 !$OMP DO COLLAPSE(2) SCHEDULE(STATIC)
1503  DO ir = sr, er
1504  DO ia = 1, na
1505  IF (ly /= 1) THEN
1506  rmom = (ly*pos(ia, ir, 2)**(ly - 1) - 2.0_dp*pos(ia, ir, 2)**(ly + 1)*zet(ipgf, iset)) &
1507  *exp(-zet(ipgf, iset)*pos(ia, ir, 4))
1508  ELSE
1509  rmom = (1.0_dp - 2.0_dp*pos(ia, ir, 2)**2*zet(ipgf, iset)) &
1510  *exp(-zet(ipgf, iset)*pos(ia, ir, 4))
1511  END IF
1512  IF (lx /= 0) rmom = rmom*pos(ia, ir, 1)**lx
1513  IF (lz /= 0) rmom = rmom*pos(ia, ir, 3)**lz
1514  dco(ia, ir, 2, start + ico) = rmom
1515 ! dco(ia, ir, 2, start + ico) = (ly*pos(ia, ir, 2)**(ly - 1) &
1516 ! - 2.0_dp*pos(ia, ir, 2)**(ly + 1)*zet(ipgf, iset)) &
1517 ! *pos(ia, ir, 1)**lx*pos(ia, ir, 3)**lz &
1518 ! *EXP(-zet(ipgf, iset)*pos(ia, ir, 4))
1519  END DO
1520  END DO
1521 !$OMP END DO NOWAIT
1522  END IF !ly == 0
1523 
1524  !z direction, special case if lz == 0
1525  IF (lz == 0) THEN
1526 !$OMP DO COLLAPSE(2) SCHEDULE(STATIC)
1527  DO ir = sr, er
1528  DO ia = 1, na
1529  rmom = -2.0_dp*pos(ia, ir, 3)*zet(ipgf, iset)*exp(-zet(ipgf, iset)*pos(ia, ir, 4))
1530  IF (lx /= 0) rmom = rmom*pos(ia, ir, 1)**lx
1531  IF (ly /= 0) rmom = rmom*pos(ia, ir, 2)**ly
1532  dco(ia, ir, 3, start + ico) = rmom
1533 ! dco(ia, ir, 3, start + ico) = -2.0_dp*pos(ia, ir, 3)*zet(ipgf, iset) &
1534 ! *pos(ia, ir, 1)**lx*pos(ia, ir, 2)**ly &
1535 ! *EXP(-zet(ipgf, iset)*pos(ia, ir, 4))
1536  END DO
1537  END DO
1538 !$OMP END DO NOWAIT
1539  ELSE
1540 !$OMP DO COLLAPSE(2) SCHEDULE(STATIC)
1541  DO ir = sr, er
1542  DO ia = 1, na
1543  IF (lz /= 1) THEN
1544  rmom = (lz*pos(ia, ir, 3)**(lz - 1) - 2.0_dp*pos(ia, ir, 3)**(lz + 1)* &
1545  zet(ipgf, iset))*exp(-zet(ipgf, iset)*pos(ia, ir, 4))
1546  ELSE
1547  rmom = (1.0_dp - 2.0_dp*pos(ia, ir, 3)**2*zet(ipgf, iset))* &
1548  exp(-zet(ipgf, iset)*pos(ia, ir, 4))
1549  END IF
1550  IF (lx /= 0) rmom = rmom*pos(ia, ir, 1)**lx
1551  IF (ly /= 0) rmom = rmom*pos(ia, ir, 2)**ly
1552  dco(ia, ir, 3, start + ico) = rmom
1553 ! dco(ia, ir, 3, start + ico) = (lz*pos(ia, ir, 3)**(lz - 1) &
1554 ! - 2.0_dp*pos(ia, ir, 3)**(lz + 1)*zet(ipgf, iset)) &
1555 ! *pos(ia, ir, 1)**lx*pos(ia, ir, 2)**ly &
1556 ! *EXP(-zet(ipgf, iset)*pos(ia, ir, 4))
1557  END DO
1558  END DO
1559 !$OMP END DO NOWAIT
1560  END IF !lz == 0
1561 
1562  END IF !gga
1563 
1564  END DO !ico
1565  END DO !ipgf
1566 
1567 !$OMP END PARALLEL
1568 
1569  !contract the co into sgf
1570  ALLOCATE (sgf(na, sr:er))
1571  IF (do_gga) ALLOCATE (dsgf(na, sr:er, 3))
1572  sgfi = first_sgf(1, iset) - 1
1573 
1574  DO isgf = 1, nsgf_set(iset)
1575  sgf = 0.0_dp
1576  IF (do_gga) dsgf = 0.0_dp
1577 
1578  DO ipgf = 1, npgf(iset)
1579  start = (ipgf - 1)*ncoset(lmax(iset))
1580  DO ico = ncoset(lmin(iset) - 1) + 1, ncoset(lmax(iset))
1581  CALL daxpy(n, ri_sphi(start + ico, sgfi + isgf), co(:, sr:er, start + ico), 1, sgf(:, sr:er), 1)
1582  END DO !ico
1583  END DO !ipgf
1584 
1585  !add the density to the grid
1586  CALL daxpy(n, ri_dcoeff(1)%array(sgfi + isgf), sgf(:, sr:er), 1, rhoa(:, sr:er, 1), 1)
1587 
1588  IF (nspins == 2) THEN
1589  CALL daxpy(n, ri_dcoeff(2)%array(sgfi + isgf), sgf(:, sr:er), 1, rhob(:, sr:er, 1), 1)
1590  END IF
1591 
1592  !deal with the gradient
1593  IF (do_gga) THEN
1594 
1595  DO ipgf = 1, npgf(iset)
1596  start = (ipgf - 1)*ncoset(lmax(iset))
1597  DO ico = ncoset(lmin(iset) - 1) + 1, ncoset(lmax(iset))
1598  DO dir = 1, 3
1599  CALL daxpy(n, ri_sphi(start + ico, sgfi + isgf), dco(:, sr:er, dir, start + ico), &
1600  1, dsgf(:, sr:er, dir), 1)
1601  END DO
1602  END DO !ico
1603  END DO !ipgf
1604 
1605  DO dir = 1, 3
1606  CALL daxpy(n, ri_dcoeff(1)%array(sgfi + isgf), dsgf(:, sr:er, dir), 1, &
1607  rho_set%drhoa(dir)%array(:, sr:er, 1), 1)
1608 
1609  IF (nspins == 2) THEN
1610  CALL daxpy(n, ri_dcoeff(2)%array(sgfi + isgf), dsgf(:, sr:er, dir), 1, &
1611  rho_set%drhob(dir)%array(:, sr:er, 1), 1)
1612  END IF
1613  END DO
1614  END IF !do_gga
1615 
1616  END DO !isgf
1617 
1618  DEALLOCATE (co, sgf)
1619  IF (do_gga) DEALLOCATE (dco, dsgf)
1620  END DO !iset
1621 
1622  !Treat spin-restricted case (copy alpha into beta)
1623  IF (nspins == 1) THEN
1624  CALL dcopy(n, rhoa(:, sr:er, 1), 1, rhob(:, sr:er, 1), 1)
1625 
1626  IF (do_gga) THEN
1627  DO dir = 1, 3
1628  CALL dcopy(n, rho_set%drhoa(dir)%array(:, sr:er, 1), 1, rho_set%drhob(dir)%array(:, sr:er, 1), 1)
1629  END DO
1630  END IF
1631  END IF
1632 
1633  CALL timestop(handle)
1634 
1635  END SUBROUTINE put_density_on_other_grid
1636 
1637 ! **************************************************************************************************
1638 !> \brief Computes the norm of the density gradient on the atomic grid
1639 !> \param rho_set ...
1640 !> \param atom_kind ...
1641 !> \param xas_atom_env ...
1642 !> \note GGA is assumed
1643 ! **************************************************************************************************
1644  SUBROUTINE compute_norm_drho(rho_set, atom_kind, xas_atom_env)
1645 
1646  TYPE(xc_rho_set_type), INTENT(INOUT) :: rho_set
1647  INTEGER, INTENT(IN) :: atom_kind
1648  TYPE(xas_atom_env_type), POINTER :: xas_atom_env
1649 
1650  INTEGER :: dir, ia, ir, n, na, nr, nspins
1651 
1652  na = xas_atom_env%grid_atom_set(atom_kind)%grid_atom%ng_sphere
1653  nr = xas_atom_env%grid_atom_set(atom_kind)%grid_atom%nr
1654  n = na*nr
1655  nspins = xas_atom_env%nspins
1656 
1657  rho_set%norm_drhoa = 0.0_dp
1658  rho_set%norm_drhob = 0.0_dp
1659  rho_set%norm_drho = 0.0_dp
1660 
1661  DO dir = 1, 3
1662  DO ir = 1, nr
1663  DO ia = 1, na
1664  rho_set%norm_drhoa(ia, ir, 1) = rho_set%norm_drhoa(ia, ir, 1) &
1665  + rho_set%drhoa(dir)%array(ia, ir, 1)**2
1666  END DO !ia
1667  END DO !ir
1668  END DO !dir
1669  rho_set%norm_drhoa = sqrt(rho_set%norm_drhoa)
1670 
1671  IF (nspins == 1) THEN
1672  !spin-restricted
1673  CALL dcopy(n, rho_set%norm_drhoa(:, :, 1), 1, rho_set%norm_drhob(:, :, 1), 1)
1674  ELSE
1675  DO dir = 1, 3
1676  DO ir = 1, nr
1677  DO ia = 1, na
1678  rho_set%norm_drhob(ia, ir, 1) = rho_set%norm_drhob(ia, ir, 1) &
1679  + rho_set%drhob(dir)%array(ia, ir, 1)**2
1680  END DO
1681  END DO
1682  END DO
1683  rho_set%norm_drhob = sqrt(rho_set%norm_drhob)
1684  END IF
1685 
1686  DO dir = 1, 3
1687  DO ir = 1, nr
1688  DO ia = 1, na
1689  rho_set%norm_drho(ia, ir, 1) = rho_set%norm_drho(ia, ir, 1) + &
1690  (rho_set%drhoa(dir)%array(ia, ir, 1) + &
1691  rho_set%drhob(dir)%array(ia, ir, 1))**2
1692  END DO
1693  END DO
1694  END DO
1695  rho_set%norm_drho = sqrt(rho_set%norm_drho)
1696 
1697  END SUBROUTINE compute_norm_drho
1698 
1699 ! **************************************************************************************************
1700 !> \brief Precomputes the spherical orbitals of the RI basis on the excited atom grids
1701 !> \param do_gga whether the gradient needs to be computed for GGA or not
1702 !> \param batch_info the parallelization info to complete with so distribution info
1703 !> \param xas_atom_env ...
1704 !> \param qs_env ...
1705 !> \note the functions are split in a purely angular part of size na and a purely radial part of
1706 !> size nr. The full function on the grid can simply be obtained with dgemm and we save space
1707 ! **************************************************************************************************
1708  SUBROUTINE precompute_so_dso(do_gga, batch_info, xas_atom_env, qs_env)
1709 
1710  LOGICAL, INTENT(IN) :: do_gga
1711  TYPE(batch_info_type) :: batch_info
1712  TYPE(xas_atom_env_type), POINTER :: xas_atom_env
1713  TYPE(qs_environment_type), POINTER :: qs_env
1714 
1715  CHARACTER(len=*), PARAMETER :: routinen = 'precompute_so_dso'
1716 
1717  INTEGER :: bo(2), dir, handle, ikind, ipgf, iset, &
1718  iso, iso_proc, l, maxso, n, na, nkind, &
1719  nr, nset, nso_proc, nsotot, starti
1720  INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf, nsgf_set
1721  INTEGER, DIMENSION(:, :), POINTER :: so_proc_info
1722  REAL(dp), ALLOCATABLE, DIMENSION(:) :: rexp
1723  REAL(dp), DIMENSION(:, :), POINTER :: dgr1, dgr2, ga, gr, slm, zet
1724  REAL(dp), DIMENSION(:, :, :), POINTER :: dga1, dga2, dslm_dxyz
1725  TYPE(grid_atom_type), POINTER :: grid_atom
1726  TYPE(gto_basis_set_type), POINTER :: ri_basis
1727  TYPE(harmonics_atom_type), POINTER :: harmonics
1728  TYPE(mp_para_env_type), POINTER :: para_env
1729  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1730 
1731  NULLIFY (qs_kind_set, harmonics, grid_atom, slm, dslm_dxyz, ri_basis, lmax, lmin, npgf)
1732  NULLIFY (nsgf_set, zet, para_env, so_proc_info)
1733 
1734  CALL timeset(routinen, handle)
1735 
1736  CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, para_env=para_env)
1737  nkind = SIZE(qs_kind_set)
1738 
1739  ALLOCATE (batch_info%so_proc_info(nkind))
1740  ALLOCATE (batch_info%nso_proc(nkind))
1741  ALLOCATE (batch_info%so_bo(2, nkind))
1742 
1743  DO ikind = 1, nkind
1744 
1745  NULLIFY (xas_atom_env%ga(ikind)%array)
1746  NULLIFY (xas_atom_env%gr(ikind)%array)
1747  NULLIFY (xas_atom_env%dga1(ikind)%array)
1748  NULLIFY (xas_atom_env%dga2(ikind)%array)
1749  NULLIFY (xas_atom_env%dgr1(ikind)%array)
1750  NULLIFY (xas_atom_env%dgr2(ikind)%array)
1751 
1752  NULLIFY (batch_info%so_proc_info(ikind)%array)
1753 
1754  IF (.NOT. any(xas_atom_env%excited_kinds == ikind)) cycle
1755 
1756  !grid info
1757  harmonics => xas_atom_env%harmonics_atom_set(ikind)%harmonics_atom
1758  grid_atom => xas_atom_env%grid_atom_set(ikind)%grid_atom
1759 
1760  na = grid_atom%ng_sphere
1761  nr = grid_atom%nr
1762  n = na*nr
1763 
1764  slm => harmonics%slm
1765  dslm_dxyz => harmonics%dslm_dxyz
1766 
1767  !basis info
1768  CALL get_qs_kind(qs_kind_set(ikind), basis_set=ri_basis, basis_type="RI_XAS")
1769  CALL get_gto_basis_set(ri_basis, lmax=lmax, npgf=npgf, zet=zet, nset=nset, &
1770  nsgf_set=nsgf_set, lmin=lmin, maxso=maxso)
1771  nsotot = maxso*nset
1772 
1773  !we split all so among the processors of the batch
1774  bo = get_limit(nsotot, batch_info%batch_size, batch_info%ipe)
1775  nso_proc = bo(2) - bo(1) + 1
1776  batch_info%so_bo(:, ikind) = bo
1777  batch_info%nso_proc(ikind) = nso_proc
1778 
1779  !store info about the so's set, pgf and index
1780  ALLOCATE (batch_info%so_proc_info(ikind)%array(3, nso_proc))
1781  so_proc_info => batch_info%so_proc_info(ikind)%array
1782  so_proc_info = -1 !default is -1 => set so value to zero
1783  DO iset = 1, nset
1784  DO ipgf = 1, npgf(iset)
1785  starti = (iset - 1)*maxso + (ipgf - 1)*nsoset(lmax(iset))
1786  DO iso = nsoset(lmin(iset) - 1) + 1, nsoset(lmax(iset))
1787 
1788  !only consider so that are on this proc
1789  IF (starti + iso < bo(1) .OR. starti + iso > bo(2)) cycle
1790  iso_proc = starti + iso - bo(1) + 1
1791  so_proc_info(1, iso_proc) = iset
1792  so_proc_info(2, iso_proc) = ipgf
1793  so_proc_info(3, iso_proc) = iso
1794 
1795  END DO
1796  END DO
1797  END DO
1798 
1799  !Put the gaussians and their gradient as purely angular or radial arrays
1800  ALLOCATE (xas_atom_env%ga(ikind)%array(na, nso_proc))
1801  ALLOCATE (xas_atom_env%gr(ikind)%array(nr, nso_proc))
1802  xas_atom_env%ga(ikind)%array = 0.0_dp; xas_atom_env%gr(ikind)%array = 0.0_dp
1803  IF (do_gga) THEN
1804  ALLOCATE (xas_atom_env%dga1(ikind)%array(na, nso_proc, 3))
1805  ALLOCATE (xas_atom_env%dgr1(ikind)%array(nr, nso_proc))
1806  ALLOCATE (xas_atom_env%dga2(ikind)%array(na, nso_proc, 3))
1807  ALLOCATE (xas_atom_env%dgr2(ikind)%array(nr, nso_proc))
1808  xas_atom_env%dga1(ikind)%array = 0.0_dp; xas_atom_env%dgr1(ikind)%array = 0.0_dp
1809  xas_atom_env%dga2(ikind)%array = 0.0_dp; xas_atom_env%dgr2(ikind)%array = 0.0_dp
1810  END IF
1811 
1812  ga => xas_atom_env%ga(ikind)%array
1813  gr => xas_atom_env%gr(ikind)%array
1814  dga1 => xas_atom_env%dga1(ikind)%array
1815  dga2 => xas_atom_env%dga2(ikind)%array
1816  dgr1 => xas_atom_env%dgr1(ikind)%array
1817  dgr2 => xas_atom_env%dgr2(ikind)%array
1818 
1819  ALLOCATE (rexp(nr))
1820 
1821  DO iso_proc = 1, nso_proc
1822  iset = so_proc_info(1, iso_proc)
1823  ipgf = so_proc_info(2, iso_proc)
1824  iso = so_proc_info(3, iso_proc)
1825  IF (iso < 0) cycle
1826 
1827  l = indso(1, iso)
1828 
1829  !The gaussian is g = r^l * Ylm * exp(-a*r^2)
1830 
1831  !radial part of the gaussian
1832  rexp(1:nr) = exp(-zet(ipgf, iset)*grid_atom%rad2(1:nr))
1833  gr(1:nr, iso_proc) = grid_atom%rad(1:nr)**l*rexp(1:nr)
1834 
1835  !angular part of the gaussian
1836  ga(1:na, iso_proc) = slm(1:na, iso)
1837 
1838  !For the gradient, devide in 2 parts: dg/dx = d/dx(r^l * Ylm) * exp(-a*r^2)
1839  ! + r^l * Ylm * d/dx(exp(-a*r^2))
1840  !Note: we make this choice of separation because of cartesian coordinates, where
1841  ! g = x^lx * y^ly * z^lz * exp(-a*r^2) and r^(l-1)*dslm_dxyz = d/dx(r^l * Ylm)
1842 
1843  IF (do_gga) THEN
1844  !radial part of the gradient => same in all three direction
1845  dgr1(1:nr, iso_proc) = grid_atom%rad(1:nr)**(l - 1)*rexp(1:nr)
1846  dgr2(1:nr, iso_proc) = -2.0_dp*zet(ipgf, iset)*grid_atom%rad(1:nr)**(l + 1)*rexp(1:nr)
1847 
1848  !angular part of the gradient
1849  DO dir = 1, 3
1850  dga1(1:na, iso_proc, dir) = dslm_dxyz(dir, 1:na, iso)
1851  dga2(1:na, iso_proc, dir) = harmonics%a(dir, 1:na)*slm(1:na, iso)
1852  END DO
1853  END IF
1854 
1855  END DO !iso_proc
1856 
1857  DEALLOCATE (rexp)
1858  END DO !ikind
1859 
1860  CALL timestop(handle)
1861 
1862  END SUBROUTINE precompute_so_dso
1863 
1864 ! **************************************************************************************************
1865 !> \brief Integrate the xc kernel as a function of r on the atomic grids for the RI_XAS basis
1866 !> \param int_fxc the global array containing the (P|fxc|Q) integrals, for all spin configurations
1867 !> \param xas_atom_env ...
1868 !> \param xas_tdp_control ...
1869 !> \param qs_env ...
1870 !> \note Note that if closed-shell, alpha-alpha term and beta-beta terms are the same
1871 !> Store the (P|fxc|Q) integrals on the processor they were computed on
1872 !> int_fxc(1)%matrix is alpha-alpha, 2: alpha-beta, 3: beta-beta
1873 ! **************************************************************************************************
1874  SUBROUTINE integrate_fxc_atoms(int_fxc, xas_atom_env, xas_tdp_control, qs_env)
1875 
1876  TYPE(cp_2d_r_p_type), DIMENSION(:, :), POINTER :: int_fxc
1877  TYPE(xas_atom_env_type), POINTER :: xas_atom_env
1878  TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
1879  TYPE(qs_environment_type), POINTER :: qs_env
1880 
1881  CHARACTER(len=*), PARAMETER :: routinen = 'integrate_fxc_atoms'
1882 
1883  INTEGER :: batch_size, dir, er, ex_bo(2), handle, i, iatom, ibatch, iex, ikind, inb, ipe, &
1884  mepos, na, natom, nb, nb_bo(2), nbatch, nbk, nex_atom, nr, num_pe, sr
1885  INTEGER, DIMENSION(2, 3) :: bounds
1886  INTEGER, DIMENSION(:), POINTER :: exat_neighbors
1887  LOGICAL :: do_gga, do_sc, do_sf
1888  TYPE(batch_info_type) :: batch_info
1889  TYPE(cp_1d_r_p_type), DIMENSION(:, :), POINTER :: ri_dcoeff
1890  TYPE(dft_control_type), POINTER :: dft_control
1891  TYPE(mp_para_env_type), POINTER :: para_env
1892  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1893  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1894  TYPE(section_vals_type), POINTER :: input, xc_functionals
1895  TYPE(xc_derivative_set_type) :: deriv_set
1896  TYPE(xc_rho_cflags_type) :: needs
1897  TYPE(xc_rho_set_type) :: rho_set
1898 
1899  NULLIFY (particle_set, qs_kind_set, dft_control, para_env, exat_neighbors)
1900  NULLIFY (input, xc_functionals)
1901 
1902  CALL timeset(routinen, handle)
1903 
1904 ! Initialize
1905  CALL get_qs_env(qs_env, particle_set=particle_set, qs_kind_set=qs_kind_set, natom=natom, &
1906  dft_control=dft_control, input=input, para_env=para_env)
1907  ALLOCATE (int_fxc(natom, 4))
1908  DO iatom = 1, natom
1909  DO i = 1, 4
1910  NULLIFY (int_fxc(iatom, i)%array)
1911  END DO
1912  END DO
1913  nex_atom = SIZE(xas_atom_env%excited_atoms)
1914  !spin conserving in the general sense here
1915  do_sc = xas_tdp_control%do_spin_cons .OR. xas_tdp_control%do_singlet .OR. xas_tdp_control%do_triplet
1916  do_sf = xas_tdp_control%do_spin_flip
1917 
1918 ! Get some info on the functionals
1919  xc_functionals => section_vals_get_subs_vals(input, "DFT%XAS_TDP%KERNEL%XC_FUNCTIONAL")
1920  ! ask for lsd in any case
1921  needs = xc_functionals_get_needs(xc_functionals, lsd=.true., calc_potential=.true.)
1922  do_gga = needs%drho_spin !because either LDA or GGA, and the former does not need gradient
1923 
1924 ! Distribute the excited atoms over batches of processors
1925 ! Then, the spherical orbital of the RI basis are distributed among the procs of the batch, making
1926 ! the GGA integration very efficient
1927  num_pe = para_env%num_pe
1928  mepos = para_env%mepos
1929 
1930  !create a batch_info_type
1931  CALL get_proc_batch_sizes(batch_size, nbatch, nex_atom, num_pe)
1932 
1933  !the batch index
1934  ibatch = mepos/batch_size
1935  !the proc index within the batch
1936  ipe = modulo(mepos, batch_size)
1937 
1938  batch_info%batch_size = batch_size
1939  batch_info%nbatch = nbatch
1940  batch_info%ibatch = ibatch
1941  batch_info%ipe = ipe
1942 
1943  !create a subcommunicator for this batch
1944  CALL batch_info%para_env%from_split(para_env, ibatch)
1945 
1946 ! Precompute the spherical orbital of the RI basis (and maybe their gradient) on the grids of the
1947 ! excited atoms. Needed for the GGA integration and to actually put the density on the grid
1948  CALL precompute_so_dso(do_gga, batch_info, xas_atom_env, qs_env)
1949 
1950  !distribute the excted atoms over the batches
1951  ex_bo = get_limit(nex_atom, nbatch, ibatch)
1952 
1953 ! Looping over the excited atoms
1954  DO iex = ex_bo(1), ex_bo(2)
1955 
1956  iatom = xas_atom_env%excited_atoms(iex)
1957  ikind = particle_set(iatom)%atomic_kind%kind_number
1958  exat_neighbors => xas_atom_env%exat_neighbors(iex)%array
1959  ri_dcoeff => xas_atom_env%ri_dcoeff(:, :, iex)
1960 
1961 ! General grid/basis info
1962  na = xas_atom_env%grid_atom_set(ikind)%grid_atom%ng_sphere
1963  nr = xas_atom_env%grid_atom_set(ikind)%grid_atom%nr
1964 
1965 ! Creating a xc_rho_set to store the density and dset for the kernel
1966  bounds(1:2, 1:3) = 1
1967  bounds(2, 1) = na
1968  bounds(2, 2) = nr
1969 
1970  CALL xc_rho_set_create(rho_set=rho_set, local_bounds=bounds, &
1971  rho_cutoff=dft_control%qs_control%eps_rho_rspace, &
1972  drho_cutoff=dft_control%qs_control%eps_rho_rspace)
1973  CALL xc_dset_create(deriv_set, local_bounds=bounds)
1974 
1975  ! allocate internals of the rho_set
1976  CALL xc_rho_set_atom_update(rho_set, needs, nspins=2, bo=bounds)
1977 
1978 ! Put the density, and possibly its gradient, on the grid (for this atom)
1979  CALL put_density_on_atomic_grid(rho_set, ri_dcoeff(iatom, :), ikind, &
1980  do_gga, batch_info, xas_atom_env, qs_env)
1981 
1982 ! Take the neighboring atom contributions to the density (and gradient)
1983 ! distribute the grid among the procs (for best load balance)
1984  nb_bo = get_limit(nr, batch_size, ipe)
1985  sr = nb_bo(1); er = nb_bo(2)
1986  DO inb = 1, SIZE(exat_neighbors)
1987 
1988  nb = exat_neighbors(inb)
1989  nbk = particle_set(nb)%atomic_kind%kind_number
1990  CALL put_density_on_other_grid(rho_set, ri_dcoeff(nb, :), nb, nbk, iatom, &
1991  ikind, sr, er, do_gga, xas_atom_env, qs_env)
1992 
1993  END DO
1994 
1995  ! make sure contributions from different procs are summed up
1996  CALL batch_info%para_env%sum(rho_set%rhoa)
1997  CALL batch_info%para_env%sum(rho_set%rhob)
1998  IF (do_gga) THEN
1999  DO dir = 1, 3
2000  CALL batch_info%para_env%sum(rho_set%drhoa(dir)%array)
2001  CALL batch_info%para_env%sum(rho_set%drhob(dir)%array)
2002  END DO
2003  END IF
2004 
2005 ! In case of GGA, also need the norm of the density gradient
2006  IF (do_gga) CALL compute_norm_drho(rho_set, ikind, xas_atom_env)
2007 
2008 ! Compute the required derivatives
2009  CALL xc_functionals_eval(xc_functionals, lsd=.true., rho_set=rho_set, deriv_set=deriv_set, &
2010  deriv_order=2)
2011 
2012  !spin-conserving (LDA part)
2013  IF (do_sc) THEN
2014  CALL integrate_sc_fxc(int_fxc, iatom, ikind, deriv_set, xas_atom_env, qs_env)
2015  END IF
2016 
2017  !spin-flip (LDA part)
2018  IF (do_sf) THEN
2019  CALL integrate_sf_fxc(int_fxc, iatom, ikind, rho_set, deriv_set, xas_atom_env, qs_env)
2020  END IF
2021 
2022  !Gradient correction (note: spin-flip only keeps the lda part, aka ALDA0)
2023  IF (do_gga .AND. do_sc) THEN
2024  CALL integrate_gga_fxc(int_fxc, iatom, ikind, batch_info, rho_set, deriv_set, &
2025  xas_atom_env, qs_env)
2026  END IF
2027 
2028 ! Clean-up
2029  CALL xc_dset_release(deriv_set)
2030  CALL xc_rho_set_release(rho_set)
2031  END DO !iex
2032 
2033  CALL release_batch_info(batch_info)
2034 
2035  !Not necessary to sync, but makes sure that any load inbalance is reported here
2036  CALL para_env%sync()
2037 
2038  CALL timestop(handle)
2039 
2040  END SUBROUTINE integrate_fxc_atoms
2041 
2042 ! **************************************************************************************************
2043 !> \brief Integrate the gradient correction part of the xc kernel on the atomic grid
2044 !> \param int_fxc the array containing the (P|fxc|Q) integrals
2045 !> \param iatom the index of the current excited atom
2046 !> \param ikind the index of the current excited kind
2047 !> \param batch_info how the so are distributed over the processor batch
2048 !> \param rho_set the variable contatinind the density and its gradient
2049 !> \param deriv_set the functional derivatives
2050 !> \param xas_atom_env ...
2051 !> \param qs_env ...
2052 !> \note Ignored in case of pure LDA, added on top of the LDA kernel in case of GGA
2053 ! **************************************************************************************************
2054  SUBROUTINE integrate_gga_fxc(int_fxc, iatom, ikind, batch_info, rho_set, deriv_set, &
2055  xas_atom_env, qs_env)
2056 
2057  TYPE(cp_2d_r_p_type), DIMENSION(:, :), POINTER :: int_fxc
2058  INTEGER, INTENT(IN) :: iatom, ikind
2059  TYPE(batch_info_type) :: batch_info
2060  TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
2061  TYPE(xc_derivative_set_type), INTENT(INOUT) :: deriv_set
2062  TYPE(xas_atom_env_type), POINTER :: xas_atom_env
2063  TYPE(qs_environment_type), POINTER :: qs_env
2064 
2065  CHARACTER(len=*), PARAMETER :: routinen = 'integrate_gga_fxc'
2066 
2067  INTEGER :: bo(2), dir, handle, i, ia, ir, jpgf, &
2068  jset, jso, l, maxso, na, nr, nset, &
2069  nsgf, nsoi, nsotot, startj, ub
2070  INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf
2071  REAL(dp), ALLOCATABLE, DIMENSION(:) :: rexp
2072  REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: int_sgf, res, so, work
2073  REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: dso
2074  REAL(dp), DIMENSION(:, :), POINTER :: dgr1, dgr2, ga, gr, ri_sphi_so, weight, &
2075  zet
2076  REAL(dp), DIMENSION(:, :, :), POINTER :: dga1, dga2
2077  TYPE(cp_2d_r_p_type), ALLOCATABLE, DIMENSION(:) :: int_so, vxc
2078  TYPE(cp_3d_r_p_type), ALLOCATABLE, DIMENSION(:) :: vxg
2079  TYPE(grid_atom_type), POINTER :: grid_atom
2080  TYPE(gto_basis_set_type), POINTER :: ri_basis
2081  TYPE(harmonics_atom_type), POINTER :: harmonics
2082  TYPE(mp_para_env_type), POINTER :: para_env
2083  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2084 
2085  NULLIFY (grid_atom, ri_basis, qs_kind_set, ga, gr, dgr1, dgr2, lmax, lmin, npgf)
2086  NULLIFY (weight, ri_sphi_so, dga1, dga2, para_env, harmonics, zet)
2087 
2088  !Strategy: we need to compute <phi_i|fxc|phij>, most of existing application of the 2nd
2089  ! functional derivative involve the response density, and the expression of the
2090  ! integral int (fxc*n^1) is well known. We substitute the spherical orbital phi_j
2091  ! in place of n^1 in the formula and thus perform the first integration. Then
2092  ! we obtain something in the form int (fxc*phi_j) = Vxc - div (Vxg) that we can
2093  ! put on the grid and treat like a potential. The second integration is done by
2094  ! using the divergence theorem and numerical integration:
2095  ! <phi_i|fxc|phi_j> = int phi_i*(Vxc - div(Vxg)) = int phi_i*Vxc + grad(phi_i).Vxg
2096  ! Note the sign change and the dot product.
2097 
2098  CALL timeset(routinen, handle)
2099 
2100  !If closed shell, only compute f_aa and f_ab (ub = 2)
2101  ub = 2
2102  IF (xas_atom_env%nspins == 2) ub = 3
2103 
2104  !Get the necessary grid info
2105  harmonics => xas_atom_env%harmonics_atom_set(ikind)%harmonics_atom
2106  grid_atom => xas_atom_env%grid_atom_set(ikind)%grid_atom
2107  na = grid_atom%ng_sphere
2108  nr = grid_atom%nr
2109  weight => grid_atom%weight
2110 
2111  !get the ri_basis info
2112  CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, para_env=para_env)
2113  CALL get_qs_kind(qs_kind_set(ikind), basis_set=ri_basis, basis_type="RI_XAS")
2114 
2115  CALL get_gto_basis_set(gto_basis_set=ri_basis, lmax=lmax, lmin=lmin, nsgf=nsgf, &
2116  maxso=maxso, npgf=npgf, nset=nset, zet=zet)
2117  nsotot = nset*maxso
2118 
2119  !Point to the precomputed so
2120  ga => xas_atom_env%ga(ikind)%array
2121  gr => xas_atom_env%gr(ikind)%array
2122  dgr1 => xas_atom_env%dgr1(ikind)%array
2123  dgr2 => xas_atom_env%dgr2(ikind)%array
2124  dga1 => xas_atom_env%dga1(ikind)%array
2125  dga2 => xas_atom_env%dga2(ikind)%array
2126 
2127  !Before integration, wanna pre-divide all relevant derivastives by the nrom of the gradient
2128  CALL divide_by_norm_drho(deriv_set, rho_set, lsd=.true.)
2129 
2130  !Wanna integrate <phi_i|fxc|phi_j>, start looping over phi_j and do the first integration, then
2131  !collect vxc and vxg and loop over phi_i for the second integration
2132  !Note: we do not use the CG coefficients because they are only useful when there is a product
2133  ! of Gaussians, which is not really the case here
2134  !Note: the spherical orbitals for phi_i are distributed among the prcos of the current batch
2135 
2136  ALLOCATE (so(na, nr))
2137  ALLOCATE (dso(na, nr, 3))
2138  ALLOCATE (rexp(nr))
2139 
2140  ALLOCATE (vxc(ub))
2141  ALLOCATE (vxg(ub))
2142  ALLOCATE (int_so(ub))
2143  DO i = 1, ub
2144  ALLOCATE (vxc(i)%array(na, nr))
2145  ALLOCATE (vxg(i)%array(na, nr, 3))
2146  ALLOCATE (int_so(i)%array(nsotot, nsotot))
2147  vxc(i)%array = 0.0_dp; vxg(i)%array = 0.0_dp; int_so(i)%array = 0.0_dp
2148  END DO
2149 
2150  DO jset = 1, nset
2151  DO jpgf = 1, npgf(jset)
2152  startj = (jset - 1)*maxso + (jpgf - 1)*nsoset(lmax(jset))
2153  DO jso = nsoset(lmin(jset) - 1) + 1, nsoset(lmax(jset))
2154  l = indso(1, jso)
2155 
2156  !put the so phi_j and its gradient on the grid
2157  !more efficient to recompute it rather than mp_bcast each chunk
2158 
2159  rexp(1:nr) = exp(-zet(jpgf, jset)*grid_atom%rad2(1:nr))
2160 !$OMP PARALLEL DO COLLAPSE(2) DEFAULT(NONE), &
2161 !$OMP SHARED(nr,na,so,dso,grid_atom,l,rexp,harmonics,jso,zet,jset,jpgf), &
2162 !$OMP PRIVATE(ir,ia,dir)
2163  DO ir = 1, nr
2164  DO ia = 1, na
2165 
2166  !so
2167  so(ia, ir) = grid_atom%rad(ir)**l*rexp(ir)*harmonics%slm(ia, jso)
2168 
2169  !dso
2170  dso(ia, ir, :) = 0.0_dp
2171  DO dir = 1, 3
2172  dso(ia, ir, dir) = dso(ia, ir, dir) &
2173  + grid_atom%rad(ir)**(l - 1)*rexp(ir)*harmonics%dslm_dxyz(dir, ia, jso) &
2174  - 2.0_dp*zet(jpgf, jset)*grid_atom%rad(ir)**(l + 1)*rexp(ir) &
2175  *harmonics%a(dir, ia)*harmonics%slm(ia, jso)
2176  END DO
2177  END DO
2178  END DO
2179 !$OMP END PARALLEL DO
2180 
2181  !Perform the first integration (analytically)
2182  CALL get_vxc_vxg(vxc, vxg, so, dso, na, nr, rho_set, deriv_set, weight)
2183 
2184  !For a given phi_j, compute the second integration with all phi_i at once
2185  !=> allows for efficient gemm to take place, especially since so are distributed
2186  nsoi = batch_info%nso_proc(ikind)
2187  bo = batch_info%so_bo(:, ikind)
2188  ALLOCATE (res(nsoi, nsoi), work(na, nsoi))
2189  res = 0.0_dp; work = 0.0_dp
2190 
2191  DO i = 1, ub
2192 
2193  !integrate so*Vxc and store in the int_so
2194  CALL dgemm('N', 'N', na, nsoi, nr, 1.0_dp, vxc(i)%array(:, :), na, &
2195  gr(:, 1:nsoi), nr, 0.0_dp, work, na)
2196  CALL dgemm('T', 'N', nsoi, nsoi, na, 1.0_dp, work, na, &
2197  ga(:, 1:nsoi), na, 0.0_dp, res, nsoi)
2198  int_so(i)%array(bo(1):bo(2), startj + jso) = get_diag(res)
2199 
2200  DO dir = 1, 3
2201 
2202  ! integrate and sum up Vxg*dso
2203  CALL dgemm('N', 'N', na, nsoi, nr, 1.0_dp, vxg(i)%array(:, :, dir), na, &
2204  dgr1(:, 1:nsoi), nr, 0.0_dp, work, na)
2205  CALL dgemm('T', 'N', nsoi, nsoi, na, 1.0_dp, work, na, &
2206  dga1(:, 1:nsoi, dir), na, 0.0_dp, res, nsoi)
2207  CALL daxpy(nsoi, 1.0_dp, get_diag(res), 1, int_so(i)%array(bo(1):bo(2), startj + jso), 1)
2208 
2209  CALL dgemm('N', 'N', na, nsoi, nr, 1.0_dp, vxg(i)%array(:, :, dir), na, &
2210  dgr2(:, 1:nsoi), nr, 0.0_dp, work, na)
2211  CALL dgemm('T', 'N', nsoi, nsoi, na, 1.0_dp, work, na, &
2212  dga2(:, 1:nsoi, dir), na, 0.0_dp, res, nsoi)
2213  CALL daxpy(nsoi, 1.0_dp, get_diag(res), 1, int_so(i)%array(bo(1):bo(2), startj + jso), 1)
2214 
2215  END DO
2216 
2217  END DO !i
2218  DEALLOCATE (res, work)
2219 
2220  END DO !jso
2221  END DO !jpgf
2222  END DO !jset
2223 
2224  !Contract into sgf and add to already computed LDA part of int_fxc
2225  ri_sphi_so => xas_atom_env%ri_sphi_so(ikind)%array
2226  ALLOCATE (int_sgf(nsgf, nsgf))
2227  DO i = 1, ub
2228  CALL batch_info%para_env%sum(int_so(i)%array)
2229  CALL contract_so2sgf(int_sgf, int_so(i)%array, ri_basis, ri_sphi_so)
2230  CALL daxpy(nsgf*nsgf, 1.0_dp, int_sgf, 1, int_fxc(iatom, i)%array, 1)
2231  END DO
2232 
2233  !Clean-up
2234  DO i = 1, ub
2235  DEALLOCATE (vxc(i)%array)
2236  DEALLOCATE (vxg(i)%array)
2237  DEALLOCATE (int_so(i)%array)
2238  END DO
2239  DEALLOCATE (vxc, vxg, int_so)
2240 
2241  CALL timestop(handle)
2242 
2243  END SUBROUTINE integrate_gga_fxc
2244 
2245 ! **************************************************************************************************
2246 !> \brief Computes the first integration of the GGA part of <phi_i|fxc|phi_j>, i.e. int fxc*phi_j.
2247 !> The result is of the form Vxc - div(Vxg). Up to 3 results are returned, correspoinding to
2248 !> f_aa, f_ab and (if open-shell) f_bb
2249 !> \param vxc ...
2250 !> \param vxg ...
2251 !> \param so the spherical orbital on the grid
2252 !> \param dso the derivative of the spherical orbital on the grid
2253 !> \param na ...
2254 !> \param nr ...
2255 !> \param rho_set ...
2256 !> \param deriv_set ...
2257 !> \param weight the grid weight
2258 !> \note This method is extremely similar to xc_calc_2nd_deriv of xc.F, but because it is a special
2259 !> case that can be further optimized and because the interface of the original routine does
2260 !> not fit this code, it has been re-written (no pw, no rho1_set but just the so, etc...)
2261 ! **************************************************************************************************
2262  SUBROUTINE get_vxc_vxg(vxc, vxg, so, dso, na, nr, rho_set, deriv_set, weight)
2263 
2264  TYPE(cp_2d_r_p_type), DIMENSION(:) :: vxc
2265  TYPE(cp_3d_r_p_type), DIMENSION(:) :: vxg
2266  REAL(dp), DIMENSION(:, :) :: so
2267  REAL(dp), DIMENSION(:, :, :) :: dso
2268  INTEGER, INTENT(IN) :: na, nr
2269  TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
2270  TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
2271  REAL(dp), DIMENSION(:, :), POINTER :: weight
2272 
2273  CHARACTER(len=*), PARAMETER :: routinen = 'get_vxc_vxg'
2274 
2275  INTEGER :: dir, handle, i, ia, ir, ub
2276  REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: dot_proda, dot_prodb, tmp
2277  REAL(dp), DIMENSION(:, :, :), POINTER :: d1e, d2e, norm_drhoa, norm_drhob
2278  TYPE(xc_derivative_type), POINTER :: deriv
2279 
2280  NULLIFY (norm_drhoa, norm_drhob, d2e, d1e, deriv)
2281 
2282  CALL timeset(routinen, handle)
2283 
2284  !Note: this routines follows the order of the terms in equaiton (A.7) of Thomas Chassaing
2285  ! thesis, except for the pure LDA terms that are dropped. The n^(1)_a and n^(1)_b
2286  ! response densities are replaced by the spherical orbital.
2287  ! The usual spin ordering is used: aa => 1, ab => 2 , bb => 3
2288 
2289  !point to the relevant components of rho_set
2290  ub = SIZE(vxc)
2291  norm_drhoa => rho_set%norm_drhoa
2292  norm_drhob => rho_set%norm_drhob
2293 
2294  !Some init
2295  DO i = 1, ub
2296  vxc(i)%array = 0.0_dp
2297  vxg(i)%array = 0.0_dp
2298  END DO
2299 
2300  ALLOCATE (tmp(na, nr), dot_proda(na, nr), dot_prodb(na, nr))
2301  dot_proda = 0.0_dp; dot_prodb = 0.0_dp
2302 
2303  !Strategy: most terms are either multiplied by drhoa or drhob => group those first and then
2304  ! multiply. Also most terms are multiplied by the dot product grad_n . grad_so, so
2305  ! precompute it as well
2306 
2307 !$OMP PARALLEL DEFAULT(NONE), &
2308 !$OMP SHARED(dot_proda,dot_prodb,tmp,vxc,vxg,deriv_set,rho_set,na,nr,norm_drhoa,norm_drhob,dso, &
2309 !$OMP so,weight,ub), &
2310 !$OMP PRIVATE(ia,ir,dir,deriv,d1e,d2e)
2311 
2312  !Precompute the very common dot products grad_na . grad_so and grad_nb . grad_so
2313  DO dir = 1, 3
2314 !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2315  DO ir = 1, nr
2316  DO ia = 1, na
2317  dot_proda(ia, ir) = dot_proda(ia, ir) + rho_set%drhoa(dir)%array(ia, ir, 1)*dso(ia, ir, dir)
2318  dot_prodb(ia, ir) = dot_prodb(ia, ir) + rho_set%drhob(dir)%array(ia, ir, 1)*dso(ia, ir, dir)
2319  END DO !ia
2320  END DO !ir
2321 !$OMP END DO NOWAIT
2322  END DO !dir
2323 
2324  !Deal with f_aa
2325 
2326  !Vxc, first term
2327  deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_norm_drhoa])
2328  IF (ASSOCIATED(deriv)) THEN
2329  CALL xc_derivative_get(deriv, deriv_data=d2e)
2330 !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2331  DO ir = 1, nr
2332  DO ia = 1, na
2333  vxc(1)%array(ia, ir) = d2e(ia, ir, 1)*dot_proda(ia, ir)
2334  END DO !ia
2335  END DO !ir
2336 !$OMP END DO NOWAIT
2337  END IF
2338 
2339  !Vxc, second term
2340  deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_norm_drho])
2341 
2342  IF (ASSOCIATED(deriv)) THEN
2343  CALL xc_derivative_get(deriv, deriv_data=d2e)
2344 !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2345  DO ir = 1, nr
2346  DO ia = 1, na
2347  vxc(1)%array(ia, ir) = vxc(1)%array(ia, ir) + d2e(ia, ir, 1)*dot_prodb(ia, ir)
2348  END DO !ia
2349  END DO !ir
2350 !$OMP END DO NOWAIT
2351  END IF
2352 
2353  !Vxc, take the grid weight into acocunt
2354 !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2355  DO ir = 1, nr
2356  DO ia = 1, na
2357  vxc(1)%array(ia, ir) = vxc(1)%array(ia, ir)*weight(ia, ir)
2358  END DO !ia
2359  END DO !ir
2360 !$OMP END DO NOWAIT
2361 
2362  !Vxg, first term (to be multiplied by drhoa)
2363  deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_norm_drhoa])
2364  IF (ASSOCIATED(deriv)) THEN
2365  CALL xc_derivative_get(deriv, deriv_data=d2e)
2366 !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2367  DO ir = 1, nr
2368  DO ia = 1, na
2369  tmp(ia, ir) = d2e(ia, ir, 1)*so(ia, ir)
2370  END DO !ia
2371  END DO !ir
2372 !$OMP END DO NOWAIT
2373  END IF
2374 
2375  !Vxg, second term (to be multiplied by drhoa)
2376  deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_norm_drho])
2377  IF (ASSOCIATED(deriv)) THEN
2378  CALL xc_derivative_get(deriv, deriv_data=d2e)
2379 !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2380  DO ir = 1, nr
2381  DO ia = 1, na
2382  tmp(ia, ir) = tmp(ia, ir) + d2e(ia, ir, 1)*dot_prodb(ia, ir)
2383  END DO !ia
2384  END DO !ir
2385 !$OMP END DO NOWAIT
2386  END IF
2387 
2388  !Vxg, third term (to be multiplied by drhoa)
2389  deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_norm_drhoa])
2390  IF (ASSOCIATED(deriv)) THEN
2391  CALL xc_derivative_get(deriv, deriv_data=d2e)
2392 !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2393  DO ir = 1, nr
2394  DO ia = 1, na
2395  tmp(ia, ir) = tmp(ia, ir) + d2e(ia, ir, 1)*dot_proda(ia, ir)
2396  END DO !ia
2397  END DO !ir
2398 !$OMP END DO NOWAIT
2399  END IF
2400 
2401  !Vxg, fourth term (to be multiplied by drhoa)
2402  deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa])
2403  IF (ASSOCIATED(deriv)) THEN
2404  CALL xc_derivative_get(deriv, deriv_data=d1e)
2405 !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2406  DO ir = 1, nr
2407  DO ia = 1, na
2408  tmp(ia, ir) = tmp(ia, ir) - d1e(ia, ir, 1)*dot_proda(ia, ir) &
2409  /max(norm_drhoa(ia, ir, 1), rho_set%drho_cutoff)**2
2410  END DO !ia
2411  END DO !ir
2412 !$OMP END DO NOWAIT
2413  END IF
2414 
2415  !put tmp*drhoa in Vxg (so that we can reuse it for drhob terms)
2416  DO dir = 1, 3
2417 !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2418  DO ir = 1, nr
2419  DO ia = 1, na
2420  vxg(1)%array(ia, ir, dir) = tmp(ia, ir)*rho_set%drhoa(dir)%array(ia, ir, 1)
2421  END DO !ia
2422  END DO !ir
2423 !$OMP END DO NOWAIT
2424  END DO !dir
2425 
2426  !Vxg, fifth term (to be multiplied by drhob)
2427  deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_norm_drho])
2428  IF (ASSOCIATED(deriv)) THEN
2429  CALL xc_derivative_get(deriv, deriv_data=d2e)
2430 !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2431  DO ir = 1, nr
2432  DO ia = 1, na
2433  tmp(ia, ir) = d2e(ia, ir, 1)*so(ia, ir)
2434  END DO !ia
2435  END DO !ir
2436 !$OMP END DO NOWAIT
2437  END IF
2438 
2439  !Vxg, sixth term (to be multiplied by drhob)
2440  deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_norm_drho])
2441  IF (ASSOCIATED(deriv)) THEN
2442  CALL xc_derivative_get(deriv, deriv_data=d2e)
2443 !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2444  DO ir = 1, nr
2445  DO ia = 1, na
2446  tmp(ia, ir) = tmp(ia, ir) + d2e(ia, ir, 1)*dot_proda(ia, ir)
2447  END DO !ia
2448  END DO !ir
2449 !$OMP END DO NOWAIT
2450  END IF
2451 
2452  !Vxg, seventh term (to be multiplied by drhob)
2453  deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_norm_drho])
2454  IF (ASSOCIATED(deriv)) THEN
2455  CALL xc_derivative_get(deriv, deriv_data=d2e)
2456 !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2457  DO ir = 1, nr
2458  DO ia = 1, na
2459  tmp(ia, ir) = tmp(ia, ir) + d2e(ia, ir, 1)*dot_prodb(ia, ir)
2460  END DO !ia
2461  END DO !ir
2462 !$OMP END DO NOWAIT
2463  END IF
2464 
2465  !put tmp*drhob in Vxg
2466  DO dir = 1, 3
2467 !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2468  DO ir = 1, nr
2469  DO ia = 1, na
2470  vxg(1)%array(ia, ir, dir) = vxg(1)%array(ia, ir, dir) + tmp(ia, ir)*rho_set%drhob(dir)%array(ia, ir, 1)
2471  END DO !ia
2472  END DO !ir
2473 !$OMP END DO NOWAIT
2474  END DO !dir
2475 
2476  !Vxg, last term
2477  deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa])
2478  IF (ASSOCIATED(deriv)) THEN
2479  CALL xc_derivative_get(deriv, deriv_data=d1e)
2480  DO dir = 1, 3
2481 !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2482  DO ir = 1, nr
2483  DO ia = 1, na
2484  vxg(1)%array(ia, ir, dir) = vxg(1)%array(ia, ir, dir) + d1e(ia, ir, 1)*dso(ia, ir, dir)
2485  END DO !ia
2486  END DO !ir
2487 !$OMP END DO NOWAIT
2488  END DO !dir
2489  END IF
2490 
2491  !Vxg, take the grid weight into account
2492  DO dir = 1, 3
2493 !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2494  DO ir = 1, nr
2495  DO ia = 1, na
2496  vxg(1)%array(ia, ir, dir) = vxg(1)%array(ia, ir, dir)*weight(ia, ir)
2497  END DO !ia
2498  END DO !ir
2499 !$OMP END DO NOWAIT
2500  END DO !dir
2501 
2502  !Deal with fab
2503 
2504  !Vxc, first term
2505  deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_norm_drhob])
2506  IF (ASSOCIATED(deriv)) THEN
2507  CALL xc_derivative_get(deriv, deriv_data=d2e)
2508 !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2509  DO ir = 1, nr
2510  DO ia = 1, na
2511  vxc(2)%array(ia, ir) = d2e(ia, ir, 1)*dot_prodb(ia, ir)
2512  END DO !ia
2513  END DO !ir
2514 !$OMP END DO NOWAIT
2515  END IF
2516 
2517  !Vxc, second term
2518  deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_norm_drho])
2519  IF (ASSOCIATED(deriv)) THEN
2520  CALL xc_derivative_get(deriv, deriv_data=d2e)
2521 !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2522  DO ir = 1, nr
2523  DO ia = 1, na
2524  vxc(2)%array(ia, ir) = vxc(2)%array(ia, ir) + d2e(ia, ir, 1)*dot_proda(ia, ir)
2525  END DO !ia
2526  END DO !ir
2527 !$OMP END DO NOWAIT
2528  END IF
2529 
2530  !Vxc, take the grid weight into acocunt
2531 !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2532  DO ir = 1, nr
2533  DO ia = 1, na
2534  vxc(2)%array(ia, ir) = vxc(2)%array(ia, ir)*weight(ia, ir)
2535  END DO !ia
2536  END DO !ir
2537 !$OMP END DO NOWAIT
2538 
2539  !Vxg, first term (to be multiplied by drhoa)
2540  deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_norm_drhoa])
2541  IF (ASSOCIATED(deriv)) THEN
2542  CALL xc_derivative_get(deriv, deriv_data=d2e)
2543 !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2544  DO ir = 1, nr
2545  DO ia = 1, na
2546  tmp(ia, ir) = d2e(ia, ir, 1)*so(ia, ir)
2547  END DO !ia
2548  END DO !ir
2549 !$OMP END DO NOWAIT
2550  END IF
2551 
2552  !Vxg, second term (to be multiplied by drhoa)
2553  deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_norm_drho])
2554  IF (ASSOCIATED(deriv)) THEN
2555  CALL xc_derivative_get(deriv, deriv_data=d2e)
2556 !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2557  DO ir = 1, nr
2558  DO ia = 1, na
2559  tmp(ia, ir) = tmp(ia, ir) + d2e(ia, ir, 1)*dot_proda(ia, ir)
2560  END DO !ia
2561  END DO !ir
2562 !$OMP END DO NOWAIT
2563  END IF
2564 
2565  !Vxg, third term (to be multiplied by drhoa)
2566  deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_norm_drhob])
2567  IF (ASSOCIATED(deriv)) THEN
2568  CALL xc_derivative_get(deriv, deriv_data=d2e)
2569 !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2570  DO ir = 1, nr
2571  DO ia = 1, na
2572  tmp(ia, ir) = tmp(ia, ir) + d2e(ia, ir, 1)*dot_prodb(ia, ir)
2573  END DO !ia
2574  END DO !ir
2575 !$OMP END DO NOWAIT
2576  END IF
2577 
2578  !put tmp*drhoa in Vxg
2579  DO dir = 1, 3
2580 !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2581  DO ir = 1, nr
2582  DO ia = 1, na
2583  vxg(2)%array(ia, ir, dir) = tmp(ia, ir)*rho_set%drhoa(dir)%array(ia, ir, 1)
2584  END DO !ia
2585  END DO !ir
2586 !$OMP END DO NOWAIT
2587  END DO !dir
2588 
2589  !Vxg, fourth term (to be multiplied by drhob)
2590  deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_norm_drho])
2591  IF (ASSOCIATED(deriv)) THEN
2592  CALL xc_derivative_get(deriv, deriv_data=d2e)
2593 !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2594  DO ir = 1, nr
2595  DO ia = 1, na
2596  tmp(ia, ir) = d2e(ia, ir, 1)*so(ia, ir)
2597  END DO !ia
2598  END DO !ir
2599 !$OMP END DO NOWAIT
2600  END IF
2601 
2602  !Vxg, fifth term (to be multiplied by drhob)
2603  deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_norm_drhob])
2604  IF (ASSOCIATED(deriv)) THEN
2605  CALL xc_derivative_get(deriv, deriv_data=d2e)
2606 !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2607  DO ir = 1, nr
2608  DO ia = 1, na
2609  tmp(ia, ir) = tmp(ia, ir) + d2e(ia, ir, 1)*dot_prodb(ia, ir)
2610  END DO !ia
2611  END DO !ir
2612 !$OMP END DO NOWAIT
2613  END IF
2614 
2615  !Vxg, sixth term (to be multiplied by drhob)
2616  deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_norm_drho])
2617  IF (ASSOCIATED(deriv)) THEN
2618  CALL xc_derivative_get(deriv, deriv_data=d2e)
2619 !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2620  DO ir = 1, nr
2621  DO ia = 1, na
2622  tmp(ia, ir) = tmp(ia, ir) + d2e(ia, ir, 1)*dot_proda(ia, ir)
2623  END DO !ia
2624  END DO !ir
2625 !$OMP END DO NOWAIT
2626  END IF
2627 
2628  !put tmp*drhob in Vxg
2629  DO dir = 1, 3
2630 !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2631  DO ir = 1, nr
2632  DO ia = 1, na
2633  vxg(2)%array(ia, ir, dir) = vxg(2)%array(ia, ir, dir) + tmp(ia, ir)*rho_set%drhob(dir)%array(ia, ir, 1)
2634  END DO
2635  END DO
2636 !$OMP END DO NOWAIT
2637  END DO
2638 
2639  !Vxg, last term
2640  deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho])
2641  IF (ASSOCIATED(deriv)) THEN
2642  CALL xc_derivative_get(deriv, deriv_data=d1e)
2643  DO dir = 1, 3
2644 !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2645  DO ir = 1, nr
2646  DO ia = 1, na
2647  vxg(2)%array(ia, ir, dir) = vxg(2)%array(ia, ir, dir) + d1e(ia, ir, 1)*dso(ia, ir, dir)
2648  END DO !ia
2649  END DO !ir
2650 !$OMP END DO NOWAIT
2651  END DO !dir
2652  END IF
2653 
2654  !Vxg, take the grid weight into account
2655  DO dir = 1, 3
2656 !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2657  DO ir = 1, nr
2658  DO ia = 1, na
2659  vxg(2)%array(ia, ir, dir) = vxg(2)%array(ia, ir, dir)*weight(ia, ir)
2660  END DO !ia
2661  END DO !ir
2662 !$OMP END DO NOWAIT
2663  END DO !dir
2664 
2665  !Deal with f_bb, if so required
2666  IF (ub == 3) THEN
2667 
2668  !Vxc, first term
2669  deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_norm_drhob])
2670  IF (ASSOCIATED(deriv)) THEN
2671  CALL xc_derivative_get(deriv, deriv_data=d2e)
2672 !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2673  DO ir = 1, nr
2674  DO ia = 1, na
2675  vxc(3)%array(ia, ir) = d2e(ia, ir, 1)*dot_prodb(ia, ir)
2676  END DO !ia
2677  END DO !ir
2678 !$OMP END DO NOWAIT
2679  END IF
2680 
2681  !Vxc, second term
2682  deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_norm_drho])
2683  IF (ASSOCIATED(deriv)) THEN
2684  CALL xc_derivative_get(deriv, deriv_data=d2e)
2685 !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2686  DO ir = 1, nr
2687  DO ia = 1, na
2688  vxc(3)%array(ia, ir) = vxc(3)%array(ia, ir) + d2e(ia, ir, 1)*dot_proda(ia, ir)
2689  END DO !i
2690  END DO !ir
2691 !$OMP END DO NOWAIT
2692  END IF
2693 
2694  !Vxc, take the grid weight into acocunt
2695 !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2696  DO ir = 1, nr
2697  DO ia = 1, na
2698  vxc(3)%array(ia, ir) = vxc(3)%array(ia, ir)*weight(ia, ir)
2699  END DO !ia
2700  END DO !ir
2701 !$OMP END DO NOWAIT
2702 
2703  !Vxg, first term (to be multiplied by drhob)
2704  deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_norm_drhob])
2705  IF (ASSOCIATED(deriv)) THEN
2706  CALL xc_derivative_get(deriv, deriv_data=d2e)
2707 !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2708  DO ir = 1, nr
2709  DO ia = 1, na
2710  tmp(ia, ir) = d2e(ia, ir, 1)*so(ia, ir)
2711  END DO !ia
2712  END DO !ir
2713 !$OMP END DO NOWAIT
2714  END IF
2715 
2716  !Vxg, second term (to be multiplied by drhob)
2717  deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob, deriv_norm_drho])
2718  IF (ASSOCIATED(deriv)) THEN
2719  CALL xc_derivative_get(deriv, deriv_data=d2e)
2720 !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2721  DO ir = 1, nr
2722  DO ia = 1, na
2723  tmp(ia, ir) = tmp(ia, ir) + d2e(ia, ir, 1)*dot_proda(ia, ir)
2724  END DO !ia
2725  END DO !ir
2726 !$OMP END DO NOWAIT
2727  END IF
2728 
2729  !Vxg, third term (to be multiplied by drhob)
2730  deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob, deriv_norm_drhob])
2731  IF (ASSOCIATED(deriv)) THEN
2732  CALL xc_derivative_get(deriv, deriv_data=d2e)
2733 !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2734  DO ir = 1, nr
2735  DO ia = 1, na
2736  tmp(ia, ir) = tmp(ia, ir) + d2e(ia, ir, 1)*dot_prodb(ia, ir)
2737  END DO !ia
2738  END DO !ir
2739 !$OMP END DO NOWAIT
2740  END IF
2741 
2742  !Vxg, fourth term (to be multiplied by drhob)
2743  deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob])
2744  IF (ASSOCIATED(deriv)) THEN
2745  CALL xc_derivative_get(deriv, deriv_data=d1e)
2746 !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2747  DO ir = 1, nr
2748  DO ia = 1, na
2749  tmp(ia, ir) = tmp(ia, ir) - d1e(ia, ir, 1)*dot_prodb(ia, ir) &
2750  /max(norm_drhob(ia, ir, 1), rho_set%drho_cutoff)**2
2751  END DO !ia
2752  END DO !ir
2753 !$OMP END DO NOWAIT
2754  END IF
2755 
2756  !put tmp*drhob in Vxg (so that we can reuse it for drhoa terms)
2757  DO dir = 1, 3
2758 !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2759  DO ir = 1, nr
2760  DO ia = 1, na
2761  vxg(3)%array(ia, ir, dir) = tmp(ia, ir)*rho_set%drhob(dir)%array(ia, ir, 1)
2762  END DO !ia
2763  END DO !ir
2764 !$OMP END DO NOWAIT
2765  END DO !dir
2766 
2767  !Vxg, fifth term (to be multiplied by drhoa)
2768  deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_norm_drho])
2769  IF (ASSOCIATED(deriv)) THEN
2770  CALL xc_derivative_get(deriv, deriv_data=d2e)
2771 !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2772  DO ir = 1, nr
2773  DO ia = 1, na
2774  tmp(ia, ir) = d2e(ia, ir, 1)*so(ia, ir)
2775  END DO !ia
2776  END DO !ir
2777 !$OMP END DO NOWAIT
2778  END IF
2779 
2780  !Vxg, sixth term (to be multiplied by drhoa)
2781  deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob, deriv_norm_drho])
2782  IF (ASSOCIATED(deriv)) THEN
2783  CALL xc_derivative_get(deriv, deriv_data=d2e)
2784 !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2785  DO ir = 1, nr
2786  DO ia = 1, na
2787  tmp(ia, ir) = tmp(ia, ir) + d2e(ia, ir, 1)*dot_prodb(ia, ir)
2788  END DO !ia
2789  END DO !ir
2790 !$OMP END DO NOWAIT
2791  END IF
2792 
2793  !Vxg, seventh term (to be multiplied by drhoa)
2794  deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_norm_drho])
2795  IF (ASSOCIATED(deriv)) THEN
2796  CALL xc_derivative_get(deriv, deriv_data=d2e)
2797 !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2798  DO ir = 1, nr
2799  DO ia = 1, na
2800  tmp(ia, ir) = tmp(ia, ir) + d2e(ia, ir, 1)*dot_proda(ia, ir)
2801  END DO !ia
2802  END DO !ir
2803 !$OMP END DO NOWAIT
2804  END IF
2805 
2806  !put tmp*drhoa in Vxg
2807  DO dir = 1, 3
2808 !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2809  DO ir = 1, nr
2810  DO ia = 1, na
2811  vxg(3)%array(ia, ir, dir) = vxg(3)%array(ia, ir, dir) + &
2812  tmp(ia, ir)*rho_set%drhoa(dir)%array(ia, ir, 1)
2813  END DO !ia
2814  END DO !ir
2815 !$OMP END DO NOWAIT
2816  END DO !dir
2817 
2818  !Vxg, last term
2819  deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob])
2820  IF (ASSOCIATED(deriv)) THEN
2821  CALL xc_derivative_get(deriv, deriv_data=d1e)
2822  DO dir = 1, 3
2823 !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2824  DO ir = 1, nr
2825  DO ia = 1, na
2826  vxg(3)%array(ia, ir, dir) = vxg(3)%array(ia, ir, dir) + d1e(ia, ir, 1)*dso(ia, ir, dir)
2827  END DO !ia
2828  END DO !ir
2829 !$OMP END DO NOWAIT
2830  END DO !dir
2831  END IF
2832 
2833  !Vxg, take the grid weight into account
2834  DO dir = 1, 3
2835 !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2836  DO ir = 1, nr
2837  DO ia = 1, na
2838  vxg(3)%array(ia, ir, dir) = vxg(3)%array(ia, ir, dir)*weight(ia, ir)
2839  END DO !ia
2840  END DO !ir
2841 !$OMP END DO NOWAIT
2842  END DO !dir
2843 
2844  END IF !f_bb
2845 
2846 !$OMP END PARALLEL
2847 
2848  CALL timestop(handle)
2849 
2850  END SUBROUTINE get_vxc_vxg
2851 
2852 ! **************************************************************************************************
2853 !> \brief Integrate the fxc kernel in the spin-conserving case, be it closed- or open-shell
2854 !> \param int_fxc the array containing the (P|fxc|Q) integrals
2855 !> \param iatom the index of the current excited atom
2856 !> \param ikind the index of the current excited kind
2857 !> \param deriv_set the set of functional derivatives
2858 !> \param xas_atom_env ...
2859 !> \param qs_env ...
2860 ! **************************************************************************************************
2861  SUBROUTINE integrate_sc_fxc(int_fxc, iatom, ikind, deriv_set, xas_atom_env, qs_env)
2862 
2863  TYPE(cp_2d_r_p_type), DIMENSION(:, :), POINTER :: int_fxc
2864  INTEGER, INTENT(IN) :: iatom, ikind
2865  TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
2866  TYPE(xas_atom_env_type), POINTER :: xas_atom_env
2867  TYPE(qs_environment_type), POINTER :: qs_env
2868 
2869  INTEGER :: i, maxso, na, nr, nset, nsotot, nspins, &
2870  ri_nsgf, ub
2871  REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: fxc, int_so
2872  REAL(dp), DIMENSION(:, :), POINTER :: ri_sphi_so
2873  TYPE(cp_3d_r_p_type), ALLOCATABLE, DIMENSION(:) :: d2e
2874  TYPE(dft_control_type), POINTER :: dft_control
2875  TYPE(grid_atom_type), POINTER :: grid_atom
2876  TYPE(gto_basis_set_type), POINTER :: ri_basis
2877  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2878  TYPE(xc_derivative_type), POINTER :: deriv
2879 
2880  NULLIFY (grid_atom, deriv, ri_basis, ri_sphi_so, qs_kind_set, dft_control)
2881 
2882  ! Initialization
2883  CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, dft_control=dft_control)
2884  grid_atom => xas_atom_env%grid_atom_set(ikind)%grid_atom
2885  na = grid_atom%ng_sphere
2886  nr = grid_atom%nr
2887  CALL get_qs_kind(qs_kind_set(ikind), basis_set=ri_basis, basis_type="RI_XAS")
2888  CALL get_gto_basis_set(ri_basis, nset=nset, maxso=maxso, nsgf=ri_nsgf)
2889  nsotot = nset*maxso
2890  ri_sphi_so => xas_atom_env%ri_sphi_so(ikind)%array
2891  nspins = dft_control%nspins
2892 
2893  ! Get the second derivatives
2894  ALLOCATE (d2e(3))
2895  deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa])
2896  CALL xc_derivative_get(deriv, deriv_data=d2e(1)%array)
2897  deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhob])
2898  CALL xc_derivative_get(deriv, deriv_data=d2e(2)%array)
2899  deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob])
2900  CALL xc_derivative_get(deriv, deriv_data=d2e(3)%array)
2901 
2902  ! Allocate some work arrays
2903  ALLOCATE (fxc(na, nr))
2904  ALLOCATE (int_so(nsotot, nsotot))
2905 
2906  ! Integrate for all three derivatives, taking the grid weight into account
2907  ! If closed shell, do not need to integrate beta-beta as it is the same as alpha-alpha
2908  ub = 2; IF (nspins == 2) ub = 3
2909  DO i = 1, ub
2910 
2911  int_so = 0.0_dp
2912  fxc(:, :) = d2e(i)%array(:, :, 1)*grid_atom%weight(:, :)
2913  CALL integrate_so_prod(int_so, fxc, ikind, xas_atom_env, qs_env)
2914 
2915  !contract into sgf. Array allocated on current processor only
2916  ALLOCATE (int_fxc(iatom, i)%array(ri_nsgf, ri_nsgf))
2917  int_fxc(iatom, i)%array = 0.0_dp
2918  CALL contract_so2sgf(int_fxc(iatom, i)%array, int_so, ri_basis, ri_sphi_so)
2919 
2920  END DO
2921 
2922  END SUBROUTINE integrate_sc_fxc
2923 
2924 ! **************************************************************************************************
2925 !> \brief Integrate the fxc kernel in the spin-flip case (open-shell assumed). The spin-flip LDA
2926 !> kernel reads: fxc = 1/(rhoa - rhob) * (dE/drhoa - dE/drhob)
2927 !> \param int_fxc the array containing the (P|fxc|Q) integrals
2928 !> \param iatom the index of the current excited atom
2929 !> \param ikind the index of the current excited kind
2930 !> \param rho_set the density in the atomic grid
2931 !> \param deriv_set the set of functional derivatives
2932 !> \param xas_atom_env ...
2933 !> \param qs_env ...
2934 ! **************************************************************************************************
2935  SUBROUTINE integrate_sf_fxc(int_fxc, iatom, ikind, rho_set, deriv_set, xas_atom_env, qs_env)
2936 
2937  TYPE(cp_2d_r_p_type), DIMENSION(:, :), POINTER :: int_fxc
2938  INTEGER, INTENT(IN) :: iatom, ikind
2939  TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
2940  TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
2941  TYPE(xas_atom_env_type), POINTER :: xas_atom_env
2942  TYPE(qs_environment_type), POINTER :: qs_env
2943 
2944  INTEGER :: ia, ir, maxso, na, nr, nset, nsotot, &
2945  ri_nsgf
2946  REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: fxc, int_so
2947  REAL(dp), DIMENSION(:, :), POINTER :: ri_sphi_so
2948  REAL(dp), DIMENSION(:, :, :), POINTER :: rhoa, rhob
2949  TYPE(cp_3d_r_p_type), ALLOCATABLE, DIMENSION(:) :: d1e, d2e
2950  TYPE(dft_control_type), POINTER :: dft_control
2951  TYPE(grid_atom_type), POINTER :: grid_atom
2952  TYPE(gto_basis_set_type), POINTER :: ri_basis
2953  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2954  TYPE(xc_derivative_type), POINTER :: deriv
2955 
2956  NULLIFY (grid_atom, deriv, ri_basis, ri_sphi_so, qs_kind_set, rhoa, rhob, dft_control)
2957 
2958  ! Initialization
2959  CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, dft_control=dft_control)
2960  grid_atom => xas_atom_env%grid_atom_set(ikind)%grid_atom
2961  na = grid_atom%ng_sphere
2962  nr = grid_atom%nr
2963  CALL get_qs_kind(qs_kind_set(ikind), basis_set=ri_basis, basis_type="RI_XAS")
2964  CALL get_gto_basis_set(ri_basis, nset=nset, maxso=maxso, nsgf=ri_nsgf)
2965  nsotot = nset*maxso
2966  ri_sphi_so => xas_atom_env%ri_sphi_so(ikind)%array
2967  rhoa => rho_set%rhoa
2968  rhob => rho_set%rhob
2969 
2970  ALLOCATE (d1e(2))
2971  deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa])
2972  CALL xc_derivative_get(deriv, deriv_data=d1e(1)%array)
2973  deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob])
2974  CALL xc_derivative_get(deriv, deriv_data=d1e(2)%array)
2975 
2976  ! In case rhoa -> rhob, take the limit, which involves the (already computed) second derivatives
2977  ALLOCATE (d2e(3))
2978  deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa])
2979  CALL xc_derivative_get(deriv, deriv_data=d2e(1)%array)
2980  deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhob])
2981  CALL xc_derivative_get(deriv, deriv_data=d2e(2)%array)
2982  deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob])
2983  CALL xc_derivative_get(deriv, deriv_data=d2e(3)%array)
2984 
2985  !Compute the kernel on the grid. Already take weight into acocunt there
2986  ALLOCATE (fxc(na, nr))
2987 !$OMP PARALLEL DO COLLAPSE(2) SCHEDULE(STATIC) DEFAULT(NONE), &
2988 !$OMP SHARED(grid_atom,fxc,d1e,d2e,dft_control,na,nr,rhoa,rhob), &
2989 !$OMP PRIVATE(ia,ir)
2990  DO ir = 1, nr
2991  DO ia = 1, na
2992 
2993  !Need to be careful not to divide by zero. Assume that if rhoa == rhob, then
2994  !take the limit fxc = 0.5* (f_aa + f_bb - 2f_ab)
2995  IF (abs(rhoa(ia, ir, 1) - rhob(ia, ir, 1)) > dft_control%qs_control%eps_rho_rspace) THEN
2996  fxc(ia, ir) = grid_atom%weight(ia, ir)/(rhoa(ia, ir, 1) - rhob(ia, ir, 1)) &
2997  *(d1e(1)%array(ia, ir, 1) - d1e(2)%array(ia, ir, 1))
2998  ELSE
2999  fxc(ia, ir) = 0.5_dp*grid_atom%weight(ia, ir)* &
3000  (d2e(1)%array(ia, ir, 1) + d2e(3)%array(ia, ir, 1) - 2._dp*d2e(2)%array(ia, ir, 1))
3001  END IF
3002 
3003  END DO
3004  END DO
3005 !$OMP END PARALLEL DO
3006 
3007  ! Integrate wrt to so
3008  ALLOCATE (int_so(nsotot, nsotot))
3009  int_so = 0.0_dp
3010  CALL integrate_so_prod(int_so, fxc, ikind, xas_atom_env, qs_env)
3011 
3012  ! Contract into sgf. Array located on current processor only
3013  ALLOCATE (int_fxc(iatom, 4)%array(ri_nsgf, ri_nsgf))
3014  int_fxc(iatom, 4)%array = 0.0_dp
3015  CALL contract_so2sgf(int_fxc(iatom, 4)%array, int_so, ri_basis, ri_sphi_so)
3016 
3017  END SUBROUTINE integrate_sf_fxc
3018 
3019 ! **************************************************************************************************
3020 !> \brief Contract spherical orbitals to spherical Gaussians (so to sgf)
3021 !> \param int_sgf the array with the sgf integrals
3022 !> \param int_so the array with the so integrals (to contract)
3023 !> \param basis the corresponding gto basis set
3024 !> \param sphi_so the contraction coefficients for the s:
3025 ! **************************************************************************************************
3026  SUBROUTINE contract_so2sgf(int_sgf, int_so, basis, sphi_so)
3027 
3028  REAL(dp), DIMENSION(:, :) :: int_sgf, int_so
3029  TYPE(gto_basis_set_type), POINTER :: basis
3030  REAL(dp), DIMENSION(:, :) :: sphi_so
3031 
3032  INTEGER :: iset, jset, maxso, nset, nsoi, nsoj, &
3033  sgfi, sgfj, starti, startj
3034  INTEGER, DIMENSION(:), POINTER :: lmax, npgf, nsgf_set
3035  INTEGER, DIMENSION(:, :), POINTER :: first_sgf
3036 
3037  NULLIFY (nsgf_set, npgf, lmax, first_sgf)
3038 
3039  CALL get_gto_basis_set(basis, nset=nset, maxso=maxso, nsgf_set=nsgf_set, first_sgf=first_sgf, &
3040  npgf=npgf, lmax=lmax)
3041 
3042  DO iset = 1, nset
3043  starti = (iset - 1)*maxso + 1
3044  nsoi = npgf(iset)*nsoset(lmax(iset))
3045  sgfi = first_sgf(1, iset)
3046 
3047  DO jset = 1, nset
3048  startj = (jset - 1)*maxso + 1
3049  nsoj = npgf(jset)*nsoset(lmax(jset))
3050  sgfj = first_sgf(1, jset)
3051 
3052  CALL ab_contract(int_sgf(sgfi:sgfi + nsgf_set(iset) - 1, sgfj:sgfj + nsgf_set(jset) - 1), &
3053  int_so(starti:starti + nsoi - 1, startj:startj + nsoj - 1), &
3054  sphi_so(:, sgfi:), sphi_so(:, sgfj:), nsoi, nsoj, &
3055  nsgf_set(iset), nsgf_set(jset))
3056  END DO !jset
3057  END DO !iset
3058 
3059  END SUBROUTINE contract_so2sgf
3060 
3061 ! **************************************************************************************************
3062 !> \brief Integrate the product of spherical gaussian orbitals with the xc kernel on the atomic grid
3063 !> \param intso the integral in terms of spherical orbitals
3064 !> \param fxc the xc kernel at each grid point
3065 !> \param ikind the kind of the atom we integrate for
3066 !> \param xas_atom_env ...
3067 !> \param qs_env ...
3068 !> \note Largely copied from gaVxcgb_noGC. Rewritten here because we need our own atomic grid,
3069 !> harmonics, basis set and we do not need the soft vxc. Could have tweaked the original, but
3070 !> it would have been messy. Also we do not need rho_atom (too big and fancy for us)
3071 !> We also go over the whole range of angular momentum l
3072 ! **************************************************************************************************
3073  SUBROUTINE integrate_so_prod(intso, fxc, ikind, xas_atom_env, qs_env)
3074 
3075  REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: intso
3076  REAL(dp), DIMENSION(:, :), INTENT(IN) :: fxc
3077  INTEGER, INTENT(IN) :: ikind
3078  TYPE(xas_atom_env_type), POINTER :: xas_atom_env
3079  TYPE(qs_environment_type), POINTER :: qs_env
3080 
3081  CHARACTER(len=*), PARAMETER :: routinen = 'integrate_so_prod'
3082 
3083  INTEGER :: handle, ia, ic, icg, ipgf1, ipgf2, iset1, iset2, iso, iso1, iso2, l, ld, lmax12, &
3084  lmin12, m1, m2, max_iso_not0, max_iso_not0_local, max_s_harm, maxl, maxso, n1, n2, na, &
3085  ngau1, ngau2, nngau1, nr, nset, size1
3086  INTEGER, ALLOCATABLE, DIMENSION(:) :: cg_n_list
3087  INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cg_list
3088  INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf
3089  REAL(dp), ALLOCATABLE, DIMENSION(:) :: g1, g2
3090  REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: gfxcg, gg, matso
3091  REAL(dp), DIMENSION(:, :), POINTER :: zet
3092  REAL(dp), DIMENSION(:, :, :), POINTER :: my_cg
3093  TYPE(grid_atom_type), POINTER :: grid_atom
3094  TYPE(gto_basis_set_type), POINTER :: ri_basis
3095  TYPE(harmonics_atom_type), POINTER :: harmonics
3096  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
3097 
3098  CALL timeset(routinen, handle)
3099 
3100  NULLIFY (grid_atom, harmonics, ri_basis, qs_kind_set, lmax, lmin, npgf, zet, my_cg)
3101 
3102 ! Initialization
3103  CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set)
3104  CALL get_qs_kind(qs_kind_set(ikind), basis_set=ri_basis, basis_type="RI_XAS")
3105  grid_atom => xas_atom_env%grid_atom_set(ikind)%grid_atom
3106  harmonics => xas_atom_env%harmonics_atom_set(ikind)%harmonics_atom
3107 
3108  CALL get_gto_basis_set(ri_basis, lmax=lmax, lmin=lmin, maxso=maxso, maxl=maxl, npgf=npgf, &
3109  nset=nset, zet=zet)
3110 
3111  na = grid_atom%ng_sphere
3112  nr = grid_atom%nr
3113  my_cg => harmonics%my_CG
3114  max_iso_not0 = harmonics%max_iso_not0
3115  max_s_harm = harmonics%max_s_harm
3116  cpassert(2*maxl .LE. indso(1, max_iso_not0))
3117 
3118  ALLOCATE (g1(nr), g2(nr), gg(nr, 0:2*maxl))
3119  ALLOCATE (gfxcg(na, 0:2*maxl))
3120  ALLOCATE (matso(nsoset(maxl), nsoset(maxl)))
3121  ALLOCATE (cg_list(2, nsoset(maxl)**2, max_s_harm), cg_n_list(max_s_harm))
3122 
3123  g1 = 0.0_dp
3124  g2 = 0.0_dp
3125  m1 = 0
3126 ! Loop over the product of so
3127  DO iset1 = 1, nset
3128  n1 = nsoset(lmax(iset1))
3129  m2 = 0
3130  DO iset2 = 1, nset
3131  CALL get_none0_cg_list(my_cg, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
3132  max_s_harm, lmax(iset1) + lmax(iset2), cg_list, cg_n_list, &
3133  max_iso_not0_local)
3134  cpassert(max_iso_not0_local .LE. max_iso_not0)
3135 
3136  n2 = nsoset(lmax(iset2))
3137  DO ipgf1 = 1, npgf(iset1)
3138  ngau1 = n1*(ipgf1 - 1) + m1
3139  size1 = nsoset(lmax(iset1)) - nsoset(lmin(iset1) - 1)
3140  nngau1 = nsoset(lmin(iset1) - 1) + ngau1
3141 
3142  g1(:) = exp(-zet(ipgf1, iset1)*grid_atom%rad2(1:nr))
3143  DO ipgf2 = 1, npgf(iset2)
3144  ngau2 = n2*(ipgf2 - 1) + m2
3145 
3146  g2(:) = exp(-zet(ipgf2, iset2)*grid_atom%rad2(1:nr))
3147  lmin12 = lmin(iset1) + lmin(iset2)
3148  lmax12 = lmax(iset1) + lmax(iset2)
3149 
3150  !get the gaussian product
3151  gg = 0.0_dp
3152  IF (lmin12 == 0) THEN
3153  gg(:, lmin12) = g1(:)*g2(:)
3154  ELSE
3155  gg(:, lmin12) = grid_atom%rad2l(1:nr, lmin12)*g1(:)*g2(:)
3156  END IF
3157 
3158  DO l = lmin12 + 1, lmax12
3159  gg(:, l) = grid_atom%rad(1:nr)*gg(:, l - 1)
3160  END DO
3161 
3162  ld = lmax12 + 1
3163  CALL dgemm('N', 'N', na, ld, nr, 1.0_dp, fxc(1:na, 1:nr), na, gg(:, 0:lmax12), &
3164  nr, 0.0_dp, gfxcg(1:na, 0:lmax12), na)
3165 
3166  !integrate
3167  matso = 0.0_dp
3168  DO iso = 1, max_iso_not0_local
3169  DO icg = 1, cg_n_list(iso)
3170  iso1 = cg_list(1, icg, iso)
3171  iso2 = cg_list(2, icg, iso)
3172  l = indso(1, iso1) + indso(1, iso2)
3173 
3174  DO ia = 1, na
3175  matso(iso1, iso2) = matso(iso1, iso2) + gfxcg(ia, l)* &
3176  my_cg(iso1, iso2, iso)*harmonics%slm(ia, iso)
3177  END DO !ia
3178  END DO !icg
3179  END DO !iso
3180 
3181  !write in integral matrix
3182  DO ic = nsoset(lmin(iset2) - 1) + 1, nsoset(lmax(iset2))
3183  iso1 = nsoset(lmin(iset1) - 1) + 1
3184  iso2 = ngau2 + ic
3185  CALL daxpy(size1, 1.0_dp, matso(iso1:, ic), 1, intso(nngau1 + 1:, iso2), 1)
3186  END DO !ic
3187 
3188  END DO !ipgf2
3189  END DO ! ipgf1
3190  m2 = m2 + maxso
3191  END DO !iset2
3192  m1 = m1 + maxso
3193  END DO !iset1
3194 
3195  CALL timestop(handle)
3196 
3197  END SUBROUTINE integrate_so_prod
3198 
3199 ! **************************************************************************************************
3200 !> \brief This routine computes the integral of a potential V wrt the derivative of the spherical
3201 !> orbitals, that is <df/dx|V|dg/dy> on the atomic grid.
3202 !> \param intso the integral in terms of the spherical orbitals (well, their derivative)
3203 !> \param V the potential (put on the grid and weighted) to integrate
3204 !> \param ikind the atomic kind for which we integrate
3205 !> \param qs_env ...
3206 !> \param soc_atom_env ...
3207 ! **************************************************************************************************
3208  SUBROUTINE integrate_so_dxdy_prod(intso, V, ikind, qs_env, soc_atom_env)
3209 
3210  REAL(dp), DIMENSION(:, :, :), INTENT(INOUT) :: intso
3211  REAL(dp), DIMENSION(:, :), INTENT(IN) :: v
3212  INTEGER, INTENT(IN) :: ikind
3213  TYPE(qs_environment_type), POINTER :: qs_env
3214  TYPE(soc_atom_env_type), OPTIONAL, POINTER :: soc_atom_env
3215 
3216  INTEGER :: i, ipgf, iset, iso, j, jpgf, jset, jso, &
3217  k, l, maxso, na, nr, nset, starti, &
3218  startj
3219  INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf
3220  REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: fga, fgr, r1, r2, work
3221  REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: a1, a2
3222  REAL(dp), DIMENSION(:, :), POINTER :: slm, zet
3223  REAL(dp), DIMENSION(:, :, :), POINTER :: dslm_dxyz
3224  TYPE(grid_atom_type), POINTER :: grid_atom
3225  TYPE(gto_basis_set_type), POINTER :: basis
3226  TYPE(harmonics_atom_type), POINTER :: harmonics
3227  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
3228 
3229  NULLIFY (grid_atom, harmonics, basis, qs_kind_set, dslm_dxyz, slm, lmin, lmax, npgf, zet)
3230 
3231  CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set)
3232  CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis, basis_type="ORB")
3233  IF (PRESENT(soc_atom_env)) THEN
3234  grid_atom => soc_atom_env%grid_atom_set(ikind)%grid_atom
3235  harmonics => soc_atom_env%harmonics_atom_set(ikind)%harmonics_atom
3236  ELSE
3237  CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis, basis_type="ORB", harmonics=harmonics, &
3238  grid_atom=grid_atom)
3239  END IF
3240 
3241  na = grid_atom%ng_sphere
3242  nr = grid_atom%nr
3243 
3244  slm => harmonics%slm
3245  dslm_dxyz => harmonics%dslm_dxyz
3246 
3247 ! Getting what we need from the orbital basis
3248  CALL get_gto_basis_set(gto_basis_set=basis, lmax=lmax, lmin=lmin, &
3249  maxso=maxso, npgf=npgf, nset=nset, zet=zet)
3250 
3251 ! Separate the functions into purely r and purely angular parts, compute them all
3252 ! and use matrix mutliplication for the integral. We use f for x derivative and g for y
3253 
3254  ! Separating the functions. Note that the radial part is the same for x and y derivatives
3255  ALLOCATE (a1(na, nset*maxso, 3), a2(na, nset*maxso, 3))
3256  ALLOCATE (r1(nr, nset*maxso), r2(nr, nset*maxso))
3257  a1 = 0.0_dp; a2 = 0.0_dp
3258  r1 = 0.0_dp; r2 = 0.0_dp
3259 
3260  DO iset = 1, nset
3261  DO ipgf = 1, npgf(iset)
3262  starti = (iset - 1)*maxso + (ipgf - 1)*nsoset(lmax(iset))
3263  DO iso = nsoset(lmin(iset) - 1) + 1, nsoset(lmax(iset))
3264  l = indso(1, iso)
3265 
3266  ! The x derivative of the spherical orbital, divided in angular and radial parts
3267  ! 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)
3268 
3269  ! the purely radial part of d/dx(r^l Y_lm) * exp(-al*r^2) (same for y)
3270  r1(1:nr, starti + iso) = grid_atom%rad(1:nr)**(l - 1)*exp(-zet(ipgf, iset)*grid_atom%rad2(1:nr))
3271 
3272  ! the purely radial part of r^l Y_lm * d/dx(exp-al*r^2) (same for y)
3273  r2(1:nr, starti + iso) = -2.0_dp*zet(ipgf, iset)*grid_atom%rad(1:nr)**(l + 1) &
3274  *exp(-zet(ipgf, iset)*grid_atom%rad2(1:nr))
3275 
3276  DO i = 1, 3
3277  ! the purely angular part of d/dx(r^l Y_lm) * exp(-al*r^2)
3278  a1(1:na, starti + iso, i) = dslm_dxyz(i, 1:na, iso)
3279 
3280  ! the purely angular part of r^l Y_lm * d/dx(exp-al*r^2)
3281  a2(1:na, starti + iso, i) = harmonics%a(i, 1:na)*slm(1:na, iso)
3282  END DO
3283 
3284  END DO !iso
3285  END DO !ipgf
3286  END DO !iset
3287 
3288  ! Do the integration in terms of so using matrix products
3289  intso = 0.0_dp
3290  ALLOCATE (fga(na, 1))
3291  ALLOCATE (fgr(nr, 1))
3292  ALLOCATE (work(na, 1))
3293  fga = 0.0_dp; fgr = 0.0_dp; work = 0.0_dp
3294 
3295  DO iset = 1, nset
3296  DO jset = 1, nset
3297  DO ipgf = 1, npgf(iset)
3298  starti = (iset - 1)*maxso + (ipgf - 1)*nsoset(lmax(iset))
3299  DO jpgf = 1, npgf(jset)
3300  startj = (jset - 1)*maxso + (jpgf - 1)*nsoset(lmax(jset))
3301 
3302  DO i = 1, 3
3303  j = mod(i, 3) + 1
3304  k = mod(i + 1, 3) + 1
3305 
3306  DO iso = nsoset(lmin(iset) - 1) + 1, nsoset(lmax(iset))
3307  DO jso = nsoset(lmin(jset) - 1) + 1, nsoset(lmax(jset))
3308 
3309  !Two component per function => 4 terms in total
3310 
3311  ! take r1*a1(j) * V * r1*a1(k)
3312  fgr(1:nr, 1) = r1(1:nr, starti + iso)*r1(1:nr, startj + jso)
3313  fga(1:na, 1) = a1(1:na, starti + iso, j)*a1(1:na, startj + jso, k)
3314 
3315  CALL dgemm('N', 'N', na, 1, nr, 1.0_dp, v, na, fgr, nr, 0.0_dp, work, na)
3316  CALL dgemm('T', 'N', 1, 1, na, 1.0_dp, work, na, fga, na, 0.0_dp, &
3317  intso(starti + iso:, startj + jso, i), 1)
3318 
3319  ! add r1*a1(j) * V * r2*a2(k)
3320  fgr(1:nr, 1) = r1(1:nr, starti + iso)*r2(1:nr, startj + jso)
3321  fga(1:na, 1) = a1(1:na, starti + iso, j)*a2(1:na, startj + jso, k)
3322 
3323  CALL dgemm('N', 'N', na, 1, nr, 1.0_dp, v, na, fgr, nr, 0.0_dp, work, na)
3324  CALL dgemm('T', 'N', 1, 1, na, 1.0_dp, work, na, fga, na, 1.0_dp, &
3325  intso(starti + iso:, startj + jso, i), 1)
3326 
3327  ! add r2*a2(j) * V * r1*a1(k)
3328  fgr(1:nr, 1) = r2(1:nr, starti + iso)*r1(1:nr, startj + jso)
3329  fga(1:na, 1) = a2(1:na, starti + iso, j)*a1(1:na, startj + jso, k)
3330 
3331  CALL dgemm('N', 'N', na, 1, nr, 1.0_dp, v, na, fgr, nr, 0.0_dp, work, na)
3332  CALL dgemm('T', 'N', 1, 1, na, 1.0_dp, work, na, fga, na, 1.0_dp, &
3333  intso(starti + iso:, startj + jso, i), 1)
3334 
3335  ! add the last term: r2*a2(j) * V * r2*a2(k)
3336  fgr(1:nr, 1) = r2(1:nr, starti + iso)*r2(1:nr, startj + jso)
3337  fga(1:na, 1) = a2(1:na, starti + iso, j)*a2(1:na, startj + jso, k)
3338 
3339  CALL dgemm('N', 'N', na, 1, nr, 1.0_dp, v, na, fgr, nr, 0.0_dp, work, na)
3340  CALL dgemm('T', 'N', 1, 1, na, 1.0_dp, work, na, fga, na, 1.0_dp, &
3341  intso(starti + iso:, startj + jso, i), 1)
3342 
3343  END DO !jso
3344  END DO !iso
3345 
3346  END DO !i
3347  END DO !jpgf
3348  END DO !ipgf
3349  END DO !jset
3350  END DO !iset
3351 
3352  DO i = 1, 3
3353  intso(:, :, i) = intso(:, :, i) - transpose(intso(:, :, i))
3354  END DO
3355 
3356  END SUBROUTINE integrate_so_dxdy_prod
3357 
3358 ! **************************************************************************************************
3359 !> \brief Computes the SOC matrix elements with respect to the ORB basis set for each atomic kind
3360 !> and put them as the block diagonal of dbcsr_matrix
3361 !> \param matrix_soc the matrix where the SOC is stored
3362 !> \param xas_atom_env ...
3363 !> \param qs_env ...
3364 !> \param soc_atom_env ...
3365 !> \note We compute: <da_dx|V\‍(4c^2-2V)|db_dy> - <da_dy|V\‍(4c^2-2V)|db_dx>, where V is a model
3366 !> potential on the atomic grid. What we get is purely imaginary
3367 ! **************************************************************************************************
3368  SUBROUTINE integrate_soc_atoms(matrix_soc, xas_atom_env, qs_env, soc_atom_env)
3369 
3370  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_soc
3371  TYPE(xas_atom_env_type), OPTIONAL, POINTER :: xas_atom_env
3372  TYPE(qs_environment_type), POINTER :: qs_env
3373  TYPE(soc_atom_env_type), OPTIONAL, POINTER :: soc_atom_env
3374 
3375  CHARACTER(len=*), PARAMETER :: routinen = 'integrate_soc_atoms'
3376 
3377  INTEGER :: blk, handle, i, iat, ikind, ir, jat, &
3378  maxso, na, nkind, nr, nset, nsgf
3379  LOGICAL :: all_potential_present
3380  REAL(dp) :: zeff
3381  REAL(dp), ALLOCATABLE, DIMENSION(:) :: vr
3382  REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: v
3383  REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: intso
3384  REAL(dp), DIMENSION(:, :), POINTER :: sphi_so
3385  REAL(dp), DIMENSION(:, :, :), POINTER :: intsgf
3386  TYPE(cp_3d_r_p_type), DIMENSION(:), POINTER :: int_soc
3387  TYPE(dbcsr_iterator_type) :: iter
3388  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
3389  TYPE(grid_atom_type), POINTER :: grid
3390  TYPE(gto_basis_set_type), POINTER :: basis
3391  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
3392  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
3393 
3394 !!DEBUG
3395  CALL timeset(routinen, handle)
3396 
3397  NULLIFY (int_soc, basis, qs_kind_set, sphi_so, matrix_s, grid)
3398  NULLIFY (particle_set)
3399 
3400  ! Initialization
3401  CALL get_qs_env(qs_env, nkind=nkind, qs_kind_set=qs_kind_set, matrix_s=matrix_s, &
3402  particle_set=particle_set)
3403 
3404  ! all_potential_present
3405  CALL get_qs_kind_set(qs_kind_set, all_potential_present=all_potential_present)
3406 
3407  ! Loop over the kinds to compute the integrals
3408  ALLOCATE (int_soc(nkind))
3409  DO ikind = 1, nkind
3410  CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis, basis_type="ORB", zeff=zeff)
3411  IF (PRESENT(soc_atom_env)) THEN
3412  grid => soc_atom_env%grid_atom_set(ikind)%grid_atom
3413  ELSE
3414  CALL get_qs_kind(qs_kind_set(ikind), grid_atom=grid)
3415  END IF
3416  CALL get_gto_basis_set(basis, nset=nset, maxso=maxso)
3417  ALLOCATE (intso(nset*maxso, nset*maxso, 3))
3418 
3419  ! compute the model potential on the grid
3420  nr = grid%nr
3421  na = grid%ng_sphere
3422 
3423  ALLOCATE (vr(nr))
3424  CALL calculate_model_potential(vr, grid, zeff)
3425 
3426  ! Compute V/(4c^2-2V) and weight it
3427  ALLOCATE (v(na, nr))
3428  v = 0.0_dp
3429  DO ir = 1, nr
3430  CALL daxpy(na, vr(ir)/(4.0_dp*c_light_au**2 - 2.0_dp*vr(ir)), grid%weight(:, ir), 1, &
3431  v(:, ir), 1)
3432  END DO
3433  DEALLOCATE (vr)
3434 
3435  ! compute the integral <da_dx|...|db_dy> in terms of so
3436  IF (PRESENT(xas_atom_env)) THEN
3437  CALL integrate_so_dxdy_prod(intso, v, ikind, qs_env)
3438  ELSE
3439  CALL integrate_so_dxdy_prod(intso, v, ikind, qs_env, soc_atom_env)
3440  END IF
3441  DEALLOCATE (v)
3442 
3443  ! contract in terms of sgf
3444  CALL get_gto_basis_set(basis, nsgf=nsgf)
3445  ALLOCATE (int_soc(ikind)%array(nsgf, nsgf, 3))
3446  intsgf => int_soc(ikind)%array
3447  IF (PRESENT(xas_atom_env)) THEN
3448  sphi_so => xas_atom_env%orb_sphi_so(ikind)%array
3449  ELSE
3450  sphi_so => soc_atom_env%orb_sphi_so(ikind)%array
3451  END IF
3452  intsgf = 0.0_dp
3453 
3454  DO i = 1, 3
3455  CALL contract_so2sgf(intsgf(:, :, i), intso(:, :, i), basis, sphi_so)
3456  END DO
3457 
3458  DEALLOCATE (intso)
3459  END DO !ikind
3460 
3461  ! Build the matrix_soc based on the matrix_s (but anti-symmetric)
3462  IF ((PRESENT(xas_atom_env)) .OR. all_potential_present) THEN
3463  DO i = 1, 3
3464  CALL dbcsr_create(matrix_soc(i)%matrix, name="SOC MATRIX", template=matrix_s(1)%matrix, &
3465  matrix_type=dbcsr_type_antisymmetric)
3466  END DO
3467  ! Iterate over its diagonal blocks and fill=it
3468  CALL dbcsr_iterator_start(iter, matrix_s(1)%matrix)
3469  DO WHILE (dbcsr_iterator_blocks_left(iter))
3470 
3471  CALL dbcsr_iterator_next_block(iter, row=iat, column=jat, blk=blk)
3472  IF (.NOT. iat == jat) cycle
3473  ikind = particle_set(iat)%atomic_kind%kind_number
3474 
3475  DO i = 1, 3
3476  CALL dbcsr_put_block(matrix_soc(i)%matrix, iat, iat, int_soc(ikind)%array(:, :, i))
3477  END DO
3478 
3479  END DO !iat
3480  CALL dbcsr_iterator_stop(iter)
3481  ELSE ! pseudopotentials here
3482  DO i = 1, 3
3483  CALL dbcsr_create(matrix_soc(i)%matrix, name="SOC MATRIX", template=matrix_s(1)%matrix, &
3484  matrix_type=dbcsr_type_no_symmetry)
3485  CALL dbcsr_set(matrix_soc(i)%matrix, 0.0_dp)
3486  CALL dbcsr_copy(matrix_soc(i)%matrix, soc_atom_env%soc_pp(i, 1)%matrix)
3487  END DO
3488  END IF
3489 
3490  DO i = 1, 3
3491  CALL dbcsr_finalize(matrix_soc(i)%matrix)
3492  END DO
3493 
3494  ! Clean-up
3495  DO ikind = 1, nkind
3496  DEALLOCATE (int_soc(ikind)%array)
3497  END DO
3498  DEALLOCATE (int_soc)
3499 
3500  CALL timestop(handle)
3501 
3502  END SUBROUTINE integrate_soc_atoms
3503 
3504 END MODULE xas_tdp_atom
pbc
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
Definition: dumpdcd.F:1203
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition: grid_common.h:117
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
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition: cp_dbcsr_operations.F:89
ai_contraction_sphi
Contraction of integrals over primitive Cartesian Gaussians based on the contraction matrix sphi whic...
Definition: ai_contraction_sphi.F:15
ai_contraction_sphi::ab_contract
subroutine, public ab_contract(abint, sab, sphi_a, sphi_b, ncoa, ncob, nsgfa, nsgfb)
contract overlap integrals (a,b) and transfer to spherical Gaussians
Definition: ai_contraction_sphi.F:52
atom_operators
Calculate the atomic operator matrices.
Definition: atom_operators.F:15
atom_operators::calculate_model_potential
subroutine, public calculate_model_potential(modpot, grid, zcore)
Calculate model potential. It is useful to guess atomic orbitals.
Definition: atom_operators.F:1187
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
cell_types
Handles all functions related to the CELL.
Definition: cell_types.F:15
cp_array_utils
various utilities that regard array of different kinds: output, allocation,... maybe it is not a good...
Definition: cp_array_utils.F:18
cp_blacs_env
methods related to the blacs parallel environment
Definition: cp_blacs_env.F:15
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition: cp_control_types.F:12
cp_dbcsr_cholesky
Interface to (sca)lapack for the Cholesky based procedures.
Definition: cp_dbcsr_cholesky.F:17
cp_dbcsr_cholesky::cp_dbcsr_cholesky_decompose
subroutine, public cp_dbcsr_cholesky_decompose(matrix, n, para_env, blacs_env)
used to replace a symmetric positive def. matrix M with its cholesky decomposition U: M = U^T * U,...
Definition: cp_dbcsr_cholesky.F:61
cp_dbcsr_cholesky::cp_dbcsr_cholesky_invert
subroutine, public cp_dbcsr_cholesky_invert(matrix, n, para_env, blacs_env, upper_to_full)
used to replace the cholesky decomposition by the inverse
Definition: cp_dbcsr_cholesky.F:114
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition: cp_dbcsr_operations.F:18
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition: cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_io_unit
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
Definition: cp_log_handling.F:515
dbt_api
This is the start of a dbt_api, all publically needed functions are exported here....
Definition: dbt_api.F:17
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition: input_constants.F:17
input_constants::do_potential_id
integer, parameter, public do_potential_id
Definition: input_constants.F:801
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
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition: kinds.F:57
lebedev
Generation of the spherical Lebedev grids. All Lebedev grids were generated with a precision of at le...
Definition: lebedev.F:57
lebedev::deallocate_lebedev_grids
subroutine, public deallocate_lebedev_grids()
...
Definition: lebedev.F:324
lebedev::lebedev_grid
type(oh_grid), dimension(nlg), target, public lebedev_grid
Definition: lebedev.F:85
lebedev::get_number_of_lebedev_grid
integer function, public get_number_of_lebedev_grid(l, n)
Get the number of the Lebedev grid, which has the requested angular momentum quantnum number l or siz...
Definition: lebedev.F:114
lebedev::init_lebedev_grids
subroutine, public init_lebedev_grids()
Load the coordinates and weights of the nonredundant Lebedev grid points.
Definition: lebedev.F:344
libint_2c_3c
2- and 3-center electron repulsion integral routines based on libint2 Currently available operators: ...
Definition: libint_2c_3c.F:14
mathlib
Collection of simple mathematical functions and subroutines.
Definition: mathlib.F:15
mathlib::invmat_symm
subroutine, public invmat_symm(a, cholesky_triangle)
returns inverse of real symmetric, positive definite matrix
Definition: mathlib.F:579
mathlib::get_diag
pure real(kind=dp) function, dimension(min(size(a, 1), size(a, 2))), public get_diag(a)
Return the diagonal elements of matrix a as a vector.
Definition: mathlib.F:497
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::nco
integer, dimension(:), allocatable, public nco
Definition: orbital_pointers.F:42
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
orbital_pointers::ncoset
integer, dimension(:), allocatable, public ncoset
Definition: orbital_pointers.F:42
orbital_pointers::indco
integer, dimension(:, :), allocatable, public indco
Definition: orbital_pointers.F:43
orbital_pointers::nso
integer, dimension(:), allocatable, public nso
Definition: orbital_pointers.F:42
orbital_transformation_matrices
Calculation of the spherical harmonics and the corresponding orbital transformation matrices.
Definition: orbital_transformation_matrices.F:17
orbital_transformation_matrices::orbtramat
type(orbtramat_type), dimension(:), pointer, public orbtramat
Definition: orbital_transformation_matrices.F:47
particle_methods
Define methods related to particle_type.
Definition: particle_methods.F:14
particle_methods::get_particle_set
subroutine, public get_particle_set(particle_set, qs_kind_set, first_sgf, last_sgf, nsgf, nmao, basis)
Get the components of a particle set.
Definition: particle_methods.F:98
particle_types
Define the data structure for the particle information.
Definition: particle_types.F:19
physcon
Definition of physical constants:
Definition: physcon.F:68
physcon::c_light_au
real(kind=dp), parameter, public c_light_au
Definition: physcon.F:90
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_grid_atom::allocate_grid_atom
subroutine, public allocate_grid_atom(grid_atom)
Initialize components of the grid_atom_type structure.
Definition: qs_grid_atom.F:73
qs_grid_atom::create_grid_atom
subroutine, public create_grid_atom(grid_atom, nr, na, llmax, ll, quadrature)
...
Definition: qs_grid_atom.F:183
qs_harmonics_atom
Definition: qs_harmonics_atom.F:8
qs_harmonics_atom::allocate_harmonics_atom
subroutine, public allocate_harmonics_atom(harmonics)
Allocate a spherical harmonics set for the atom grid.
Definition: qs_harmonics_atom.F:66
qs_harmonics_atom::get_maxl_cg
subroutine, public get_maxl_cg(harmonics, orb_basis, llmax, max_s_harm)
...
Definition: qs_harmonics_atom.F:403
qs_harmonics_atom::create_harmonics_atom
subroutine, public create_harmonics_atom(harmonics, my_CG, na, llmax, maxs, max_s_harm, ll, wa, azi, pol)
...
Definition: qs_harmonics_atom.F:158
qs_integral_utils
Some utility functions for the calculation of integrals.
Definition: qs_integral_utils.F:14
qs_integral_utils::basis_set_list_setup
subroutine, public basis_set_list_setup(basis_set_list, basis_type, qs_kind_set)
Set up an easy accessible list of the basis sets for all kinds.
Definition: qs_integral_utils.F:149
qs_kind_types
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
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_kind_types::get_qs_kind_set
subroutine, public get_qs_kind_set(qs_kind_set, all_potential_present, tnadd_potential_present, gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, basis_rcut, basis_type, total_zeff_corr)
Get attributes of an atomic kind set.
Definition: qs_kind_types.F:864
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition: qs_neighbor_list_types.F:17
qs_neighbor_list_types::release_neighbor_list_sets
subroutine, public release_neighbor_list_sets(nlists)
releases an array of neighbor_list_sets
Definition: qs_neighbor_list_types.F:1097
qs_overlap
Calculation of overlap matrix, its derivatives and forces.
Definition: qs_overlap.F:19
qs_overlap::build_overlap_matrix_simple
subroutine, public build_overlap_matrix_simple(ks_env, matrix_s, basis_set_list_a, basis_set_list_b, sab_nl)
Calculation of the overlap matrix over Cartesian Gaussian functions.
Definition: qs_overlap.F:558
qs_rho_types
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition: qs_rho_types.F:18
qs_rho_types::qs_rho_get
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition: qs_rho_types.F:229
qs_tddfpt2_soc_types
Definition: qs_tddfpt2_soc_types.F:7
spherical_harmonics
Calculate spherical harmonics.
Definition: spherical_harmonics.F:27
spherical_harmonics::clebsch_gordon_init
subroutine, public clebsch_gordon_init(l)
...
Definition: spherical_harmonics.F:238
spherical_harmonics::clebsch_gordon_deallocate
subroutine, public clebsch_gordon_deallocate()
...
Definition: spherical_harmonics.F:220
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
xas_tdp_atom
This module deals with all the integrals done on local atomic grids in xas_tdp. This is mostly used t...
Definition: xas_tdp_atom.F:15
xas_tdp_atom::calculate_density_coeffs
subroutine, public calculate_density_coeffs(xas_atom_env, qs_env)
Compute the coefficients to project the density on a partial RI_XAS basis.
Definition: xas_tdp_atom.F:860
xas_tdp_atom::compute_sphi_so
subroutine, public compute_sphi_so(ikind, basis_type, sphi_so, qs_env)
Computes the contraction coefficients to go from spherical orbitals to sgf for a given atomic kind an...
Definition: xas_tdp_atom.F:450
xas_tdp_atom::init_xas_atom_env
subroutine, public init_xas_atom_env(xas_atom_env, xas_tdp_env, xas_tdp_control, qs_env, ltddfpt)
Initializes a xas_atom_env type given the qs_enxas_atom_env, qs_envv.
Definition: xas_tdp_atom.F:159
xas_tdp_atom::integrate_soc_atoms
subroutine, public integrate_soc_atoms(matrix_soc, xas_atom_env, qs_env, soc_atom_env)
Computes the SOC matrix elements with respect to the ORB basis set for each atomic kind and put them ...
Definition: xas_tdp_atom.F:3369
xas_tdp_atom::truncate_radial_grid
subroutine, public truncate_radial_grid(grid_atom, max_radius)
Reduces the radial extension of an atomic grid such that it only covers a given radius.
Definition: xas_tdp_atom.F:403
xas_tdp_atom::init_xas_atom_grid_harmo
subroutine, public init_xas_atom_grid_harmo(xas_atom_env, grid_info, do_xc, qs_env)
Initializes the atomic grids and harmonics for the RI atomic calculations.
Definition: xas_tdp_atom.F:237
xas_tdp_atom::integrate_fxc_atoms
subroutine, public integrate_fxc_atoms(int_fxc, xas_atom_env, xas_tdp_control, qs_env)
Integrate the xc kernel as a function of r on the atomic grids for the RI_XAS basis.
Definition: xas_tdp_atom.F:1875
xas_tdp_integrals
3-center integrals machinery for the XAS_TDP method
Definition: xas_tdp_integrals.F:13
xas_tdp_integrals::build_xas_tdp_3c_nl
subroutine, public build_xas_tdp_3c_nl(ac_list, basis_a, basis_c, op_type, qs_env, excited_atoms, x_range, ext_dist2d)
Builds a neighbor lists set taylored for 3-center integral within XAS TDP, such that only excited ato...
Definition: xas_tdp_integrals.F:765
xas_tdp_integrals::build_xas_tdp_ovlp_nl
subroutine, public build_xas_tdp_ovlp_nl(ab_list, basis_a, basis_b, qs_env, excited_atoms, ext_dist2d)
Builds a neighbor lists set for overlaping 2-center S_ab, where b is restricted on a a given list of ...
Definition: xas_tdp_integrals.F:665
xas_tdp_integrals::create_pqx_tensor
subroutine, public create_pqx_tensor(pq_X, ab_nl, ac_nl, matrix_dist, blk_size_1, blk_size_2, blk_size_3, only_bc_same_center)
Prepares a tensor to hold 3-center integrals (pq|X), where p,q are distributed according to the given...
Definition: xas_tdp_integrals.F:128
xas_tdp_integrals::fill_pqx_tensor
subroutine, public fill_pqx_tensor(pq_X, ab_nl, ac_nl, basis_set_list_a, basis_set_list_b, basis_set_list_c, potential_parameter, qs_env, only_bc_same_center, eps_screen)
Fills the given 3 dimensional (pq|X) tensor with integrals.
Definition: xas_tdp_integrals.F:277
xas_tdp_types
Define XAS TDP control type and associated create, release, etc subroutines, as well as XAS TDP envir...
Definition: xas_tdp_types.F:13
xas_tdp_types::get_proc_batch_sizes
subroutine, public get_proc_batch_sizes(batch_size, nbatch, nex_atom, nprocs)
Uses heuristics to determine a good batching of the processros for fxc integration.
Definition: xas_tdp_types.F:1427
xas_tdp_types::release_batch_info
subroutine, public release_batch_info(batch_info)
Releases a batch_info type.
Definition: xas_tdp_types.F:1400
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_derivative_desc
Module with functions to handle derivative descriptors. derivative description are strings have the f...
Definition: xc_derivative_desc.F:26
xc_derivative_desc::deriv_norm_drho
integer, parameter, public deriv_norm_drho
Definition: xc_derivative_desc.F:35
xc_derivative_desc::deriv_norm_drhoa
integer, parameter, public deriv_norm_drhoa
Definition: xc_derivative_desc.F:35
xc_derivative_desc::deriv_rhob
integer, parameter, public deriv_rhob
Definition: xc_derivative_desc.F:35
xc_derivative_desc::deriv_rhoa
integer, parameter, public deriv_rhoa
Definition: xc_derivative_desc.F:35
xc_derivative_desc::deriv_norm_drhob
integer, parameter, public deriv_norm_drhob
Definition: xc_derivative_desc.F:35
xc_derivative_set_types
represent a group ofunctional derivatives
Definition: xc_derivative_set_types.F:14
xc_derivative_set_types::xc_dset_get_derivative
type(xc_derivative_type) function, pointer, public xc_dset_get_derivative(derivative_set, description, allocate_deriv)
returns the requested xc_derivative
Definition: xc_derivative_set_types.F:66
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_derivative_types
Provides types for the management of the xc-functionals and their derivatives.
Definition: xc_derivative_types.F:12
xc_derivative_types::xc_derivative_get
subroutine, public xc_derivative_get(deriv, split_desc, order, deriv_data, accept_null_data)
returns various information on the given derivative
Definition: xc_derivative_types.F:103
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_derivatives::xc_functionals_eval
subroutine, public xc_functionals_eval(functionals, lsd, rho_set, deriv_set, deriv_order)
...
Definition: xc_derivatives.F:565
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
xc
Exchange and Correlation functional calculations.
Definition: xc.F:17
xc::divide_by_norm_drho
subroutine, public divide_by_norm_drho(deriv_set, rho_set, lsd)
divides derivatives from deriv_set by norm_drho
Definition: xc.F:5412