!--------------------------------------------------------------------------------------------------! ! CP2K: A general program to perform molecular dynamics simulations ! ! Copyright (C) 2000 - 2020 CP2K developers group ! !--------------------------------------------------------------------------------------------------! ! ************************************************************************************************** !> \brief This module deals with all the integrals done on local atomic grids in xas_tdp. This is !> mostly used to compute the xc kernel matrix elements wrt two RI basis elements (centered !> on the same excited atom)
, where the kernel fxc is purely a function of the !> ground state density and r. This is also used to compute the SOC matrix elements in the !> orbital basis ! ************************************************************************************************** MODULE xas_tdp_atom USE ai_contraction_sphi, ONLY: ab_contract USE atom_operators, ONLY: calculate_model_potential USE basis_set_types, ONLY: get_gto_basis_set,& gto_basis_set_p_type,& gto_basis_set_type USE cell_types, ONLY: cell_type,& pbc USE cp_array_utils, ONLY: cp_1d_i_p_type,& cp_1d_r_p_type,& cp_2d_r_p_type,& cp_3d_r_p_type USE cp_blacs_env, ONLY: cp_blacs_env_type USE cp_control_types, ONLY: dft_control_type,& qs_control_type USE cp_dbcsr_cholesky, ONLY: cp_dbcsr_cholesky_decompose,& cp_dbcsr_cholesky_invert USE cp_dbcsr_operations, ONLY: cp_dbcsr_dist2d_to_dist,& dbcsr_deallocate_matrix_set USE cp_log_handling, ONLY: cp_logger_get_default_io_unit USE cp_para_env, ONLY: cp_para_env_create USE cp_para_types, ONLY: cp_para_env_type USE dbcsr_api, ONLY: & dbcsr_complete_redistribute, dbcsr_create, dbcsr_distribution_get, dbcsr_distribution_new, & dbcsr_distribution_release, dbcsr_distribution_type, dbcsr_filter, dbcsr_finalize, & dbcsr_get_block_p, dbcsr_get_stored_coordinates, dbcsr_iterator_blocks_left, & dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, & dbcsr_p_type, dbcsr_put_block, dbcsr_release, dbcsr_replicate_all, dbcsr_type USE dbcsr_tensor_api, ONLY: dbcsr_t_destroy,& dbcsr_t_get_block,& dbcsr_t_iterator_blocks_left,& dbcsr_t_iterator_next_block,& dbcsr_t_iterator_start,& dbcsr_t_iterator_stop,& dbcsr_t_iterator_type,& dbcsr_t_type USE distribution_2d_types, ONLY: distribution_2d_release,& distribution_2d_type USE input_constants, ONLY: do_potential_id USE input_section_types, ONLY: section_vals_get_subs_vals,& section_vals_type USE kinds, ONLY: default_string_length,& dp USE lebedev, ONLY: deallocate_lebedev_grids,& get_number_of_lebedev_grid,& init_lebedev_grids,& lebedev_grid USE libint_2c_3c, ONLY: libint_potential_type USE mathconstants, ONLY: dfac,& pi USE mathlib, ONLY: get_diag USE memory_utilities, ONLY: reallocate USE message_passing, ONLY: mp_comm_split_direct,& mp_irecv,& mp_isend,& mp_sum,& mp_sync,& mp_waitall USE orbital_pointers, ONLY: indco,& indso,& nco,& ncoset,& nso,& nsoset USE orbital_transformation_matrices, ONLY: orbtramat USE particle_methods, ONLY: get_particle_set USE particle_types, ONLY: particle_type USE physcon, ONLY: c_light_au USE qs_environment_types, ONLY: get_qs_env,& qs_environment_type USE qs_grid_atom, ONLY: allocate_grid_atom,& create_grid_atom,& grid_atom_type USE qs_harmonics_atom, ONLY: allocate_harmonics_atom,& create_harmonics_atom,& get_maxl_CG,& get_none0_cg_list,& harmonics_atom_type USE qs_integral_utils, ONLY: basis_set_list_setup USE qs_kind_types, ONLY: get_qs_kind,& get_qs_kind_set,& qs_kind_type USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type,& release_neighbor_list_sets USE qs_overlap, ONLY: build_overlap_matrix_simple USE qs_rho_types, ONLY: qs_rho_get,& qs_rho_type USE spherical_harmonics, ONLY: clebsch_gordon,& clebsch_gordon_deallocate,& clebsch_gordon_init USE util, ONLY: get_limit,& sort_unique USE xas_tdp_integrals, ONLY: build_xas_tdp_3c_nl,& build_xas_tdp_ovlp_nl,& create_pqX_tensor,& fill_pqx_tensor,& get_opt_3c_dist2d USE xas_tdp_types, ONLY: batch_info_type,& get_proc_batch_sizes,& release_batch_info,& xas_atom_env_type,& xas_tdp_control_type,& xas_tdp_env_type USE xc, ONLY: divide_by_norm_drho USE xc_atom, ONLY: xc_rho_set_atom_update USE xc_derivative_set_types, ONLY: xc_derivative_set_type,& xc_dset_create,& xc_dset_get_derivative,& xc_dset_release USE xc_derivative_types, ONLY: xc_derivative_get,& xc_derivative_type USE xc_derivatives, ONLY: xc_functionals_eval,& xc_functionals_get_needs USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type USE xc_rho_set_types, ONLY: xc_rho_set_create,& xc_rho_set_release,& xc_rho_set_type !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num #include "./base/base_uses.f90" IMPLICIT NONE PRIVATE CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xas_tdp_atom' PUBLIC :: init_xas_atom_env, integrate_fxc_atoms, integrate_soc_atoms CONTAINS ! ************************************************************************************************** !> \brief Initializes a xas_atom_env type given the qs_enxas_atom_env, qs_envv !> \param xas_atom_env the xas_atom_env to initialize !> \param xas_tdp_env ... !> \param xas_tdp_control ... !> \param qs_env ... ! ************************************************************************************************** SUBROUTINE init_xas_atom_env(xas_atom_env, xas_tdp_env, xas_tdp_control, qs_env) TYPE(xas_atom_env_type), POINTER :: xas_atom_env TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control TYPE(qs_environment_type), POINTER :: qs_env CHARACTER(len=*), PARAMETER :: routineN = 'init_xas_atom_env', & routineP = moduleN//':'//routineN INTEGER :: handle, ikind, natom, nex_atoms, & nex_kinds, nkind, nspins LOGICAL :: do_soc, do_xc TYPE(particle_type), DIMENSION(:), POINTER :: particle_set TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set NULLIFY (qs_kind_set, particle_set) CALL timeset(routineN, handle) ! Initializing the type CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, natom=natom, particle_set=particle_set) nkind = SIZE(qs_kind_set) nex_kinds = xas_tdp_env%nex_kinds nex_atoms = xas_tdp_env%nex_atoms do_xc = xas_tdp_control%do_xc do_soc = xas_tdp_control%do_soc nspins = 1; IF (xas_tdp_control%do_uks .OR. xas_tdp_control%do_roks) nspins = 2 xas_atom_env%nspins = nspins xas_atom_env%ri_radius = xas_tdp_control%ri_radius ALLOCATE (xas_atom_env%grid_atom_set(nkind)) ALLOCATE (xas_atom_env%harmonics_atom_set(nkind)) ALLOCATE (xas_atom_env%ri_sphi_so(nkind)) ALLOCATE (xas_atom_env%orb_sphi_so(nkind)) IF (do_xc) THEN ALLOCATE (xas_atom_env%gr(nkind)) ALLOCATE (xas_atom_env%ga(nkind)) ALLOCATE (xas_atom_env%dgr1(nkind)) ALLOCATE (xas_atom_env%dgr2(nkind)) ALLOCATE (xas_atom_env%dga1(nkind)) ALLOCATE (xas_atom_env%dga2(nkind)) END IF xas_atom_env%excited_atoms => xas_tdp_env%ex_atom_indices xas_atom_env%excited_kinds => xas_tdp_env%ex_kind_indices ! Allocate and initialize the atomic grids and harmonics CALL init_xas_atom_grid_harmo(xas_atom_env, xas_tdp_control%grid_info, do_xc, qs_env) ! Compute the contraction coefficients for spherical orbitals DO ikind = 1, nkind NULLIFY (xas_atom_env%orb_sphi_so(ikind)%array, xas_atom_env%ri_sphi_so(ikind)%array) IF (do_soc) THEN CALL compute_sphi_so(ikind, "ORB", xas_atom_env%orb_sphi_so(ikind)%array, qs_env) END IF IF (ANY(xas_atom_env%excited_kinds == ikind) .AND. do_xc) THEN CALL compute_sphi_so(ikind, "RI_XAS", xas_atom_env%ri_sphi_so(ikind)%array, qs_env) END IF END DO !ikind ! Compute the coefficients to expand the density in the RI_XAS basis, if requested IF (do_xc) CALL calculate_density_coeffs(xas_atom_env, qs_env) CALL timestop(handle) END SUBROUTINE init_xas_atom_env ! ************************************************************************************************** !> \brief Initializes the atomic grids and harmonics for the RI atomic calculations !> \param xas_atom_env ... !> \param grid_info ... !> \param do_xc Whether the xc kernel will ne computed on the atomic grids. If not, the harmonics !> are built for the orbital basis for all kinds. !> \param qs_env ... !> \note Largely inspired by init_rho_atom subroutine ! ************************************************************************************************** SUBROUTINE init_xas_atom_grid_harmo(xas_atom_env, grid_info, do_xc, qs_env) TYPE(xas_atom_env_type), POINTER :: xas_atom_env CHARACTER(len=default_string_length), & DIMENSION(:, :), POINTER :: grid_info LOGICAL, INTENT(IN) :: do_xc TYPE(qs_environment_type), POINTER :: qs_env CHARACTER(len=*), PARAMETER :: routineN = 'init_xas_atom_grid_harmo', & routineP = moduleN//':'//routineN CHARACTER(LEN=default_string_length) :: kind_name INTEGER :: igrid, ikind, il, iso, iso1, iso2, l1, l1l2, l2, la, lc1, lc2, lcleb, ll, llmax, & lp, m1, m2, max_s_harm, max_s_set, maxl, maxlgto, maxs, mm, mp, na, nr, quadrature, stat REAL(dp) :: kind_radius REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: rga REAL(dp), DIMENSION(:, :, :), POINTER :: my_CG TYPE(dft_control_type), POINTER :: dft_control TYPE(grid_atom_type), POINTER :: grid_atom TYPE(gto_basis_set_type), POINTER :: tmp_basis TYPE(harmonics_atom_type), POINTER :: harmonics TYPE(qs_control_type), POINTER :: qs_control TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set NULLIFY (my_CG, qs_kind_set, tmp_basis, grid_atom, harmonics, qs_control, dft_control) ! Initialization of some integer for the CG coeff generation CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, dft_control=dft_control) IF (do_xc) THEN CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto, basis_type="RI_XAS") ELSE CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto, basis_type="ORB") END IF qs_control => dft_control%qs_control !maximum expansion llmax = 2*maxlgto max_s_harm = nsoset(llmax) max_s_set = nsoset(maxlgto) lcleb = llmax ! Allocate and compute the CG coeffs (copied from init_rho_atom) CALL clebsch_gordon_init(lcleb) CALL reallocate(my_CG, 1, max_s_set, 1, max_s_set, 1, max_s_harm) ALLOCATE (rga(lcleb, 2)) DO lc1 = 0, maxlgto DO iso1 = nsoset(lc1 - 1) + 1, nsoset(lc1) l1 = indso(1, iso1) m1 = indso(2, iso1) DO lc2 = 0, maxlgto DO iso2 = nsoset(lc2 - 1) + 1, nsoset(lc2) l2 = indso(1, iso2) m2 = indso(2, iso2) CALL clebsch_gordon(l1, m1, l2, m2, rga) IF (l1 + l2 > llmax) THEN l1l2 = llmax ELSE l1l2 = l1 + l2 END IF mp = m1 + m2 mm = m1 - m2 IF (m1*m2 < 0 .OR. (m1*m2 == 0 .AND. (m1 < 0 .OR. m2 < 0))) THEN mp = -ABS(mp) mm = -ABS(mm) ELSE mp = ABS(mp) mm = ABS(mm) END IF DO lp = MOD(l1 + l2, 2), l1l2, 2 il = lp/2 + 1 IF (ABS(mp) <= lp) THEN IF (mp >= 0) THEN iso = nsoset(lp - 1) + lp + 1 + mp ELSE iso = nsoset(lp - 1) + lp + 1 - ABS(mp) END IF my_CG(iso1, iso2, iso) = rga(il, 1) ENDIF IF (mp /= mm .AND. ABS(mm) <= lp) THEN IF (mm >= 0) THEN iso = nsoset(lp - 1) + lp + 1 + mm ELSE iso = nsoset(lp - 1) + lp + 1 - ABS(mm) END IF my_CG(iso1, iso2, iso) = rga(il, 2) ENDIF END DO ENDDO ! iso2 ENDDO ! lc2 ENDDO ! iso1 ENDDO ! lc1 DEALLOCATE (rga) CALL clebsch_gordon_deallocate() ! Create the Lebedev grids and compute the spherical harmonics CALL init_lebedev_grids() quadrature = qs_control%gapw_control%quadrature DO ikind = 1, SIZE(xas_atom_env%grid_atom_set) ! Allocate the grid and the harmonics for this kind NULLIFY (xas_atom_env%grid_atom_set(ikind)%grid_atom) NULLIFY (xas_atom_env%harmonics_atom_set(ikind)%harmonics_atom) CALL allocate_grid_atom(xas_atom_env%grid_atom_set(ikind)%grid_atom) CALL allocate_harmonics_atom(xas_atom_env%harmonics_atom_set(ikind)%harmonics_atom) NULLIFY (grid_atom, harmonics) grid_atom => xas_atom_env%grid_atom_set(ikind)%grid_atom harmonics => xas_atom_env%harmonics_atom_set(ikind)%harmonics_atom ! Initialize some integers CALL get_qs_kind(qs_kind_set(ikind), ngrid_rad=nr, ngrid_ang=na, name=kind_name) !take the grid dimension given as input, if none, take the GAPW ones above IF (ANY(grid_info == kind_name)) THEN DO igrid = 1, SIZE(grid_info, 1) IF (grid_info(igrid, 1) == kind_name) THEN !hack to convert string into integer READ (grid_info(igrid, 2), *, iostat=stat) na IF (stat .NE. 0) CPABORT("The 'na' value for the GRID keyword must be an integer") READ (grid_info(igrid, 3), *, iostat=stat) nr IF (stat .NE. 0) CPABORT("The 'nr' value for the GRID keyword must be an integer") EXIT END IF END DO ELSE IF (do_xc .AND. ANY(xas_atom_env%excited_kinds == ikind)) THEN !need good integration grids for the xc kernel, but taking the default GAPW grid CPWARN("The default (and possibly small) GAPW grid is being used for one excited KIND") END IF ll = get_number_of_lebedev_grid(n=na) na = lebedev_grid(ll)%n la = lebedev_grid(ll)%l grid_atom%ng_sphere = na grid_atom%nr = nr ! If this is an excited kind, create the harmonics with the RI_XAS basis, otherwise the ORB IF (ANY(xas_atom_env%excited_kinds == ikind) .AND. do_xc) THEN CALL get_qs_kind(qs_kind_set(ikind), basis_set=tmp_basis, basis_type="RI_XAS") ELSE CALL get_qs_kind(qs_kind_set(ikind), basis_set=tmp_basis, basis_type="ORB") END IF CALL get_gto_basis_set(gto_basis_set=tmp_basis, maxl=maxl, kind_radius=kind_radius) CALL create_grid_atom(grid_atom, nr, na, llmax, ll, quadrature) CALL truncate_radial_grid(grid_atom, kind_radius) maxs = nsoset(maxl) CALL create_harmonics_atom(harmonics, & my_CG, na, llmax, maxs, max_s_harm, ll, grid_atom%wa, & grid_atom%azi, grid_atom%pol) CALL get_maxl_CG(harmonics, tmp_basis, llmax, max_s_harm) END DO CALL deallocate_lebedev_grids() DEALLOCATE (my_CG) END SUBROUTINE init_xas_atom_grid_harmo ! ************************************************************************************************** !> \brief Reduces the radial extension of an atomic grid such that it only covers a given radius !> \param grid_atom the atomic grid from which we truncate the radial part !> \param max_radius the maximal radial extension of the resulting grid !> \note Since the RI density used for
is only of quality close to the atom, and the
!> integrand only non-zero within the radius of the gaussian P,Q. One can reduce the grid
!> extansion to the largest radius of the RI basis set
! **************************************************************************************************
SUBROUTINE truncate_radial_grid(grid_atom, max_radius)
TYPE(grid_atom_type), POINTER :: grid_atom
REAL(dp), INTENT(IN) :: max_radius
CHARACTER(len=*), PARAMETER :: routineN = 'truncate_radial_grid', &
routineP = moduleN//':'//routineN
INTEGER :: first_ir, ir, llmax_p1, na, new_nr, nr
nr = grid_atom%nr
na = grid_atom%ng_sphere
llmax_p1 = SIZE(grid_atom%rad2l, 2) - 1
! Find the index corresponding to the limiting radius (small ir => large radius)
DO ir = 1, nr
IF (grid_atom%rad(ir) < max_radius) THEN
first_ir = ir
EXIT
END IF
END DO
new_nr = nr - first_ir + 1
! Reallcoate everything that depends on r
grid_atom%nr = new_nr
grid_atom%rad(1:new_nr) = grid_atom%rad(first_ir:nr)
grid_atom%rad2(1:new_nr) = grid_atom%rad2(first_ir:nr)
grid_atom%wr(1:new_nr) = grid_atom%wr(first_ir:nr)
grid_atom%weight(:, 1:new_nr) = grid_atom%weight(:, first_ir:nr)
grid_atom%rad2l(1:new_nr, :) = grid_atom%rad2l(first_ir:nr, :)
grid_atom%oorad2l(1:new_nr, :) = grid_atom%oorad2l(first_ir:nr, :)
CALL reallocate(grid_atom%rad, 1, new_nr)
CALL reallocate(grid_atom%rad2, 1, new_nr)
CALL reallocate(grid_atom%wr, 1, new_nr)
CALL reallocate(grid_atom%weight, 1, na, 1, new_nr)
CALL reallocate(grid_atom%rad2l, 1, new_nr, 0, llmax_p1)
CALL reallocate(grid_atom%oorad2l, 1, new_nr, 0, llmax_p1)
END SUBROUTINE truncate_radial_grid
! **************************************************************************************************
!> \brief Computes the contraction coefficients to go from spherical orbitals to sgf for a given
!> atomic kind and a given basis type.
!> \param ikind the kind for which we compute the coefficients
!> \param basis_type the type of basis for which we compute
!> \param sphi_so where the new contraction coefficients are stored (not yet allocated)
!> \param qs_env ...
! **************************************************************************************************
SUBROUTINE compute_sphi_so(ikind, basis_type, sphi_so, qs_env)
INTEGER, INTENT(IN) :: ikind
CHARACTER(len=*), INTENT(IN) :: basis_type
REAL(dp), DIMENSION(:, :), POINTER :: sphi_so
TYPE(qs_environment_type), POINTER :: qs_env
CHARACTER(len=*), PARAMETER :: routineN = 'compute_sphi_so', &
routineP = moduleN//':'//routineN
INTEGER :: ico, ipgf, iset, iso, l, lx, ly, lz, &
maxso, nset, sgfi, start_c, start_s
INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf, nsgf_set
INTEGER, DIMENSION(:, :), POINTER :: first_sgf
REAL(dp) :: factor
REAL(dp), DIMENSION(:, :), POINTER :: sphi
TYPE(gto_basis_set_type), POINTER :: basis
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
NULLIFY (basis, lmax, lmin, npgf, nsgf_set, qs_kind_set, first_sgf, sphi)
CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set)
CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis, basis_type=basis_type)
CALL get_gto_basis_set(basis, lmax=lmax, nset=nset, npgf=npgf, maxso=maxso, lmin=lmin, &
nsgf_set=nsgf_set, sphi=sphi, first_sgf=first_sgf)
ALLOCATE (sphi_so(maxso, SUM(nsgf_set)))
sphi_so = 0.0_dp
DO iset = 1, nset
sgfi = first_sgf(1, iset)
DO ipgf = 1, npgf(iset)
start_s = (ipgf - 1)*nsoset(lmax(iset))
start_c = (ipgf - 1)*ncoset(lmax(iset))
DO l = lmin(iset), lmax(iset)
DO iso = 1, nso(l)
DO ico = 1, nco(l)
lx = indco(1, ico + ncoset(l - 1))
ly = indco(2, ico + ncoset(l - 1))
lz = indco(3, ico + ncoset(l - 1))
factor = orbtramat(l)%s2c(iso, ico) &
*SQRT(4.0_dp*pi/dfac(2*l + 1)*dfac(2*lx - 1)*dfac(2*ly - 1)*dfac(2*lz - 1))
CALL daxpy(nsgf_set(iset), factor, &
sphi(start_c + ncoset(l - 1) + ico, sgfi:sgfi + nsgf_set(iset) - 1), 1, &
sphi_so(start_s + nsoset(l - 1) + iso, sgfi:sgfi + nsgf_set(iset) - 1), 1)
END DO !ico
END DO !iso
END DO !l
END DO !ipgf
END DO !iset
END SUBROUTINE compute_sphi_so
! **************************************************************************************************
!> \brief Find the neighbors of a given set of atoms based on the non-zero blocks of a provided
!> overlap matrix. Optionally returns an array containing the indices of all involved atoms
!> (the given subset plus all their neighbors, without repetition) AND/OR an array of arrays
!> providing the indices of the neighbors of each input atom.
!> \param base_atoms the set of atoms for which we search neighbors
!> \param mat_s the overlap matrix used to find neighbors
!> \param radius the cutoff radius after which atoms are not considered neighbors
!> \param qs_env ...
!> \param all_neighbors the array uniquely contatining all indices of all atoms involved
!> \param neighbor_set array of arrays containing the neighbors of all given atoms
! **************************************************************************************************
SUBROUTINE find_neighbors(base_atoms, mat_s, radius, qs_env, all_neighbors, neighbor_set)
INTEGER, DIMENSION(:), INTENT(INOUT) :: base_atoms
TYPE(dbcsr_type), INTENT(IN) :: mat_s
REAL(dp) :: radius
TYPE(qs_environment_type), POINTER :: qs_env
INTEGER, DIMENSION(:), OPTIONAL, POINTER :: all_neighbors
TYPE(cp_1d_i_p_type), DIMENSION(:), OPTIONAL, &
POINTER :: neighbor_set
CHARACTER(len=*), PARAMETER :: routineN = 'find_neighbors', routineP = moduleN//':'//routineN
INTEGER :: blk, i, iat, ibase, iblk, jblk, mepos, &
natom, nb, nbase
INTEGER, ALLOCATABLE, DIMENSION(:) :: blk_to_base, inb, who_is_there
INTEGER, ALLOCATABLE, DIMENSION(:, :) :: n_neighbors
LOGICAL, ALLOCATABLE, DIMENSION(:) :: is_base_atom
REAL(dp) :: dist2, rad2, ri(3), rij(3), rj(3)
TYPE(cell_type), POINTER :: cell
TYPE(cp_1d_i_p_type), DIMENSION(:), POINTER :: my_neighbor_set
TYPE(cp_para_env_type), POINTER :: para_env
TYPE(dbcsr_iterator_type) :: iter
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
NULLIFY (particle_set, para_env, my_neighbor_set, cell)
! Initialization
CALL get_qs_env(qs_env, para_env=para_env, natom=natom, particle_set=particle_set, cell=cell)
mepos = para_env%mepos
nbase = SIZE(base_atoms)
!work with the neighbor_set structure, see at the end if we keep it
ALLOCATE (my_neighbor_set(nbase))
rad2 = radius**2
ALLOCATE (blk_to_base(natom), is_base_atom(natom))
blk_to_base = 0; is_base_atom = .FALSE.
DO ibase = 1, nbase
blk_to_base(base_atoms(ibase)) = ibase
is_base_atom(base_atoms(ibase)) = .TRUE.
END DO
! First loop over S => count the number of neighbors
ALLOCATE (n_neighbors(nbase, 0:para_env%num_pe - 1))
n_neighbors = 0
CALL dbcsr_iterator_start(iter, mat_s)
DO WHILE (dbcsr_iterator_blocks_left(iter))
CALL dbcsr_iterator_next_block(iter, row=iblk, column=jblk, blk=blk)
!avoid self-neighbors
IF (iblk == jblk) CYCLE
!test distance
ri = pbc(particle_set(iblk)%r, cell)
rj = pbc(particle_set(jblk)%r, cell)
rij = pbc(ri, rj, cell)
dist2 = SUM(rij**2)
IF (dist2 > rad2) CYCLE
IF (is_base_atom(iblk)) THEN
ibase = blk_to_base(iblk)
n_neighbors(ibase, mepos) = n_neighbors(ibase, mepos) + 1
END IF
IF (is_base_atom(jblk)) THEN
ibase = blk_to_base(jblk)
n_neighbors(ibase, mepos) = n_neighbors(ibase, mepos) + 1
END IF
END DO !iter
CALL dbcsr_iterator_stop(iter)
CALL mp_sum(n_neighbors, para_env%group)
! Allocate the neighbor_set arrays at the correct length
DO ibase = 1, nbase
ALLOCATE (my_neighbor_set(ibase)%array(SUM(n_neighbors(ibase, :))))
my_neighbor_set(ibase)%array = 0
END DO
! Loop a second time over S, this time fill the neighbors details
CALL dbcsr_iterator_start(iter, mat_s)
ALLOCATE (inb(nbase))
inb = 1
DO WHILE (dbcsr_iterator_blocks_left(iter))
CALL dbcsr_iterator_next_block(iter, row=iblk, column=jblk, blk=blk)
IF (iblk == jblk) CYCLE
!test distance
ri = pbc(particle_set(iblk)%r, cell)
rj = pbc(particle_set(jblk)%r, cell)
rij = pbc(ri, rj, cell)
dist2 = SUM(rij**2)
IF (dist2 > rad2) CYCLE
IF (is_base_atom(iblk)) THEN
ibase = blk_to_base(iblk)
my_neighbor_set(ibase)%array(SUM(n_neighbors(ibase, 0:mepos - 1)) + inb(ibase)) = jblk
inb(ibase) = inb(ibase) + 1
END IF
IF (is_base_atom(jblk)) THEN
ibase = blk_to_base(jblk)
my_neighbor_set(ibase)%array(SUM(n_neighbors(ibase, 0:mepos - 1)) + inb(ibase)) = iblk
inb(ibase) = inb(ibase) + 1
END IF
END DO !iter
CALL dbcsr_iterator_stop(iter)
! Make sure the info is shared among the procs
DO ibase = 1, nbase
CALL mp_sum(my_neighbor_set(ibase)%array, para_env%group)
END DO
! Gather all indices if asked for it
IF (PRESENT(all_neighbors)) THEN
ALLOCATE (who_is_there(natom))
who_is_there = 0
DO ibase = 1, nbase
who_is_there(base_atoms(ibase)) = 1
DO nb = 1, SIZE(my_neighbor_set(ibase)%array)
who_is_there(my_neighbor_set(ibase)%array(nb)) = 1
END DO
END DO
ALLOCATE (all_neighbors(natom))
i = 0
DO iat = 1, natom
IF (who_is_there(iat) == 1) THEN
i = i + 1
all_neighbors(i) = iat
END IF
END DO
CALL reallocate(all_neighbors, 1, i)
END IF
! If not asked for the neighbor set, deallocate it
IF (PRESENT(neighbor_set)) THEN
neighbor_set => my_neighbor_set
ELSE
DO ibase = 1, nbase
DEALLOCATE (my_neighbor_set(ibase)%array)
END DO
DEALLOCATE (my_neighbor_set)
END IF
END SUBROUTINE find_neighbors
! **************************************************************************************************
!> \brief Returns the RI inverse overlap for a subset of the RI_XAS matrix spaning a given
!> excited atom and its neighbors.
!> \param ri_sinv the inverse overlap as a dbcsr matrix
!> \param whole_s the whole RI overlap matrix
!> \param neighbors the indeces of the excited atom and their neighbors
!> \param idx_to_nb array telling where any atom can be found in neighbors (if there at all)
!> \param basis_set_ri the RI basis set list for all kinds
!> \param qs_env ...
!> \note It is assumed that the neighbors are sorted, the output matrix is assumed to be small and
!> is replicated on all processors
! **************************************************************************************************
SUBROUTINE get_exat_ri_sinv(ri_sinv, whole_s, neighbors, idx_to_nb, basis_set_ri, qs_env)
TYPE(dbcsr_type) :: ri_sinv, whole_s
INTEGER, DIMENSION(:), INTENT(IN) :: neighbors, idx_to_nb
TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_ri
TYPE(qs_environment_type), POINTER :: qs_env
CHARACTER(len=*), PARAMETER :: routineN = 'get_exat_ri_sinv', &
routineP = moduleN//':'//routineN
INTEGER :: blk, dest, group, handle, iat, ikind, &
inb, ir, is, jat, jnb, natom, nnb, &
npcols, nprows, source, tag
INTEGER, ALLOCATABLE, DIMENSION(:) :: recv_req, send_req
INTEGER, DIMENSION(:), POINTER :: col_dist, nsgf, row_dist
INTEGER, DIMENSION(:, :), POINTER :: pgrid
LOGICAL :: found_risinv, found_whole
LOGICAL, ALLOCATABLE, DIMENSION(:) :: is_neighbor
REAL(dp) :: ri(3), rij(3), rj(3)
REAL(dp), ALLOCATABLE, DIMENSION(:) :: radius
REAL(dp), DIMENSION(:, :), POINTER :: block_risinv, block_whole
TYPE(cell_type), POINTER :: cell
TYPE(cp_2d_r_p_type), ALLOCATABLE, DIMENSION(:) :: recv_buff, send_buff
TYPE(cp_blacs_env_type), POINTER :: blacs_env
TYPE(cp_para_env_type), POINTER :: para_env
TYPE(dbcsr_distribution_type) :: sinv_dist
TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
TYPE(dbcsr_iterator_type) :: iter
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
NULLIFY (pgrid, dbcsr_dist, row_dist, col_dist, nsgf, particle_set, block_whole, block_risinv)
NULLIFY (cell, para_env, blacs_env)
CALL timeset(routineN, handle)
CALL get_qs_env(qs_env, dbcsr_dist=dbcsr_dist, particle_set=particle_set, natom=natom, &
para_env=para_env, blacs_env=blacs_env, cell=cell)
nnb = SIZE(neighbors)
ALLOCATE (nsgf(nnb), is_neighbor(natom), radius(nnb))
is_neighbor = .FALSE.
DO inb = 1, nnb
iat = neighbors(inb)
ikind = particle_set(iat)%atomic_kind%kind_number
CALL get_gto_basis_set(basis_set_ri(ikind)%gto_basis_set, nsgf=nsgf(inb), kind_radius=radius(inb))
is_neighbor(iat) = .TRUE.
END DO
!Create the ri_sinv matrix based on some arbitrary dbcsr_dist
CALL dbcsr_distribution_get(dbcsr_dist, group=group, pgrid=pgrid, nprows=nprows, npcols=npcols)
ALLOCATE (row_dist(nnb), col_dist(nnb))
DO inb = 1, nnb
row_dist(inb) = MODULO(nprows - inb, nprows)
col_dist(inb) = MODULO(npcols - inb, npcols)
END DO
CALL dbcsr_distribution_new(sinv_dist, group=group, pgrid=pgrid, row_dist=row_dist, &
col_dist=col_dist)
CALL dbcsr_create(matrix=ri_sinv, name="RI_SINV", matrix_type="S", dist=sinv_dist, &
row_blk_size=nsgf, col_blk_size=nsgf)
!reserving the blocks in the correct pattern
DO inb = 1, nnb
ri = pbc(particle_set(neighbors(inb))%r, cell)
DO jnb = inb, nnb
!do the atom overlap ?
rj = pbc(particle_set(neighbors(jnb))%r, cell)
rij = pbc(ri, rj, cell)
IF (SUM(rij**2) > (radius(inb) + radius(jnb))**2) CYCLE
CALL dbcsr_get_stored_coordinates(ri_sinv, inb, jnb, blk)
IF (para_env%mepos == blk) THEN
ALLOCATE (block_risinv(nsgf(inb), nsgf(jnb)))
block_risinv = 0.0_dp
CALL dbcsr_put_block(ri_sinv, inb, jnb, block_risinv)
DEALLOCATE (block_risinv)
END IF
END DO
END DO
CALL dbcsr_finalize(ri_sinv)
CALL dbcsr_distribution_release(sinv_dist)
DEALLOCATE (row_dist, col_dist)
!prepare the send and recv buffers we will need for change of dist between the two matrices
!worst case scenario: all neighbors are on same procs => need to send nnb**2 messages
ALLOCATE (send_buff(nnb**2), recv_buff(nnb**2))
ALLOCATE (send_req(nnb**2), recv_req(nnb**2))
is = 0; ir = 0
!Loop over the whole RI overlap matrix and pick the blocks we need
CALL dbcsr_iterator_start(iter, whole_s)
DO WHILE (dbcsr_iterator_blocks_left(iter))
CALL dbcsr_iterator_next_block(iter, row=iat, column=jat, blk=blk)
CALL dbcsr_get_block_p(whole_s, iat, jat, block_whole, found_whole)
!only interested in neighbors
IF (.NOT. found_whole) CYCLE
IF (.NOT. is_neighbor(iat)) CYCLE
IF (.NOT. is_neighbor(jat)) CYCLE
inb = idx_to_nb(iat)
jnb = idx_to_nb(jat)
!If blocks are on the same proc for both matrices, simply copy
CALL dbcsr_get_block_p(ri_sinv, inb, jnb, block_risinv, found_risinv)
IF (found_risinv) THEN
CALL dcopy(nsgf(inb)*nsgf(jnb), block_whole, 1, block_risinv, 1)
ELSE
!send the block with unique tag to the proc where inb,jnb is in ri_sinv
CALL dbcsr_get_stored_coordinates(ri_sinv, inb, jnb, dest)
is = is + 1
send_buff(is)%array => block_whole
tag = natom*iat + jat
CALL mp_isend(msgin=send_buff(is)%array, dest=dest, comm=group, request=send_req(is), tag=tag)
END IF
END DO !dbcsr iter
CALL dbcsr_iterator_stop(iter)
!Loop over ri_sinv and receive all those blocks
CALL dbcsr_iterator_start(iter, ri_sinv)
DO WHILE (dbcsr_iterator_blocks_left(iter))
CALL dbcsr_iterator_next_block(iter, row=inb, column=jnb, blk=blk)
CALL dbcsr_get_block_p(ri_sinv, inb, jnb, block_risinv, found_risinv)
IF (.NOT. found_risinv) CYCLE
iat = neighbors(inb)
jat = neighbors(jnb)
!If blocks are on the same proc on both matrices do nothing
CALL dbcsr_get_stored_coordinates(whole_s, iat, jat, source)
IF (para_env%mepos == source) CYCLE
tag = natom*iat + jat
ir = ir + 1
recv_buff(ir)%array => block_risinv
CALL mp_irecv(msgout=recv_buff(ir)%array, source=source, request=recv_req(ir), &
tag=tag, comm=group)
END DO
CALL dbcsr_iterator_stop(iter)
!make sure that all comm is over before proceeding
CALL mp_waitall(send_req(1:is))
CALL mp_waitall(recv_req(1:ir))
!Invert
CALL cp_dbcsr_cholesky_decompose(ri_sinv, para_env=para_env, blacs_env=blacs_env)
CALL cp_dbcsr_cholesky_invert(ri_sinv, para_env=para_env, blacs_env=blacs_env, upper_to_full=.TRUE.)
CALL dbcsr_filter(ri_sinv, 1.E-10_dp) !make sure ri_sinv is sparse coming out of fm routines
CALL dbcsr_replicate_all(ri_sinv)
!clean-up
DEALLOCATE (nsgf)
CALL timestop(handle)
END SUBROUTINE get_exat_ri_sinv
! **************************************************************************************************
!> \brief Compute the coefficients to project the density on a partial RI_XAS basis
!> \param xas_atom_env ...
!> \param qs_env ...
!> \note The density is n = sum_ab P_ab*phi_a*phi_b, the RI basis covers the products of orbital sgfs
!> => n = sum_ab sum_cd P_ab (phi_a phi_b xi_c) S_cd^-1 xi_d
!> = sum_d coeff_d xi_d , where xi are the RI basis func.
!> In this case, with the partial RI projection, the RI basis is restricted to an excited atom
!> and its neighbors at a time. Leads to smaller overlap matrix to invert and less 3-center
!> overlap to compute. The procedure is repeated for each excited atom
!> This is a test implementation for tensors
! **************************************************************************************************
SUBROUTINE calculate_density_coeffs(xas_atom_env, qs_env)
TYPE(xas_atom_env_type), POINTER :: xas_atom_env
TYPE(qs_environment_type), POINTER :: qs_env
CHARACTER(len=*), PARAMETER :: routineN = 'calculate_density_coeffs', &
routineP = moduleN//':'//routineN
INTEGER :: blk, exat, handle, i, iat, iatom, iex, &
inb, ind(3), ispin, jatom, jnb, katom, &
natom, nex, nkind, nnb, nspins, &
output_unit, ri_at
INTEGER, ALLOCATABLE, DIMENSION(:) :: blk_size_orb, blk_size_ri, idx_to_nb, &
neighbors
INTEGER, DIMENSION(:), POINTER :: all_ri_atoms
LOGICAL :: pmat_found, pmat_foundt, sinv_found, &
sinv_foundt, tensor_found, unique
REAL(dp) :: factor, prefac
REAL(dp), ALLOCATABLE, DIMENSION(:) :: work2
REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: work1
REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: t_block
REAL(dp), DIMENSION(:, :), POINTER :: pmat_block, pmat_blockt, sinv_block, &
sinv_blockt
TYPE(cp_blacs_env_type), POINTER :: blacs_env
TYPE(cp_para_env_type), POINTER :: para_env
TYPE(dbcsr_distribution_type) :: opt_dbcsr_dist
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: overlap, rho_ao, rho_redist
TYPE(dbcsr_t_iterator_type) :: iter
TYPE(dbcsr_t_type) :: pqX
TYPE(dbcsr_type) :: ri_sinv
TYPE(dbcsr_type), POINTER :: ri_mats
TYPE(distribution_2d_type), POINTER :: opt_dist2d
TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_orb, basis_set_ri
TYPE(libint_potential_type) :: pot
TYPE(neighbor_list_set_p_type), DIMENSION(:), &
POINTER :: ab_list, ac_list, sab_ri
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(qs_rho_type), POINTER :: rho
NULLIFY (qs_kind_set, basis_set_ri, basis_set_orb, ac_list, rho, rho_ao, sab_ri, ri_mats)
NULLIFY (particle_set, para_env, all_ri_atoms, overlap, pmat_blockt)
NULLIFY (blacs_env, pmat_block, ab_list, opt_dist2d, rho_redist, sinv_block, sinv_blockt)
!Idea: We don't do a full RI here as it would be too expensive in many ways (inversion of a
! large matrix, many 3-center overlaps, density coefficients for all atoms, etc...)
! Instead, we go excited atom by excited atom and only do a RI expansion on its basis
! and that of its closest neighbors (defined through RI_RADIUS), such that we only have
! very small matrices to invert and only a few (abc) overlp integrals with c on the
! excited atom its neighbors. This is physically sound since we only need the density
! well defined on the excited atom grid as we do (aI|P)*(P|Q)^-1*(Q|fxc|R)*(R|S)^-1*(S|Jb)
CALL timeset(routineN, handle)
CALL get_qs_env(qs_env, nkind=nkind, qs_kind_set=qs_kind_set, rho=rho, &
natom=natom, particle_set=particle_set, &
para_env=para_env, blacs_env=blacs_env)
nspins = xas_atom_env%nspins
nex = SIZE(xas_atom_env%excited_atoms)
CALL qs_rho_get(rho, rho_ao=rho_ao)
! Create the needed neighbor list and basis set lists.
ALLOCATE (basis_set_ri(nkind))
ALLOCATE (basis_set_orb(nkind))
CALL basis_set_list_setup(basis_set_ri, "RI_XAS", qs_kind_set)
CALL basis_set_list_setup(basis_set_orb, "ORB", qs_kind_set)
! Compute the RI overlap matrix on the whole system
CALL build_xas_tdp_ovlp_nl(sab_ri, basis_set_ri, basis_set_ri, qs_env)
CALL build_overlap_matrix_simple(qs_env%ks_env, overlap, basis_set_ri, basis_set_ri, sab_ri)
ri_mats => overlap(1)%matrix
CALL release_neighbor_list_sets(sab_ri)
! Get the neighbors of the excited atoms (= all the atoms where density coeffs are needed)
CALL find_neighbors(xas_atom_env%excited_atoms, ri_mats, xas_atom_env%ri_radius, &
qs_env, all_neighbors=all_ri_atoms, neighbor_set=xas_atom_env%exat_neighbors)
!keep in mind that double occupation is included in rho_ao in case of closed-shell
factor = 0.5_dp; IF (nspins == 2) factor = 1.0_dp
! Allocate space for the projected density coefficients. On all ri_atoms.
! Note: the sub-region where we project the density changes from excited atom to excited atom
! => need different sets of RI coeffs
ALLOCATE (blk_size_ri(natom))
CALL get_particle_set(particle_set, qs_kind_set, nsgf=blk_size_ri, basis=basis_set_ri)
ALLOCATE (xas_atom_env%ri_dcoeff(natom, nspins, nex))
DO iex = 1, nex
DO ispin = 1, nspins
DO iat = 1, natom
NULLIFY (xas_atom_env%ri_dcoeff(iat, ispin, iex)%array)
IF ((.NOT. ANY(xas_atom_env%exat_neighbors(iex)%array == iat)) &
.AND. (.NOT. xas_atom_env%excited_atoms(iex) == iat)) CYCLE
ALLOCATE (xas_atom_env%ri_dcoeff(iat, ispin, iex)%array(blk_size_ri(iat)))
xas_atom_env%ri_dcoeff(iat, ispin, iex)%array = 0.0_dp
END DO
END DO
END DO
output_unit = cp_logger_get_default_io_unit()
IF (output_unit > 0) THEN
WRITE (output_unit, FMT="(/,T7,A,/,T7,A)") &
"Excited atom, natoms in RI_REGION:", &
"---------------------------------"
END IF
!We go atom by atom, first compute the optimal distribution for the integrals, then
!the integrals themselves that we put into a tensor, then the contraction with the
!density
ALLOCATE (blk_size_orb(natom))
CALL get_particle_set(particle_set, qs_kind_set, nsgf=blk_size_orb, basis=basis_set_orb)
DO iex = 1, nex
!get neighbors of current atom
exat = xas_atom_env%excited_atoms(iex)
nnb = 1 + SIZE(xas_atom_env%exat_neighbors(iex)%array)
ALLOCATE (neighbors(nnb))
neighbors(1) = exat
neighbors(2:nnb) = xas_atom_env%exat_neighbors(iex)%array(:)
CALL sort_unique(neighbors, unique)
!link the atoms to their position in neighbors
ALLOCATE (idx_to_nb(natom))
idx_to_nb = 0
DO inb = 1, nnb
idx_to_nb(neighbors(inb)) = inb
END DO
IF (output_unit > 0) THEN
WRITE (output_unit, FMT="(T7,I12,I21)") &
exat, nnb
END IF
!Get the optimal distribution_2d for the overlap integrals (abc), centers c on the current
!excited atom and its neighbors defined by RI_RADIUS
CALL build_xas_tdp_ovlp_nl(ab_list, basis_set_orb, basis_set_orb, qs_env)
CALL build_xas_tdp_3c_nl(ac_list, basis_set_orb, basis_set_ri, do_potential_id, &
qs_env, excited_atoms=neighbors)
CALL get_opt_3c_dist2d(opt_dist2d, ab_list, ac_list, basis_set_orb, basis_set_orb, &
basis_set_ri, qs_env)
CALL release_neighbor_list_sets(ab_list)
CALL release_neighbor_list_sets(ac_list)
!get the ab and ac neighbor lists based on the optimized distribution
CALL build_xas_tdp_ovlp_nl(ab_list, basis_set_orb, basis_set_orb, qs_env, ext_dist2d=opt_dist2d)
CALL build_xas_tdp_3c_nl(ac_list, basis_set_orb, basis_set_ri, do_potential_id, &
qs_env, excited_atoms=neighbors, ext_dist2d=opt_dist2d)
!get the AO density matrix in the optimized distribution
CALL cp_dbcsr_dist2d_to_dist(opt_dist2d, opt_dbcsr_dist)
ALLOCATE (rho_redist(nspins))
DO ispin = 1, nspins
ALLOCATE (rho_redist(ispin)%matrix)
CALL dbcsr_create(rho_redist(ispin)%matrix, template=rho_ao(ispin)%matrix, dist=opt_dbcsr_dist)
CALL dbcsr_complete_redistribute(rho_ao(ispin)%matrix, rho_redist(ispin)%matrix)
END DO
!Compute the 3-center overlap integrals
pot%potential_type = do_potential_id
CALL create_pqX_tensor(pqX, ab_list, ac_list, opt_dbcsr_dist, blk_size_orb, blk_size_orb, &
blk_size_ri)
CALL fill_pqX_tensor(pqX, ab_list, ac_list, basis_set_orb, basis_set_orb, basis_set_ri, &
pot, qs_env)
!Compute the RI inverse overlap matrix on the reduced RI basis that spans the excited
!atom and its neighbors, ri_sinv is replicated over all procs
CALL get_exat_ri_sinv(ri_sinv, ri_mats, neighbors, idx_to_nb, basis_set_ri, qs_env)
!Do the actual contraction: coeff_y = sum_pq sum_x P_pq (phi_p phi_q xi_x) S_xy^-1
!TODO: re-enable OMP parallelization once/if the tensor iterator becomes thread safe
!!$OMP PARALLEL DEFAULT(NONE) &
!!$OMP SHARED (pqX,nspins,blk_size_ri,nnb,neighbors,ri_sinv,factor,xas_atom_env,iex,rho_redist,idx_to_nb)&
!!$OMP PRIVATE(iter,ind,blk,t_block,tensor_found,iatom,jatom,inb,prefac,ispin,work1,pmat_block, &
!!$OMP pmat_blockt,pmat_found,pmat_foundt,i,jnb,ri_at,sinv_block,sinv_blockt,sinv_found,&
!!$OMP sinv_foundt,work2,katom)
CALL dbcsr_t_iterator_start(iter, pqX)
DO WHILE (dbcsr_t_iterator_blocks_left(iter))
CALL dbcsr_t_iterator_next_block(iter, ind, blk)
!!$OMP CRITICAL (get_t_block)
CALL dbcsr_t_get_block(pqX, ind, t_block, tensor_found)
!!$OMP END CRITICAL (get_t_block)
iatom = ind(1)
jatom = ind(2)
katom = ind(3)
inb = idx_to_nb(katom)
!non-diagonal elements need to be counted twice
prefac = 2.0_dp
IF (iatom == jatom) prefac = 1.0_dp
DO ispin = 1, nspins
!rho_redist is symmetric, block can be in either location
CALL dbcsr_get_block_p(rho_redist(ispin)%matrix, iatom, jatom, pmat_block, pmat_found)
CALL dbcsr_get_block_p(rho_redist(ispin)%matrix, jatom, iatom, pmat_blockt, pmat_foundt)
IF ((.NOT. pmat_found) .AND. (.NOT. pmat_foundt)) CYCLE
ALLOCATE (work1(blk_size_ri(katom), 1))
work1 = 0.0_dp
!first contraction with the density matrix
IF (pmat_found) THEN
DO i = 1, blk_size_ri(katom)
work1(i, 1) = prefac*SUM(pmat_block(:, :)*t_block(:, :, i))
END DO
ELSE
DO i = 1, blk_size_ri(katom)
work1(i, 1) = prefac*SUM(TRANSPOSE(pmat_blockt(:, :))*t_block(:, :, i))
END DO
END IF
!loop over neighbors
DO jnb = 1, nnb
ri_at = neighbors(jnb)
!ri_sinv is a symmetric matrix => actual block is one of the two
CALL dbcsr_get_block_p(ri_sinv, inb, jnb, sinv_block, sinv_found)
CALL dbcsr_get_block_p(ri_sinv, jnb, inb, sinv_blockt, sinv_foundt)
IF ((.NOT. sinv_found) .AND. (.NOT. sinv_foundt)) CYCLE
!second contraction with the inverse RI overlap
ALLOCATE (work2(SIZE(xas_atom_env%ri_dcoeff(ri_at, ispin, iex)%array)))
work2 = 0.0_dp
IF (sinv_found) THEN
DO i = 1, blk_size_ri(katom)
CALL daxpy(SIZE(work2), factor*work1(i, 1), sinv_block(i, :), 1, work2(:), 1)
END DO
ELSE
DO i = 1, blk_size_ri(katom)
CALL daxpy(SIZE(work2), factor*work1(i, 1), sinv_blockt(:, i), 1, work2(:), 1)
END DO
END IF
!!$OMP CRITICAL (ri_dcoeff_update)
xas_atom_env%ri_dcoeff(ri_at, ispin, iex)%array(:) = &
xas_atom_env%ri_dcoeff(ri_at, ispin, iex)%array(:) + work2(:)
!!$OMP END CRITICAL (ri_dcoeff_update)
DEALLOCATE (work2)
END DO !jnb
DEALLOCATE (work1)
END DO
DEALLOCATE (t_block)
END DO !iter
CALL dbcsr_t_iterator_stop(iter)
!!$OMP END PARALLEL
!clean-up
CALL dbcsr_release(ri_sinv)
CALL dbcsr_t_destroy(pqX)
CALL distribution_2d_release(opt_dist2d)
CALL dbcsr_deallocate_matrix_set(rho_redist)
CALL dbcsr_distribution_release(opt_dbcsr_dist)
CALL release_neighbor_list_sets(ab_list)
CALL release_neighbor_list_sets(ac_list)
DEALLOCATE (neighbors, idx_to_nb)
END DO !iex
!making sure all procs have the same info
DO iex = 1, nex
DO ispin = 1, nspins
DO iat = 1, natom
IF ((.NOT. ANY(xas_atom_env%exat_neighbors(iex)%array == iat)) &
.AND. (.NOT. xas_atom_env%excited_atoms(iex) == iat)) CYCLE
CALL mp_sum(xas_atom_env%ri_dcoeff(iat, ispin, iex)%array, para_env%group)
END DO !iat
END DO !ispin
END DO !iex
! clean-up
CALL dbcsr_deallocate_matrix_set(overlap)
DEALLOCATE (basis_set_ri, basis_set_orb, all_ri_atoms)
CALL timestop(handle)
END SUBROUTINE calculate_density_coeffs
! **************************************************************************************************
!> \brief Evaluates the density on a given atomic grid
!> \param rho_set where the densities are stored
!> \param ri_dcoeff the arrays containing the RI density coefficients of this atom, for each spin
!> \param atom_kind the kind of the atom in question
!> \param do_gga whether the gradient of the density should also be put on the grid
!> \param batch_info how the so are distributed
!> \param xas_atom_env ...
!> \param qs_env ...
!> \note The density is expressed as n = sum_d coeff_d*xi_d. Knowing the coordinate of each grid
!> grid point, one can simply evaluate xi_d(r)
! **************************************************************************************************
SUBROUTINE put_density_on_atomic_grid(rho_set, ri_dcoeff, atom_kind, do_gga, batch_info, &
xas_atom_env, qs_env)
TYPE(xc_rho_set_type), POINTER :: rho_set
TYPE(cp_1d_r_p_type), DIMENSION(:) :: ri_dcoeff
INTEGER, INTENT(IN) :: atom_kind
LOGICAL, INTENT(IN) :: do_gga
TYPE(batch_info_type) :: batch_info
TYPE(xas_atom_env_type), POINTER :: xas_atom_env
TYPE(qs_environment_type), POINTER :: qs_env
CHARACTER(len=*), PARAMETER :: routineN = 'put_density_on_atomic_grid', &
routineP = moduleN//':'//routineN
INTEGER :: dir, handle, ipgf, iset, iso, iso_proc, &
ispin, maxso, n, na, nr, nset, nsgfi, &
nsoi, nspins, sgfi, starti
INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf, nsgf_set
INTEGER, DIMENSION(:, :), POINTER :: first_sgf
REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: so, work1, work2
REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: dso
REAL(dp), DIMENSION(:, :), POINTER :: dgr1, dgr2, ga, gr, ri_sphi_so, zet
REAL(dp), DIMENSION(:, :, :), POINTER :: dga1, dga2, rhoa, rhob
TYPE(cp_1d_r_p_type), DIMENSION(:), POINTER :: ri_dcoeff_so
TYPE(grid_atom_type), POINTER :: grid_atom
TYPE(gto_basis_set_type), POINTER :: ri_basis
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
NULLIFY (grid_atom, ri_basis, qs_kind_set, lmax, npgf, zet, nsgf_set, ri_sphi_so)
NULLIFY (lmin, first_sgf, rhoa, rhob, ga, gr, dgr1, dgr2, dga1, dga2, ri_dcoeff_so)
CALL timeset(routineN, handle)
! Strategy: it makes sense to evaluate the spherical orbital on the grid (because of symmetry)
! From there, one can directly contract into sgf using ri_sphi_so and then take the weight
! The spherical orbital were precomputed and split in a purely radial and a purely
! angular part. The full values on each grid point are obtain through gemm
! Generalities
CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set)
CALL get_qs_kind(qs_kind_set(atom_kind), basis_set=ri_basis, basis_type="RI_XAS")
CALL get_gto_basis_set(ri_basis, lmax=lmax, npgf=npgf, zet=zet, nset=nset, nsgf_set=nsgf_set, &
first_sgf=first_sgf, lmin=lmin, maxso=maxso)
! Get the grid and the info we need from it
grid_atom => xas_atom_env%grid_atom_set(atom_kind)%grid_atom
na = grid_atom%ng_sphere
nr = grid_atom%nr
n = na*nr
nspins = xas_atom_env%nspins
ri_sphi_so => xas_atom_env%ri_sphi_so(atom_kind)%array
! Point to the rho_set densities
rhoa => rho_set%rhoa
rhob => rho_set%rhob
rhoa = 0.0_dp; rhob = 0.0_dp;
IF (do_gga) THEN
DO dir = 1, 3
rho_set%drhoa(dir)%array = 0.0_dp
rho_set%drhob(dir)%array = 0.0_dp
END DO
END IF
! Point to the precomputed SO
ga => xas_atom_env%ga(atom_kind)%array
gr => xas_atom_env%gr(atom_kind)%array
IF (do_gga) THEN
dga1 => xas_atom_env%dga1(atom_kind)%array
dga2 => xas_atom_env%dga2(atom_kind)%array
dgr1 => xas_atom_env%dgr1(atom_kind)%array
dgr2 => xas_atom_env%dgr2(atom_kind)%array
END IF
! Need to express the ri_dcoeffs in terms of so (and not sgf)
ALLOCATE (ri_dcoeff_so(nspins))
DO ispin = 1, nspins
ALLOCATE (ri_dcoeff_so(ispin)%array(nset*maxso))
ri_dcoeff_so(ispin)%array = 0.0_dp
!for a given so, loop over sgf and sum
DO iset = 1, nset
sgfi = first_sgf(1, iset)
nsoi = npgf(iset)*nsoset(lmax(iset))
nsgfi = nsgf_set(iset)
ALLOCATE (work1(nsoi, 1), work2(nsgfi, 1))
work2(:, 1) = ri_dcoeff(ispin)%array(sgfi:sgfi + nsgfi - 1)
CALL dgemm('N', 'N', nsoi, 1, nsgfi, 1.0_dp, ri_sphi_so(1:nsoi, sgfi:sgfi + nsgfi - 1), &
nsoi, work2, nsgfi, 0.0_dp, work1, nsoi)
ri_dcoeff_so(ispin)%array((iset - 1)*maxso + 1:(iset - 1)*maxso + nsoi) = work1(:, 1)
DEALLOCATE (work1, work2)
END DO
END DO
!allocate space to store the spherical orbitals on the grid
ALLOCATE (so(na, nr))
IF (do_gga) ALLOCATE (dso(na, nr, 3))
! Loop over the spherical orbitals on this proc
DO iso_proc = 1, batch_info%nso_proc(atom_kind)
iset = batch_info%so_proc_info(atom_kind)%array(1, iso_proc)
ipgf = batch_info%so_proc_info(atom_kind)%array(2, iso_proc)
iso = batch_info%so_proc_info(atom_kind)%array(3, iso_proc)
IF (iso < 0) CYCLE
starti = (iset - 1)*maxso + (ipgf - 1)*nsoset(lmax(iset))
!the spherical orbital on the grid
CALL dgemm('N', 'T', na, nr, 1, 1.0_dp, ga(:, iso_proc:iso_proc), na, &
gr(:, iso_proc:iso_proc), nr, 0.0_dp, so(:, :), na)
!the gradient on the grid
IF (do_gga) THEN
DO dir = 1, 3
CALL dgemm('N', 'T', na, nr, 1, 1.0_dp, dga1(:, iso_proc:iso_proc, dir), na, &
dgr1(:, iso_proc:iso_proc), nr, 0.0_dp, dso(:, :, dir), na)
CALL dgemm('N', 'T', na, nr, 1, 1.0_dp, dga2(:, iso_proc:iso_proc, dir), na, &
dgr2(:, iso_proc:iso_proc), nr, 1.0_dp, dso(:, :, dir), na)
END DO
END IF
!put the so on the grid with the approriate coefficients and sum
CALL daxpy(n, ri_dcoeff_so(1)%array(starti + iso), so(:, :), 1, rhoa(:, :, 1), 1)
IF (nspins == 2) THEN
CALL daxpy(n, ri_dcoeff_so(2)%array(starti + iso), so(:, :), 1, rhob(:, :, 1), 1)
END IF
IF (do_gga) THEN
!put the gradient of the so on the grid with correspond RI coeff
DO dir = 1, 3
CALL daxpy(n, ri_dcoeff_so(1)%array(starti + iso), dso(:, :, dir), &
1, rho_set%drhoa(dir)%array(:, :, 1), 1)
IF (nspins == 2) THEN
CALL daxpy(n, ri_dcoeff_so(2)%array(starti + iso), dso(:, :, dir), &
1, rho_set%drhob(dir)%array(:, :, 1), 1)
END IF
END DO !dir
END IF !do_gga
END DO
! Treat spin restricted case (=> copy alpha into beta)
IF (nspins == 1) THEN
CALL dcopy(n, rhoa(:, :, 1), 1, rhob(:, :, 1), 1)
IF (do_gga) THEN
DO dir = 1, 3
CALL dcopy(n, rho_set%drhoa(dir)%array(:, :, 1), 1, rho_set%drhob(dir)%array(:, :, 1), 1)
END DO
END IF
END IF
! Note: sum over procs is done outside
! clean-up
DO ispin = 1, nspins
DEALLOCATE (ri_dcoeff_so(ispin)%array)
END DO
DEALLOCATE (ri_dcoeff_so)
CALL timestop(handle)
END SUBROUTINE put_density_on_atomic_grid
! **************************************************************************************************
!> \brief Adds the density of a given source atom with source kind (with ri_dcoeff) on the atomic
!> grid belonging to another target atom of target kind. The evaluations of the basis
!> function first requires the evaluation of the x,y,z coordinates on each grid point of
!> target atom wrt to the position of source atom
!> \param rho_set where the densities are stored
!> \param ri_dcoeff the arrays containing the RI density coefficient of source_iat, for each spin
!> \param source_iat the index of the source atom
!> \param source_ikind the kind of the source atom
!> \param target_iat the index of the target atom
!> \param target_ikind the kind of the target atom
!> \param sr starting r index for the local grid
!> \param er ending r index for the local grid
!> \param do_gga whether the gradient of the density is needed
!> \param xas_atom_env ...
!> \param qs_env ...
! **************************************************************************************************
SUBROUTINE put_density_on_other_grid(rho_set, ri_dcoeff, source_iat, source_ikind, target_iat, &
target_ikind, sr, er, do_gga, xas_atom_env, qs_env)
TYPE(xc_rho_set_type), POINTER :: rho_set
TYPE(cp_1d_r_p_type), DIMENSION(:) :: ri_dcoeff
INTEGER, INTENT(IN) :: source_iat, source_ikind, target_iat, &
target_ikind, sr, er
LOGICAL, INTENT(IN) :: do_gga
TYPE(xas_atom_env_type), POINTER :: xas_atom_env
TYPE(qs_environment_type), POINTER :: qs_env
CHARACTER(len=*), PARAMETER :: routineN = 'put_density_on_other_grid', &
routineP = moduleN//':'//routineN
INTEGER :: dir, handle, ia, ico, ipgf, ir, iset, &
isgf, lx, ly, lz, n, na, nr, nset, &
nspins, sgfi, start
INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf, nsgf_set
INTEGER, DIMENSION(:, :), POINTER :: first_sgf
REAL(dp) :: rmom
REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: sgf
REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: co, dsgf, pos
REAL(dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: dco
REAL(dp), DIMENSION(3) :: rs, rst, rt
REAL(dp), DIMENSION(:, :), POINTER :: ri_sphi, zet
REAL(dp), DIMENSION(:, :, :), POINTER :: rhoa, rhob
TYPE(cell_type), POINTER :: cell
TYPE(grid_atom_type), POINTER :: grid_atom
TYPE(gto_basis_set_type), POINTER :: ri_basis
TYPE(harmonics_atom_type), POINTER :: harmonics
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
NULLIFY (qs_kind_set, ri_basis, lmax, npgf, nsgf_set, lmin, zet, first_sgf, grid_atom)
NULLIFY (harmonics, rhoa, rhob, particle_set, cell, ri_sphi)
!Same logic as the put_density_on_own_grid routine. Loop over orbitals, put them on the grid,
!contract into sgf and daxpy with coeff. Notable difference: use cartesian orbitals instead of
!spherical, since the center of the grid is not the origin and thus, spherical symmetry can't
!be exploited so well
CALL timeset(routineN, handle)
!Generalities
CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, particle_set=particle_set, cell=cell)
!want basis of the source atom
CALL get_qs_kind(qs_kind_set(source_ikind), basis_set=ri_basis, basis_type="RI_XAS")
CALL get_gto_basis_set(ri_basis, lmax=lmax, npgf=npgf, zet=zet, nset=nset, nsgf_set=nsgf_set, &
first_sgf=first_sgf, lmin=lmin, sphi=ri_sphi)
! Want the grid and harmonics of the target atom
grid_atom => xas_atom_env%grid_atom_set(target_ikind)%grid_atom
harmonics => xas_atom_env%harmonics_atom_set(target_ikind)%harmonics_atom
na = grid_atom%ng_sphere
nr = er - sr + 1
n = na*nr
nspins = xas_atom_env%nspins
! Point to the rho_set densities
rhoa => rho_set%rhoa
rhob => rho_set%rhob
! Need the source-target position vector
rs = pbc(particle_set(source_iat)%r, cell)
rt = pbc(particle_set(target_iat)%r, cell)
rst = pbc(rs, rt, cell)
! Precompute the positions on the target grid
ALLOCATE (pos(na, sr:er, 4))
!$OMP PARALLEL DO COLLAPSE(2) SCHEDULE(STATIC) DEFAULT(NONE), &
!$OMP SHARED(na,sr,er,pos,harmonics,grid_atom,rst), &
!$OMP PRIVATE(ia,ir)
DO ir = sr, er
DO ia = 1, na
pos(ia, ir, 1:3) = harmonics%a(:, ia)*grid_atom%rad(ir) + rst
pos(ia, ir, 4) = pos(ia, ir, 1)**2 + pos(ia, ir, 2)**2 + pos(ia, ir, 3)**2
END DO
END DO
!$OMP END PARALLEL DO
! Loop over the cartesian gaussian functions and evaluate them
DO iset = 1, nset
!allocate space to store the cartesian orbtial on the grid
ALLOCATE (co(na, sr:er, npgf(iset)*ncoset(lmax(iset))))
IF (do_gga) ALLOCATE (dco(na, sr:er, 3, npgf(iset)*ncoset(lmax(iset))))
!$OMP PARALLEL DEFAULT(NONE), &
!$OMP SHARED(co,npgf,ncoset,lmax,lmin,indco,pos,zet,iset,na,sr,er,do_gga,dco), &
!$OMP PRIVATE(ipgf,start,ico,lx,ly,lz,ia,ir,rmom)
!$OMP DO COLLAPSE(2) SCHEDULE(STATIC)
DO ir = sr, er
DO ia = 1, na
co(ia, ir, :) = 0.0_dp
IF (do_gga) THEN
dco(ia, ir, :, :) = 0.0_dp
END IF
END DO
END DO
!$OMP END DO NOWAIT
DO ipgf = 1, npgf(iset)
start = (ipgf - 1)*ncoset(lmax(iset))
!loop over the cartesian orbitals
DO ico = ncoset(lmin(iset) - 1) + 1, ncoset(lmax(iset))
lx = indco(1, ico)
ly = indco(2, ico)
lz = indco(3, ico)
! compute g = x**lx * y**ly * z**lz * exp(-zet * r**2)
!$OMP DO COLLAPSE(2) SCHEDULE(STATIC)
DO ir = sr, er
DO ia = 1, na
rmom = EXP(-zet(ipgf, iset)*pos(ia, ir, 4))
IF (lx /= 0) rmom = rmom*pos(ia, ir, 1)**lx
IF (ly /= 0) rmom = rmom*pos(ia, ir, 2)**ly
IF (lz /= 0) rmom = rmom*pos(ia, ir, 3)**lz
co(ia, ir, start + ico) = rmom
!co(ia, ir, start + ico) = pos(ia, ir, 1)**lx*pos(ia, ir, 2)**ly*pos(ia, ir, 3)**lz &
! *EXP(-zet(ipgf, iset)*pos(ia, ir, 4))
END DO
END DO
!$OMP END DO NOWAIT
IF (do_gga) THEN
!the gradient: dg_x = lx*x**(lx-1) * y**ly * z**lz * exp(-zet * r**2)
! -2*zet* x**(lx+1) * y**ly * z**lz * exp(-zet * r**2)
! = (lx*x**(lx-1) - 2*zet*x**(lx+1)) * y**ly * z**lz * exp(-zet * r**2)
!x direction, special case if lx == 0
IF (lx == 0) THEN
!$OMP DO COLLAPSE(2) SCHEDULE(STATIC)
DO ir = sr, er
DO ia = 1, na
rmom = -2.0_dp*pos(ia, ir, 1)*zet(ipgf, iset)*EXP(-zet(ipgf, iset)*pos(ia, ir, 4))
IF (ly /= 0) rmom = rmom*pos(ia, ir, 2)**ly
IF (lz /= 0) rmom = rmom*pos(ia, ir, 3)**lz
dco(ia, ir, 1, start + ico) = rmom
! dco(ia, ir, 1, start + ico) = -2.0_dp*pos(ia, ir, 1)*zet(ipgf, iset) &
! *pos(ia, ir, 2)**ly*pos(ia, ir, 3)**lz &
! *EXP(-zet(ipgf, iset)*pos(ia, ir, 4))
END DO
END DO
!$OMP END DO NOWAIT
ELSE
!$OMP DO COLLAPSE(2) SCHEDULE(STATIC)
DO ir = sr, er
DO ia = 1, na
IF (lx /= 1) THEN
rmom = (lx*pos(ia, ir, 1)**(lx - 1) - 2.0_dp*pos(ia, ir, 1)**(lx + 1)* &
zet(ipgf, iset))*EXP(-zet(ipgf, iset)*pos(ia, ir, 4))
ELSE
rmom = (1.0_dp - 2.0_dp*pos(ia, ir, 1)**2*zet(ipgf, iset))* &
EXP(-zet(ipgf, iset)*pos(ia, ir, 4))
END IF
IF (ly /= 0) rmom = rmom*pos(ia, ir, 2)**ly
IF (lz /= 0) rmom = rmom*pos(ia, ir, 3)**lz
dco(ia, ir, 1, start + ico) = rmom
! dco(ia, ir, 1, start + ico) = (lx*pos(ia, ir, 1)**(lx - 1) &
! - 2.0_dp*pos(ia, ir, 1)**(lx + 1)*zet(ipgf, iset)) &
! *pos(ia, ir, 2)**ly*pos(ia, ir, 3)**lz &
! *EXP(-zet(ipgf, iset)*pos(ia, ir, 4))
END DO
END DO
!$OMP END DO NOWAIT
END IF !lx == 0
!y direction, special case if ly == 0
IF (ly == 0) THEN
!$OMP DO COLLAPSE(2) SCHEDULE(STATIC)
DO ir = sr, er
DO ia = 1, na
rmom = -2.0_dp*pos(ia, ir, 2)*zet(ipgf, iset)*EXP(-zet(ipgf, iset)*pos(ia, ir, 4))
IF (lx /= 0) rmom = rmom*pos(ia, ir, 1)**lx
IF (lz /= 0) rmom = rmom*pos(ia, ir, 3)**lz
dco(ia, ir, 2, start + ico) = rmom
! dco(ia, ir, 2, start + ico) = -2.0_dp*pos(ia, ir, 2)*zet(ipgf, iset) &
! *pos(ia, ir, 1)**lx*pos(ia, ir, 3)**lz &
! *EXP(-zet(ipgf, iset)*pos(ia, ir, 4))
END DO
END DO
!$OMP END DO NOWAIT
ELSE
!$OMP DO COLLAPSE(2) SCHEDULE(STATIC)
DO ir = sr, er
DO ia = 1, na
IF (ly /= 1) THEN
rmom = (ly*pos(ia, ir, 2)**(ly - 1) - 2.0_dp*pos(ia, ir, 2)**(ly + 1)*zet(ipgf, iset)) &
*EXP(-zet(ipgf, iset)*pos(ia, ir, 4))
ELSE
rmom = (1.0_dp - 2.0_dp*pos(ia, ir, 2)**2*zet(ipgf, iset)) &
*EXP(-zet(ipgf, iset)*pos(ia, ir, 4))
END IF
IF (lx /= 0) rmom = rmom*pos(ia, ir, 1)**lx
IF (lz /= 0) rmom = rmom*pos(ia, ir, 3)**lz
dco(ia, ir, 2, start + ico) = rmom
! dco(ia, ir, 2, start + ico) = (ly*pos(ia, ir, 2)**(ly - 1) &
! - 2.0_dp*pos(ia, ir, 2)**(ly + 1)*zet(ipgf, iset)) &
! *pos(ia, ir, 1)**lx*pos(ia, ir, 3)**lz &
! *EXP(-zet(ipgf, iset)*pos(ia, ir, 4))
END DO
END DO
!$OMP END DO NOWAIT
END IF !ly == 0
!z direction, special case if lz == 0
IF (lz == 0) THEN
!$OMP DO COLLAPSE(2) SCHEDULE(STATIC)
DO ir = sr, er
DO ia = 1, na
rmom = -2.0_dp*pos(ia, ir, 3)*zet(ipgf, iset)*EXP(-zet(ipgf, iset)*pos(ia, ir, 4))
IF (lx /= 0) rmom = rmom*pos(ia, ir, 1)**lx
IF (ly /= 0) rmom = rmom*pos(ia, ir, 2)**ly
dco(ia, ir, 3, start + ico) = rmom
! dco(ia, ir, 3, start + ico) = -2.0_dp*pos(ia, ir, 3)*zet(ipgf, iset) &
! *pos(ia, ir, 1)**lx*pos(ia, ir, 2)**ly &
! *EXP(-zet(ipgf, iset)*pos(ia, ir, 4))
END DO
END DO
!$OMP END DO NOWAIT
ELSE
!$OMP DO COLLAPSE(2) SCHEDULE(STATIC)
DO ir = sr, er
DO ia = 1, na
IF (lz /= 1) THEN
rmom = (lz*pos(ia, ir, 3)**(lz - 1) - 2.0_dp*pos(ia, ir, 3)**(lz + 1)* &
zet(ipgf, iset))*EXP(-zet(ipgf, iset)*pos(ia, ir, 4))
ELSE
rmom = (1.0_dp - 2.0_dp*pos(ia, ir, 3)**2*zet(ipgf, iset))* &
EXP(-zet(ipgf, iset)*pos(ia, ir, 4))
END IF
IF (lx /= 0) rmom = rmom*pos(ia, ir, 1)**lx
IF (ly /= 0) rmom = rmom*pos(ia, ir, 2)**ly
dco(ia, ir, 3, start + ico) = rmom
! dco(ia, ir, 3, start + ico) = (lz*pos(ia, ir, 3)**(lz - 1) &
! - 2.0_dp*pos(ia, ir, 3)**(lz + 1)*zet(ipgf, iset)) &
! *pos(ia, ir, 1)**lx*pos(ia, ir, 2)**ly &
! *EXP(-zet(ipgf, iset)*pos(ia, ir, 4))
END DO
END DO
!$OMP END DO NOWAIT
END IF !lz == 0
END IF !gga
END DO !ico
END DO !ipgf
!$OMP END PARALLEL
!contract the co into sgf
ALLOCATE (sgf(na, sr:er))
IF (do_gga) ALLOCATE (dsgf(na, sr:er, 3))
sgfi = first_sgf(1, iset) - 1
DO isgf = 1, nsgf_set(iset)
sgf = 0.0_dp
IF (do_gga) dsgf = 0.0_dp
DO ipgf = 1, npgf(iset)
start = (ipgf - 1)*ncoset(lmax(iset))
DO ico = ncoset(lmin(iset) - 1) + 1, ncoset(lmax(iset))
CALL daxpy(n, ri_sphi(start + ico, sgfi + isgf), co(:, sr:er, start + ico), 1, sgf(:, sr:er), 1)
END DO !ico
END DO !ipgf
!add the density to the grid
CALL daxpy(n, ri_dcoeff(1)%array(sgfi + isgf), sgf(:, sr:er), 1, rhoa(:, sr:er, 1), 1)
IF (nspins == 2) THEN
CALL daxpy(n, ri_dcoeff(2)%array(sgfi + isgf), sgf(:, sr:er), 1, rhob(:, sr:er, 1), 1)
END IF
!deal with the gradient
IF (do_gga) THEN
DO ipgf = 1, npgf(iset)
start = (ipgf - 1)*ncoset(lmax(iset))
DO ico = ncoset(lmin(iset) - 1) + 1, ncoset(lmax(iset))
DO dir = 1, 3
CALL daxpy(n, ri_sphi(start + ico, sgfi + isgf), dco(:, sr:er, dir, start + ico), &
1, dsgf(:, sr:er, dir), 1)
END DO
END DO !ico
END DO !ipgf
DO dir = 1, 3
CALL daxpy(n, ri_dcoeff(1)%array(sgfi + isgf), dsgf(:, sr:er, dir), 1, &
rho_set%drhoa(dir)%array(:, sr:er, 1), 1)
IF (nspins == 2) THEN
CALL daxpy(n, ri_dcoeff(2)%array(sgfi + isgf), dsgf(:, sr:er, dir), 1, &
rho_set%drhob(dir)%array(:, sr:er, 1), 1)
END IF
END DO
END IF !do_gga
END DO !isgf
DEALLOCATE (co, sgf)
IF (do_gga) DEALLOCATE (dco, dsgf)
END DO !iset
!Treat spin-restricted case (copy alpha into beta)
IF (nspins == 1) THEN
CALL dcopy(n, rhoa(:, sr:er, 1), 1, rhob(:, sr:er, 1), 1)
IF (do_gga) THEN
DO dir = 1, 3
CALL dcopy(n, rho_set%drhoa(dir)%array(:, sr:er, 1), 1, rho_set%drhob(dir)%array(:, sr:er, 1), 1)
END DO
END IF
END IF
CALL timestop(handle)
END SUBROUTINE put_density_on_other_grid
! **************************************************************************************************
!> \brief Computes the norm of the density gradient on the atomic grid
!> \param rho_set ...
!> \param atom_kind ...
!> \param xas_atom_env ...
!> \note GGA is assumed
! **************************************************************************************************
SUBROUTINE compute_norm_drho(rho_set, atom_kind, xas_atom_env)
TYPE(xc_rho_set_type), POINTER :: rho_set
INTEGER, INTENT(IN) :: atom_kind
TYPE(xas_atom_env_type), POINTER :: xas_atom_env
CHARACTER(len=*), PARAMETER :: routineN = 'compute_norm_drho', &
routineP = moduleN//':'//routineN
INTEGER :: dir, ia, ir, n, na, nr, nspins
na = xas_atom_env%grid_atom_set(atom_kind)%grid_atom%ng_sphere
nr = xas_atom_env%grid_atom_set(atom_kind)%grid_atom%nr
n = na*nr
nspins = xas_atom_env%nspins
rho_set%norm_drhoa = 0.0_dp
rho_set%norm_drhob = 0.0_dp
rho_set%norm_drho = 0.0_dp
DO dir = 1, 3
DO ir = 1, nr
DO ia = 1, na
rho_set%norm_drhoa(ia, ir, 1) = rho_set%norm_drhoa(ia, ir, 1) &
+ rho_set%drhoa(dir)%array(ia, ir, 1)**2
END DO !ia
END DO !ir
END DO !dir
rho_set%norm_drhoa = SQRT(rho_set%norm_drhoa)
IF (nspins == 1) THEN
!spin-restricted
CALL dcopy(n, rho_set%norm_drhoa(:, :, 1), 1, rho_set%norm_drhob(:, :, 1), 1)
ELSE
DO dir = 1, 3
DO ir = 1, nr
DO ia = 1, na
rho_set%norm_drhob(ia, ir, 1) = rho_set%norm_drhob(ia, ir, 1) &
+ rho_set%drhob(dir)%array(ia, ir, 1)**2
END DO
END DO
END DO
rho_set%norm_drhob = SQRT(rho_set%norm_drhob)
END IF
DO dir = 1, 3
DO ir = 1, nr
DO ia = 1, na
rho_set%norm_drho(ia, ir, 1) = rho_set%norm_drho(ia, ir, 1) + &
(rho_set%drhoa(dir)%array(ia, ir, 1) + &
rho_set%drhob(dir)%array(ia, ir, 1))**2
END DO
END DO
END DO
rho_set%norm_drho = SQRT(rho_set%norm_drho)
END SUBROUTINE compute_norm_drho
! **************************************************************************************************
!> \brief Precomputes the spherical orbitals of the RI basis on the excited atom grids
!> \param do_gga whether the gradient needs to be computed for GGA or not
!> \param batch_info the parallelization info to complete with so distribution info
!> \param xas_atom_env ...
!> \param qs_env ...
!> \note the functions are split in a purely angular part of size na and a purely radial part of
!> size nr. The full function on the grid can simply be obtained with dgemm and we save space
! **************************************************************************************************
SUBROUTINE precompute_so_dso(do_gga, batch_info, xas_atom_env, qs_env)
LOGICAL, INTENT(IN) :: do_gga
TYPE(batch_info_type) :: batch_info
TYPE(xas_atom_env_type), POINTER :: xas_atom_env
TYPE(qs_environment_type), POINTER :: qs_env
CHARACTER(len=*), PARAMETER :: routineN = 'precompute_so_dso', &
routineP = moduleN//':'//routineN
INTEGER :: bo(2), dir, handle, ikind, ipgf, iset, &
iso, iso_proc, l, maxso, n, na, nkind, &
nr, nset, nso_proc, nsotot, starti
INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf, nsgf_set
INTEGER, DIMENSION(:, :), POINTER :: so_proc_info
REAL(dp), ALLOCATABLE, DIMENSION(:) :: rexp
REAL(dp), DIMENSION(:, :), POINTER :: dgr1, dgr2, ga, gr, slm, zet
REAL(dp), DIMENSION(:, :, :), POINTER :: dga1, dga2, dslm_dxyz
TYPE(cp_para_env_type), POINTER :: para_env
TYPE(grid_atom_type), POINTER :: grid_atom
TYPE(gto_basis_set_type), POINTER :: ri_basis
TYPE(harmonics_atom_type), POINTER :: harmonics
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
NULLIFY (qs_kind_set, harmonics, grid_atom, slm, dslm_dxyz, ri_basis, lmax, lmin, npgf)
NULLIFY (nsgf_set, zet, para_env, so_proc_info)
CALL timeset(routineN, handle)
CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, para_env=para_env)
nkind = SIZE(qs_kind_set)
ALLOCATE (batch_info%so_proc_info(nkind))
ALLOCATE (batch_info%nso_proc(nkind))
ALLOCATE (batch_info%so_bo(2, nkind))
DO ikind = 1, nkind
NULLIFY (xas_atom_env%ga(ikind)%array)
NULLIFY (xas_atom_env%gr(ikind)%array)
NULLIFY (xas_atom_env%dga1(ikind)%array)
NULLIFY (xas_atom_env%dga2(ikind)%array)
NULLIFY (xas_atom_env%dgr1(ikind)%array)
NULLIFY (xas_atom_env%dgr2(ikind)%array)
NULLIFY (batch_info%so_proc_info(ikind)%array)
IF (.NOT. ANY(xas_atom_env%excited_kinds == ikind)) CYCLE
!grid info
harmonics => xas_atom_env%harmonics_atom_set(ikind)%harmonics_atom
grid_atom => xas_atom_env%grid_atom_set(ikind)%grid_atom
na = grid_atom%ng_sphere
nr = grid_atom%nr
n = na*nr
slm => harmonics%slm
dslm_dxyz => harmonics%dslm_dxyz
!basis info
CALL get_qs_kind(qs_kind_set(ikind), basis_set=ri_basis, basis_type="RI_XAS")
CALL get_gto_basis_set(ri_basis, lmax=lmax, npgf=npgf, zet=zet, nset=nset, &
nsgf_set=nsgf_set, lmin=lmin, maxso=maxso)
nsotot = maxso*nset
!we split all so among the processors of the batch
bo = get_limit(nsotot, batch_info%batch_size, batch_info%ipe)
nso_proc = bo(2) - bo(1) + 1
batch_info%so_bo(:, ikind) = bo
batch_info%nso_proc(ikind) = nso_proc
!store info about the so's set, pgf and index
ALLOCATE (batch_info%so_proc_info(ikind)%array(3, nso_proc))
so_proc_info => batch_info%so_proc_info(ikind)%array
so_proc_info = -1 !default is -1 => set so value to zero
DO iset = 1, nset
DO ipgf = 1, npgf(iset)
starti = (iset - 1)*maxso + (ipgf - 1)*nsoset(lmax(iset))
DO iso = nsoset(lmin(iset) - 1) + 1, nsoset(lmax(iset))
!only consider so that are on this proc
IF (starti + iso < bo(1) .OR. starti + iso > bo(2)) CYCLE
iso_proc = starti + iso - bo(1) + 1
so_proc_info(1, iso_proc) = iset
so_proc_info(2, iso_proc) = ipgf
so_proc_info(3, iso_proc) = iso
END DO
END DO
END DO
!Put the gaussians and their gradient as purely angular or radial arrays
ALLOCATE (xas_atom_env%ga(ikind)%array(na, nso_proc))
ALLOCATE (xas_atom_env%gr(ikind)%array(nr, nso_proc))
xas_atom_env%ga(ikind)%array = 0.0_dp; xas_atom_env%gr(ikind)%array = 0.0_dp
IF (do_gga) THEN
ALLOCATE (xas_atom_env%dga1(ikind)%array(na, nso_proc, 3))
ALLOCATE (xas_atom_env%dgr1(ikind)%array(nr, nso_proc))
ALLOCATE (xas_atom_env%dga2(ikind)%array(na, nso_proc, 3))
ALLOCATE (xas_atom_env%dgr2(ikind)%array(nr, nso_proc))
xas_atom_env%dga1(ikind)%array = 0.0_dp; xas_atom_env%dgr1(ikind)%array = 0.0_dp
xas_atom_env%dga2(ikind)%array = 0.0_dp; xas_atom_env%dgr2(ikind)%array = 0.0_dp
END IF
ga => xas_atom_env%ga(ikind)%array
gr => xas_atom_env%gr(ikind)%array
dga1 => xas_atom_env%dga1(ikind)%array
dga2 => xas_atom_env%dga2(ikind)%array
dgr1 => xas_atom_env%dgr1(ikind)%array
dgr2 => xas_atom_env%dgr2(ikind)%array
ALLOCATE (rexp(nr))
DO iso_proc = 1, nso_proc
iset = so_proc_info(1, iso_proc)
ipgf = so_proc_info(2, iso_proc)
iso = so_proc_info(3, iso_proc)
IF (iso < 0) CYCLE
l = indso(1, iso)
!radial part of the gaussian
rexp(1:nr) = EXP(-zet(ipgf, iset)*grid_atom%rad2(1:nr))
gr(1:nr, iso_proc) = grid_atom%rad(1:nr)**l*rexp(1:nr)
!angular part of the gaussian
ga(1:na, iso_proc) = slm(1:na, iso)
IF (do_gga) THEN
!radial part of the gradient => same in all three direction
dgr1(1:nr, iso_proc) = grid_atom%rad(1:nr)**(l - 1)*rexp(1:nr)
dgr2(1:nr, iso_proc) = -2.0_dp*zet(ipgf, iset)*grid_atom%rad(1:nr)**(l + 1)*rexp(1:nr)
!angular part of the gradient
DO dir = 1, 3
dga1(1:na, iso_proc, dir) = dslm_dxyz(dir, 1:na, iso)
dga2(1:na, iso_proc, dir) = harmonics%a(dir, 1:na)*slm(1:na, iso)
END DO
END IF
END DO !iso_proc
END DO !ikind
CALL timestop(handle)
END SUBROUTINE precompute_so_dso
! **************************************************************************************************
!> \brief Integrate the xc kernel as a function of r on the atomic grids for the RI_XAS basis
!> \param int_fxc the global array containing the (P|fxc|Q) integrals, for all spin configurations
!> \param xas_atom_env ...
!> \param xas_tdp_control ...
!> \param qs_env ...
!> \note Note that if closed-shell, alpha-alpha term and beta-beta terms are the same
!> Store the (P|fxc|Q) integrals on the processor they were computed on
!> int_fxc(1)%matrix is alpha-alpha, 2: alpha-beta, 3: beta-beta
! **************************************************************************************************
SUBROUTINE integrate_fxc_atoms(int_fxc, xas_atom_env, xas_tdp_control, qs_env)
TYPE(cp_2d_r_p_type), DIMENSION(:, :), POINTER :: int_fxc
TYPE(xas_atom_env_type), POINTER :: xas_atom_env
TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
TYPE(qs_environment_type), POINTER :: qs_env
CHARACTER(len=*), PARAMETER :: routineN = 'integrate_fxc_atoms', &
routineP = moduleN//':'//routineN
INTEGER :: batch_size, dir, er, ex_bo(2), handle, i, iatom, ibatch, iex, ikind, inb, ipe, &
mepos, na, natom, nb, nb_bo(2), nbatch, nbk, nex_atom, nr, num_pe, sr, subcomm
INTEGER, DIMENSION(2, 3) :: bounds
INTEGER, DIMENSION(:), POINTER :: exat_neighbors
LOGICAL :: do_gga, do_sc, do_sf
TYPE(batch_info_type) :: batch_info
TYPE(cp_1d_r_p_type), DIMENSION(:, :), POINTER :: ri_dcoeff
TYPE(cp_para_env_type), POINTER :: para_env
TYPE(dft_control_type), POINTER :: dft_control
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(section_vals_type), POINTER :: input, xc_functionals
TYPE(xc_derivative_set_type), POINTER :: deriv_set
TYPE(xc_rho_cflags_type) :: needs
TYPE(xc_rho_set_type), POINTER :: rho_set
NULLIFY (particle_set, qs_kind_set, dft_control, deriv_set, para_env, exat_neighbors)
NULLIFY (rho_set, input, xc_functionals)
CALL timeset(routineN, handle)
! Initialize
CALL get_qs_env(qs_env, particle_set=particle_set, qs_kind_set=qs_kind_set, natom=natom, &
dft_control=dft_control, input=input, para_env=para_env)
ALLOCATE (int_fxc(natom, 4))
DO iatom = 1, natom
DO i = 1, 4
NULLIFY (int_fxc(iatom, i)%array)
END DO
END DO
nex_atom = SIZE(xas_atom_env%excited_atoms)
!spin conserving in the general sense here
do_sc = xas_tdp_control%do_spin_cons .OR. xas_tdp_control%do_singlet .OR. xas_tdp_control%do_triplet
do_sf = xas_tdp_control%do_spin_flip
! Get some info on the functionals
xc_functionals => section_vals_get_subs_vals(input, "DFT%XAS_TDP%KERNEL%XC_FUNCTIONAL")
! ask for lsd in any case
needs = xc_functionals_get_needs(xc_functionals, lsd=.TRUE., add_basic_components=.TRUE.)
do_gga = needs%drho_spin !because either LDA or GGA, and the former does not need gradient
! Distribute the excited atoms over batches of processors
! Then, the spherical orbital of the RI basis are distributed among the procs of the batch, making
! the GGA integration very efficient
num_pe = para_env%num_pe
mepos = para_env%mepos
!create a batch_info_type
CALL get_proc_batch_sizes(batch_size, nbatch, nex_atom, num_pe)
!the batch index
ibatch = mepos/batch_size
!the proc index within the batch
ipe = MODULO(mepos, batch_size)
batch_info%batch_size = batch_size
batch_info%nbatch = nbatch
batch_info%ibatch = ibatch
batch_info%ipe = ipe
!create a subcommunicator for this batch
CALL mp_comm_split_direct(para_env%group, subcomm, ibatch)
NULLIFY (batch_info%para_env)
CALL cp_para_env_create(batch_info%para_env, subcomm)
! Precompute the spherical orbital of the RI basis (and maybe their gradient) on the grids of the
! excited atoms. Needed for the GGA integration and to actually put the density on the grid
CALL precompute_so_dso(do_gga, batch_info, xas_atom_env, qs_env)
!distribute the excted atoms over the batches
ex_bo = get_limit(nex_atom, nbatch, ibatch)
! Looping over the excited atoms
DO iex = ex_bo(1), ex_bo(2)
iatom = xas_atom_env%excited_atoms(iex)
ikind = particle_set(iatom)%atomic_kind%kind_number
exat_neighbors => xas_atom_env%exat_neighbors(iex)%array
ri_dcoeff => xas_atom_env%ri_dcoeff(:, :, iex)
! General grid/basis info
na = xas_atom_env%grid_atom_set(ikind)%grid_atom%ng_sphere
nr = xas_atom_env%grid_atom_set(ikind)%grid_atom%nr
! Creating a xc_rho_set to store the density and dset for the kernel
bounds(1:2, 1:3) = 1
bounds(2, 1) = na
bounds(2, 2) = nr
CALL xc_rho_set_create(rho_set=rho_set, local_bounds=bounds, &
rho_cutoff=dft_control%qs_control%eps_rho_rspace, &
drho_cutoff=dft_control%qs_control%eps_rho_rspace)
CALL xc_dset_create(deriv_set, local_bounds=bounds)
! allocate internals of the rho_set
CALL xc_rho_set_atom_update(rho_set, needs, nspins=2, bo=bounds)
! Put the density, and possibly its gradient, on the grid (for this atom)
CALL put_density_on_atomic_grid(rho_set, ri_dcoeff(iatom, :), ikind, &
do_gga, batch_info, xas_atom_env, qs_env)
! Take the neighboring atom contributions to the density (and gradient)
! distribute the grid among the procs (for best load balance)
nb_bo = get_limit(nr, batch_size, ipe)
sr = nb_bo(1); er = nb_bo(2)
DO inb = 1, SIZE(exat_neighbors)
nb = exat_neighbors(inb)
nbk = particle_set(nb)%atomic_kind%kind_number
CALL put_density_on_other_grid(rho_set, ri_dcoeff(nb, :), nb, nbk, iatom, &
ikind, sr, er, do_gga, xas_atom_env, qs_env)
END DO
! make sure contributions from different procs are summed up
CALL mp_sum(rho_set%rhoa, batch_info%para_env%group)
CALL mp_sum(rho_set%rhob, batch_info%para_env%group)
IF (do_gga) THEN
DO dir = 1, 3
CALL mp_sum(rho_set%drhoa(dir)%array, batch_info%para_env%group)
CALL mp_sum(rho_set%drhob(dir)%array, batch_info%para_env%group)
END DO
END IF
! In case of GGA, also need the norm of the density gradient
IF (do_gga) CALL compute_norm_drho(rho_set, ikind, xas_atom_env)
! Compute the required derivatives
CALL xc_functionals_eval(xc_functionals, lsd=.TRUE., rho_set=rho_set, deriv_set=deriv_set, &
deriv_order=2)
!spin-conserving (LDA part)
IF (do_sc) THEN
CALL integrate_sc_fxc(int_fxc, iatom, ikind, deriv_set, xas_atom_env, qs_env)
END IF
!spin-flip (LDA part)
IF (do_sf) THEN
CALL integrate_sf_fxc(int_fxc, iatom, ikind, rho_set, deriv_set, xas_atom_env, qs_env)
END IF
!Gradient correction (note: spin-flip only keeps the lda part, aka ALDA0)
IF (do_gga .AND. do_sc) THEN
CALL integrate_gga_fxc(int_fxc, iatom, ikind, batch_info, rho_set, deriv_set, &
xas_atom_env, qs_env)
END IF
! Clean-up
CALL xc_dset_release(deriv_set)
CALL xc_rho_set_release(rho_set)
END DO !iex
CALL release_batch_info(batch_info)
!Not necessary to sync, but makes sure that any load inbalance is reported here
CALL mp_sync(para_env%group)
CALL timestop(handle)
END SUBROUTINE integrate_fxc_atoms
! **************************************************************************************************
!> \brief Integrate the gradient correction part of the xc kernel on the atomic grid
!> \param int_fxc the array containing the (P|fxc|Q) integrals
!> \param iatom the index of the current excited atom
!> \param ikind the index of the current excited kind
!> \param batch_info how the so are distributed over the processor batch
!> \param rho_set the variable contatinind the density and its gradient
!> \param deriv_set the functional derivatives
!> \param xas_atom_env ...
!> \param qs_env ...
!> \note Ignored in case of pure LDA, added on top of the LDA kernel in case of GGA
! **************************************************************************************************
SUBROUTINE integrate_gga_fxc(int_fxc, iatom, ikind, batch_info, rho_set, deriv_set, &
xas_atom_env, qs_env)
TYPE(cp_2d_r_p_type), DIMENSION(:, :), POINTER :: int_fxc
INTEGER, INTENT(IN) :: iatom, ikind
TYPE(batch_info_type) :: batch_info
TYPE(xc_rho_set_type), POINTER :: rho_set
TYPE(xc_derivative_set_type), POINTER :: deriv_set
TYPE(xas_atom_env_type), POINTER :: xas_atom_env
TYPE(qs_environment_type), POINTER :: qs_env
CHARACTER(len=*), PARAMETER :: routineN = 'integrate_gga_fxc', &
routineP = moduleN//':'//routineN
INTEGER :: bo(2), dir, handle, i, ia, ir, jpgf, &
jset, jso, l, maxso, na, nr, nset, &
nsgf, nsoi, nsotot, startj, ub
INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf
REAL(dp), ALLOCATABLE, DIMENSION(:) :: rexp
REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: int_sgf, res, so, work
REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: dso
REAL(dp), DIMENSION(:, :), POINTER :: dgr1, dgr2, ga, gr, ri_sphi_so, weight, &
zet
REAL(dp), DIMENSION(:, :, :), POINTER :: dga1, dga2
TYPE(cp_2d_r_p_type), ALLOCATABLE, DIMENSION(:) :: int_so, vxc
TYPE(cp_3d_r_p_type), ALLOCATABLE, DIMENSION(:) :: vxg
TYPE(cp_para_env_type), POINTER :: para_env
TYPE(grid_atom_type), POINTER :: grid_atom
TYPE(gto_basis_set_type), POINTER :: ri_basis
TYPE(harmonics_atom_type), POINTER :: harmonics
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
NULLIFY (grid_atom, ri_basis, qs_kind_set, ga, gr, dgr1, dgr2, lmax, lmin, npgf)
NULLIFY (weight, ri_sphi_so, dga1, dga2, para_env, harmonics, zet)
!Strategy: we need to compute