ri_rs_utils.F

Go to the documentation of this file.

!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2026 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!
 
! **************************************************************************************************
!> \brief Geometry and AO evaluation for real-space resolution of the identity.
! **************************************************************************************************
MODULE ri_rs_utils
   USE atomic_kind_types,               ONLY: get_atomic_kind_set
   USE basis_set_types,                 ONLY: gto_basis_set_p_type,&
                                              gto_basis_set_type
   USE cell_types,                      ONLY: cell_type,&
                                              get_cell,&
                                              pbc
   USE kinds,                           ONLY: dp
   USE orbital_pointers,                ONLY: indco,&
                                              ncoset
   USE particle_types,                  ONLY: particle_type
   USE physcon,                         ONLY: angstrom
   USE qs_kind_types,                   ONLY: get_qs_kind,&
                                              qs_kind_type
   USE ri_rs_types,                     ONLY: ri_rs_env
#include "./base/base_uses.f90"
 
   IMPLICIT NONE
   PRIVATE
 
   PUBLIC :: evaluate_ao_basis_on_points, evaluate_ao_on_points, filter_grid_to_voronoi, &
             get_rirs_cluster_atoms, precompute_ri_rs_radii
 
   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ri_rs_utils'
 
CONTAINS
 
! **************************************************************************************************
!> \brief Compute per-atom AO and RI basis radii from the most diffuse Gaussian
!>        primitive in the AO ("ORB") and RI auxiliary ("RI_AUX") basis sets.
!>        Stores per-atom screening radii in ri_rs and prints a per-kind table.
!>        Radius: r_kind = sqrt(-log(eps)/alpha_min_kind), with eps = eps_filter.
!> \param ri_rs RI-RS geometry references and screening radii
!> \param basis_set_AO AO basis sets, one per atomic kind
!> \param basis_set_RI auxiliary basis sets, one per atomic kind
!> \param eps screening threshold
!> \param unit_nr output unit
! **************************************************************************************************
   SUBROUTINE precompute_ri_rs_radii(ri_rs, basis_set_AO, basis_set_RI, eps, unit_nr)
      TYPE(ri_rs_env), INTENT(INOUT)                     :: ri_rs
      TYPE(gto_basis_set_p_type), INTENT(IN)             :: basis_set_ao(:), basis_set_ri(:)
      REAL(kind=dp), INTENT(IN)                          :: eps
      INTEGER, INTENT(IN)                                :: unit_nr
 
      CHARACTER(LEN=*), PARAMETER :: routinen = 'precompute_ri_rs_radii'
      REAL(kind=dp), PARAMETER :: min_exponent_for_radius = 1.0e-3_dp
 
      INTEGER                                            :: handle, i, iatom, ikind, j, natom, nkind
      INTEGER, ALLOCATABLE, DIMENSION(:)                 :: kind_of
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: alpha_min_ao_kind, alpha_min_ri_kind
      REAL(kind=dp), DIMENSION(:, :), POINTER            :: zet_ao, zet_ri
 
      CALL timeset(routinen, handle)
 
      cpassert(ASSOCIATED(ri_rs%atomic_kind_set))
      cpassert(ASSOCIATED(ri_rs%particle_set))
      cpassert(0.0_dp < eps .AND. eps < 1.0_dp)
      nkind = SIZE(ri_rs%atomic_kind_set)
      natom = SIZE(ri_rs%particle_set)
 
      ALLOCATE (alpha_min_ao_kind(nkind), alpha_min_ri_kind(nkind))
      alpha_min_ao_kind = huge(1.0_dp)
      alpha_min_ri_kind = huge(1.0_dp)
 
      DO ikind = 1, nkind
         zet_ao => basis_set_ao(ikind)%gto_basis_set%zet
         zet_ri => basis_set_ri(ikind)%gto_basis_set%zet
         DO i = 1, SIZE(zet_ao, 1)
            DO j = 1, SIZE(zet_ao, 2)
               IF (zet_ao(i, j) > min_exponent_for_radius) THEN
                  alpha_min_ao_kind(ikind) = min(alpha_min_ao_kind(ikind), zet_ao(i, j))
               END IF
            END DO
         END DO
         DO i = 1, SIZE(zet_ri, 1)
            DO j = 1, SIZE(zet_ri, 2)
               IF (zet_ri(i, j) > min_exponent_for_radius) THEN
                  alpha_min_ri_kind(ikind) = min(alpha_min_ri_kind(ikind), zet_ri(i, j))
               END IF
            END DO
         END DO
      END DO
 
      CALL get_atomic_kind_set(atomic_kind_set=ri_rs%atomic_kind_set, kind_of=kind_of)
 
      ALLOCATE (ri_rs%radius_ao_per_atom(natom))
      ALLOCATE (ri_rs%radius_ri_per_atom(natom))
      DO iatom = 1, natom
         ikind = kind_of(iatom)
         ri_rs%radius_ao_per_atom(iatom) = sqrt(-log(eps)/alpha_min_ao_kind(ikind))
         ri_rs%radius_ri_per_atom(iatom) = sqrt(-log(eps)/alpha_min_ri_kind(ikind))
      END DO
 
      IF (unit_nr > 0) THEN
         WRITE (unit_nr, '(T2,A)') 'RI-RS basis radii (Å):'
         WRITE (unit_nr, '(T4,A6,2X,A4,2A14)') 'Kind', 'Elem', 'r_AO (Å)', 'r_RI (Å)'
         DO ikind = 1, nkind
            WRITE (unit_nr, '(T4,I6,2X,A4,2F14.4)') &
               ikind, &
               ri_rs%atomic_kind_set(ikind)%element_symbol, &
               sqrt(-log(eps)/alpha_min_ao_kind(ikind))*angstrom, &
               sqrt(-log(eps)/alpha_min_ri_kind(ikind))*angstrom
         END DO
         WRITE (unit_nr, '(A)') ' '
      END IF
 
      DEALLOCATE (alpha_min_ao_kind, alpha_min_ri_kind, kind_of)
 
      CALL timestop(handle)
 
   END SUBROUTINE precompute_ri_rs_radii
 
!> \brief Retain source-grid points inside the Voronoi volume of one atom.
!> \param points Source-grid coordinates on entry and Voronoi-filtered coordinates on return.
!> \param icenter_atom Atom whose Voronoi volume is retained.
!> \param particle_set Molecular atom positions.
!> \param mask Optional membership mask; if present, leave points unchanged.
!> \param atom_indices Optional subset of nuclei defining the Voronoi partition.
! **************************************************************************************************
   SUBROUTINE filter_grid_to_voronoi(points, icenter_atom, particle_set, mask, atom_indices)
      REAL(kind=dp), ALLOCATABLE, INTENT(INOUT)          :: points(:, :)
      INTEGER, INTENT(IN)                                :: icenter_atom
      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
      LOGICAL, INTENT(OUT), OPTIONAL                     :: mask(:)
      INTEGER, INTENT(IN), OPTIONAL                      :: atom_indices(:)
 
      INTEGER                                            :: iatom, iatom_index, ipoint, n_atoms, &
                                                            n_keep
      LOGICAL, ALLOCATABLE                               :: keep(:)
      REAL(kind=dp)                                      :: displacement(3), distance2, max_radius2, &
                                                            tolerance
      REAL(kind=dp), ALLOCATABLE                         :: filtered_points(:, :)
 
      ALLOCATE (keep(SIZE(points, 2)), source=.true.)
      max_radius2 = maxval(sum(points**2, dim=1))
      n_atoms = SIZE(particle_set)
      IF (PRESENT(atom_indices)) n_atoms = SIZE(atom_indices)
      DO iatom_index = 1, n_atoms
         IF (PRESENT(atom_indices)) THEN
            iatom = atom_indices(iatom_index)
            cpassert(iatom >= 1 .AND. iatom <= SIZE(particle_set))
         ELSE
            iatom = iatom_index
         END IF
         IF (iatom == icenter_atom) cycle
         displacement(:) = particle_set(iatom)%r - particle_set(icenter_atom)%r
         distance2 = sum(displacement**2)
         ! |R_B-R_A| > 2 max_l|r_l-R_A| cannot cut this finite point set.
         IF (distance2 > 4.0_dp*max_radius2) cycle
         tolerance = 32.0_dp*epsilon(1.0_dp)*max(1.0_dp, distance2)
         DO ipoint = 1, SIZE(points, 2)
            IF (.NOT. keep(ipoint)) cycle
            ! |r_l-R_B|² - |r_l-R_A|² = |R_B-R_A|² - 2(r_l-R_A)·(R_B-R_A).
            IF (2.0_dp*dot_product(points(:, ipoint), displacement) > distance2 + tolerance) THEN
               keep(ipoint) = .false.
            ELSE IF (iatom < icenter_atom) THEN
               IF (abs(2.0_dp*dot_product(points(:, ipoint), displacement) - distance2) <= tolerance) THEN
                  keep(ipoint) = .false.
               END IF
            END IF
         END DO
      END DO
      IF (PRESENT(mask)) THEN
         cpassert(SIZE(mask) == SIZE(keep))
         mask(:) = keep
         RETURN
      END IF
      n_keep = count(keep)
      ALLOCATE (filtered_points(3, n_keep))
      IF (n_keep > 0) filtered_points(:, :) = reshape(pack(points, spread(keep, 1, 3)), [3, n_keep])
      CALL move_alloc(filtered_points, points)
   END SUBROUTINE filter_grid_to_voronoi
 
! **************************************************************************************************
!> \brief Form C_A from nuclei within the specified radius, in global atom order.
!> \param particle_set Molecular nuclei.
!> \param cell Simulation cell.
!> \param icenter_atom Central atom A.
!> \param radius Cluster radius.
!> \param atom_indices Cluster atom indices including A.
! **************************************************************************************************
   SUBROUTINE get_rirs_cluster_atoms(particle_set, cell, icenter_atom, radius, atom_indices)
      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
      TYPE(cell_type), POINTER                           :: cell
      INTEGER, INTENT(IN)                                :: icenter_atom
      REAL(kind=dp), INTENT(IN)                          :: radius
      INTEGER, ALLOCATABLE, INTENT(OUT)                  :: atom_indices(:)
 
      INTEGER                                            :: iatom, ncluster
      INTEGER, ALLOCATABLE                               :: work(:)
 
      ALLOCATE (work(SIZE(particle_set)))
      ncluster = 0
      DO iatom = 1, SIZE(particle_set)
         IF (sum(pbc(particle_set(iatom)%r - particle_set(icenter_atom)%r, cell)**2) <= radius**2) THEN
            ncluster = ncluster + 1
            work(ncluster) = iatom
         END IF
      END DO
      ALLOCATE (atom_indices(ncluster), source=work(:ncluster))
   END SUBROUTINE get_rirs_cluster_atoms
 
! **************************************************************************************************
!> \brief Evaluate contracted spherical AOs, and optionally their grid-coordinate derivatives.
!>
!> For atom A and displacement d(alpha)=r_l(alpha)-R_A(alpha), the value is
!>
!>   ϕ_μ(r_l) = Σ_(p,c) S_(pc,μ) d_x^lx d_y^ly d_z^lz exp(-ζ_p |d|^2).
!>
!> The optional derivative is evaluated analytically as
!>
!>   d ϕ_μ(r_l)/d r_(l,α)
!>     = sum_(p,c) S_(pc,mu) exp(-zeta_p |d|^2)
!>       [d polynomial_c/d d_alpha - 2 zeta_p d_alpha polynomial_c].
!>
!> \param phi AO values, accumulated into phi(l,mu).
!> \param grid_points Cartesian grid points, indexed (alpha,l).
!> \param iatom Source atom whose contracted AOs are evaluated.
!> \param particle_set Molecular particles.
!> \param qs_kind_set Quickstep atomic kinds containing the ORB bases.
!> \param cell Simulation cell used for the minimum-image displacement.
!> \param dphi Optional AO derivatives, accumulated into dphi(alpha,l,mu).
!> \param cutoff_squared Optional squared AO cutoff radius.
! **************************************************************************************************
   SUBROUTINE evaluate_ao_on_points(phi, grid_points, iatom, particle_set, qs_kind_set, cell, &
                                    dphi, cutoff_squared)
      REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT)      :: phi
      REAL(kind=dp), DIMENSION(:, :), INTENT(IN)         :: grid_points
      INTEGER, INTENT(IN)                                :: iatom
      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
      TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
      TYPE(cell_type), POINTER                           :: cell
      REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT), &
         OPTIONAL                                        :: dphi
      REAL(kind=dp), INTENT(IN), OPTIONAL                :: cutoff_squared
 
      CHARACTER(len=*), PARAMETER :: routinen = 'evaluate_ao_on_points'
 
      INTEGER                                            :: handle, ikind
      TYPE(gto_basis_set_type), POINTER                  :: basis
 
      CALL timeset(routinen, handle)
      cpassert(SIZE(grid_points, 1) == 3)
      cpassert(SIZE(phi, 1) == SIZE(grid_points, 2))
      IF (PRESENT(dphi)) THEN
         cpassert(SIZE(dphi, 1) == 3)
         cpassert(SIZE(dphi, 2) == SIZE(phi, 1))
         cpassert(SIZE(dphi, 3) == SIZE(phi, 2))
      END IF
 
      ikind = particle_set(iatom)%atomic_kind%kind_number
      CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis, basis_type='ORB')
      IF (.NOT. ASSOCIATED(basis)) THEN
         CALL timestop(handle)
         RETURN
      END IF
      CALL evaluate_ao_basis_on_points(phi, grid_points, basis, &
                                       particle_set(iatom)%r, cell, dphi, cutoff_squared)
      CALL timestop(handle)
   END SUBROUTINE evaluate_ao_on_points
 
! **************************************************************************************************
!> \brief Evaluate one explicitly supplied contracted Gaussian basis on Cartesian points.
!> \param phi AO values, accumulated into phi(l,mu).
!> \param grid_points Cartesian grid points, indexed (alpha,l).
!> \param basis Contracted Gaussian basis to evaluate.
!> \param source_position Centre of the basis.
!> \param cell Simulation cell defining minimum-image or image-summed displacements.
!> \param dphi Optional AO derivatives, accumulated into dphi(alpha,l,mu).
!> \param cutoff_squared Optional squared AO cutoff radius.
!> \param periodic_images Sum images -1:1 along periodic axes and only image 0 along other axes.
! **************************************************************************************************
   SUBROUTINE evaluate_ao_basis_on_points(phi, grid_points, basis, source_position, cell, &
                                          dphi, cutoff_squared, periodic_images)
      REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT)      :: phi
      REAL(kind=dp), DIMENSION(:, :), INTENT(IN)         :: grid_points
      TYPE(gto_basis_set_type), POINTER                  :: basis
      REAL(kind=dp), DIMENSION(3), INTENT(IN)            :: source_position
      TYPE(cell_type), POINTER                           :: cell
      REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT), &
         OPTIONAL                                        :: dphi
      REAL(kind=dp), INTENT(IN), OPTIONAL                :: cutoff_squared
      LOGICAL, INTENT(IN), OPTIONAL                      :: periodic_images
 
      INTEGER :: first_sgf, ialpha, ico, iend_co, image_extent(3), ipgf, ipoint, irow, iset, isgf, &
         ishell, istart_co, ix, iy, iz, l, last_sgf, lx, ly, lz, n_cart_total
      LOGICAL                                            :: sum_images
      REAL(kind=dp) :: exponent, exponential, hmat(3, 3), polynomial, polynomial_derivative(3), &
         radius2, relative(3), relative_raw(3), weight
 
      cpassert(ASSOCIATED(basis))
      cpassert(SIZE(grid_points, 1) == 3)
      cpassert(SIZE(phi, 1) == SIZE(grid_points, 2))
      IF (PRESENT(dphi)) THEN
         cpassert(SIZE(dphi, 1) == 3)
         cpassert(SIZE(dphi, 2) == SIZE(phi, 1))
         cpassert(SIZE(dphi, 3) == SIZE(phi, 2))
      END IF
 
      sum_images = .false.
      IF (PRESENT(periodic_images)) sum_images = periodic_images
      image_extent = 0
      hmat = 0.0_dp
      IF (sum_images) THEN
         image_extent = cell%perd
         CALL get_cell(cell, h=hmat)
      END IF
 
      !$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(STATIC) &
      !$OMP SHARED(phi, dphi, grid_points, cell, basis, source_position, cutoff_squared, ncoset, indco) &
      !$OMP SHARED(sum_images, image_extent, hmat) &
      !$OMP PRIVATE(ipoint, relative, radius2, iset, n_cart_total, ishell, l, istart_co, &
      !$OMP         iend_co, first_sgf, last_sgf, ipgf, exponent, exponential, isgf, ico, &
      !$OMP         irow, weight, lx, ly, lz, polynomial, polynomial_derivative, ialpha, &
      !$OMP         relative_raw, ix, iy, iz)
      DO ipoint = 1, SIZE(grid_points, 2)
         relative_raw = grid_points(:, ipoint) - source_position
         IF (.NOT. sum_images) relative_raw = pbc(relative_raw, cell)
         DO ix = -image_extent(1), image_extent(1)
         DO iy = -image_extent(2), image_extent(2)
         DO iz = -image_extent(3), image_extent(3)
            relative = relative_raw
            IF (sum_images) relative = relative_raw - matmul(hmat, real([ix, iy, iz], dp))
            radius2 = dot_product(relative, relative)
            IF (PRESENT(cutoff_squared)) THEN
               IF (radius2 > cutoff_squared) cycle
            END IF
 
            DO iset = 1, basis%nset
               n_cart_total = ncoset(basis%lmax(iset))
               DO ishell = 1, basis%nshell(iset)
                  l = basis%l(ishell, iset)
                  istart_co = ncoset(l - 1) + 1
                  iend_co = ncoset(l)
                  first_sgf = basis%first_sgf(ishell, iset)
                  last_sgf = basis%last_sgf(ishell, iset)
                  DO ipgf = 1, basis%npgf(iset)
                     exponent = basis%zet(ipgf, iset)
                     exponential = exp(-exponent*radius2)
                     DO isgf = first_sgf, last_sgf
                        DO ico = istart_co, iend_co
                           irow = (ipgf - 1)*n_cart_total + ico
                           weight = basis%sphi(irow, isgf)
                           lx = indco(1, ico)
                           ly = indco(2, ico)
                           lz = indco(3, ico)
                           polynomial = relative(1)**lx*relative(2)**ly*relative(3)**lz
                           phi(ipoint, isgf) = phi(ipoint, isgf) + weight*polynomial*exponential
 
                           IF (PRESENT(dphi)) THEN
                              polynomial_derivative = 0.0_dp
                              IF (lx > 0) polynomial_derivative(1) = real(lx, dp)*relative(1)**(lx - 1)* &
                                                                     relative(2)**ly*relative(3)**lz
                              IF (ly > 0) polynomial_derivative(2) = real(ly, dp)*relative(1)**lx* &
                                                                     relative(2)**(ly - 1)*relative(3)**lz
                              IF (lz > 0) polynomial_derivative(3) = real(lz, dp)*relative(1)**lx* &
                                                                     relative(2)**ly*relative(3)**(lz - 1)
                              DO ialpha = 1, 3
                                 dphi(ialpha, ipoint, isgf) = dphi(ialpha, ipoint, isgf) + weight*exponential* &
                                                              (polynomial_derivative(ialpha) - &
                                                               2.0_dp*exponent*relative(ialpha)*polynomial)
                              END DO
                           END IF
                        END DO
                     END DO
                  END DO
               END DO
            END DO
         END DO
         END DO
         END DO
      END DO
      !$OMP END PARALLEL DO
   END SUBROUTINE evaluate_ao_basis_on_points
 
END MODULE ri_rs_utils