hfx_pair_list_methods.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 Routines for optimizing load balance between processes in HFX calculations
!> \par History
!>      04.2008 created [Manuel Guidon]
!>      11.2019 fixed initial value for potential_id (A. Bussy)
!> \author Manuel Guidon
! **************************************************************************************************
MODULE hfx_pair_list_methods
   USE cell_types,                      ONLY: cell_type,&
                                              pbc
   USE hfx_types,                       ONLY: &
        hfx_basis_type, hfx_block_range_type, hfx_cell_type, hfx_pgf_list, hfx_pgf_product_list, &
        hfx_potential_type, hfx_screen_coeff_type, pair_list_type, pair_set_list_type
   USE input_constants,                 ONLY: do_potential_tshpsc,&
                                              do_potential_id,&
                                              do_potential_mix_cl_trunc,&
                                              do_potential_short,&
                                              do_potential_truncated
   USE kinds,                           ONLY: dp
   USE libint_wrapper,                  ONLY: prim_data_f_size
   USE mathconstants,                   ONLY: pi
   USE mp2_types,                       ONLY: pair_list_type_mp2
   USE particle_types,                  ONLY: particle_type
   USE r12_fundamentals,                ONLY: compute_r12_fundamentals
#include "./base/base_uses.f90"
 
   IMPLICIT NONE
   PRIVATE
 
   PUBLIC :: build_pair_list, &
             build_pair_list_mp2, &
             build_pair_list_pgf, &
             build_pgf_product_list, &
             build_atomic_pair_list, &
             pgf_product_list_size, &
             build_pair_list_pbc_pgf, &
             bra_t, allocate_bra, min_image_gap_sq
 
   ! an initial estimate for the size of the product list
   INTEGER, SAVE :: pgf_product_list_size = 128
 
   ! Store screened primitive gaussian function pairs between two shells (A and B)
   TYPE :: bra_t
      REAL(dp), ALLOCATABLE, DIMENSION(:, :)    :: pgf_scr
      INTEGER, ALLOCATABLE, DIMENSION(:, :)     :: pgf_idx, cell_idx
      INTEGER :: cell_cnt = 0, pgf_cnt = 0
   END TYPE bra_t
 
CONTAINS
 
! **************************************************************************************************
!> \brief Allocates all arrays within a Bra_t structure.
!>
!> This subroutine ensures that a given `Bra_t` pointer has allocated memory
!> for storing primitive Gaussian pair indices, screening data, and cell indices.
!> Existing allocations are safely deallocated first.
!>
!> \param[in,out] bra         Pointer to the Bra_t structure to allocate.
!> \param[in]     max_npgf    Maximum number of primitives per shell.
!> \param[in]     max_ncell   Maximum number of neighbor cells.
!>
!> \note The arrays are allocated as follows:
!>       - `pgf_idx(3, max_npgf^2 * max_ncell)`
!>       - `pgf_scr(5, max_npgf^2 * max_ncell)`
!>       - `cell_idx(3, max_ncell)`
!>
! **************************************************************************************************
   SUBROUTINE allocate_bra(bra, max_npgf, max_ncell)
      TYPE(bra_t), POINTER                               :: bra
      INTEGER                                            :: max_npgf, max_ncell
 
      IF (ALLOCATED(bra%pgf_idx)) DEALLOCATE (bra%pgf_idx)  ! TODO add if n_prm changed
      IF (ALLOCATED(bra%pgf_scr)) DEALLOCATE (bra%pgf_scr)  ! TODO add if n_prm changed
      IF (ALLOCATED(bra%cell_idx)) DEALLOCATE (bra%cell_idx) ! TODO add if n_prm changed
      ALLOCATE (bra%pgf_idx(3, max_npgf*max_npgf*max_ncell))
      ALLOCATE (bra%pgf_scr(5, max_npgf*max_npgf*max_ncell))
      ALLOCATE (bra%cell_idx(3, max_ncell))
   END SUBROUTINE allocate_bra
 
! **************************************************************************************************
!> \brief Builds a screened list of primitives from centers A and B, intersecting with another shell
!>
!> This subroutine populates the `Bra_t` structure with indices and screening data
!> for primitive Gaussian pairs (ipgf,jpgf,n3) that pass the Schwarz screening criterion
!>
!> \param[in]  npgfa              Number of primitives in shell A
!> \param[in]  npgfb              Number of primitives in shell B
!> \param[out] bra                Pointer to Bra_t structure to populate
!> \param[in]  screen1            Screening coeffiecients for the AB pair based on most diffuse pgf
!> \param[in]  screen2            Screening coeffiecients for the CD pair based on most diffuse pgf
!> \param[in]  pgf                Screening coeeficents for primitive gaussians
!> \param[in]  log10_pmax         Log10 of maximum integral prefactor
!> \param[in]  log10_eps_schwarz  Log10 of Schwarz screening threshold
!> \param[in]  ra                 Cartesian coordinates of center A
!> \param[in]  rb                 Cartesian coordinates of center B
!> \param[out] nelements          Total number of valid primitives found
!> \param[in]  neighbor_cells     Array of lattice vectors
!> \param[in]  do_periodic        Logical flag for periodic boundary conditions on B
!>
!> \note Loops over all cells, shifting the B shell position accordingly
!>       Only pairs satisfying the screening condition are stored.
!>
!> \note Each valid cell contributes a block of pairs stored contiguously in
!>       bra%pgf_idx, with cell metadata stored in bra%cell_idx
!> \note screening primitive pairs (ipgf, jpgf) using both coarse and fine screening thresholds
! **************************************************************************************************
   SUBROUTINE build_pair_list_pbc_pgf(npgfa, npgfb, bra, screen1, screen2, pgf, &
                                      log10_pmax, log10_eps_schwarz, ra, rb, nelements, &
                                      neighbor_cells, do_periodic)
 
      INTEGER, INTENT(IN)                                :: npgfa, npgfb
      TYPE(bra_t), POINTER                               :: bra
      REAL(dp), INTENT(IN)                               :: screen1(2), screen2(2)
      TYPE(hfx_screen_coeff_type), DIMENSION(:, :), &
         POINTER                                         :: pgf
      REAL(dp), INTENT(IN)                               :: log10_pmax, log10_eps_schwarz, ra(3), &
                                                            rb(3)
      INTEGER, INTENT(OUT)                               :: nelements
      TYPE(hfx_cell_type), DIMENSION(:), POINTER         :: neighbor_cells
      LOGICAL, INTENT(IN)                                :: do_periodic
 
      INTEGER                                            :: cell_cnt, cell_idx, element_cnt, &
                                                            element_idx, element_off, ipgf, jpgf
      REAL(dp)                                           :: ab(3), im_b(3), pgf_max, rab2
 
      element_idx = 0
      element_off = 0
      cell_cnt = 0
      DO cell_idx = 1, SIZE(neighbor_cells)
         ! move B to this cell while keeping A fixed
         ! NOTE rb has already been moved by build_pair_list
         !      so that when cell_r is (000) AB is in the pbc cell
         IF (do_periodic) THEN
            im_b = rb + neighbor_cells(cell_idx)%cell_r(:)
         ELSE
            im_b = rb
         END IF
         ab = ra - im_b
         rab2 = ab(1)**2 + ab(2)**2 + ab(3)**2
         ! First screening on the most diffuse AB gaussians along with the most diffuse CD gaussian
         IF (screen1(1)*rab2 + screen1(2) + screen2(2) + log10_pmax < log10_eps_schwarz) cycle
 
         element_cnt = 0
         DO ipgf = 1, npgfa
            DO jpgf = 1, npgfb
 
               ! Second screening on this actual AB pair and the most diffuse CD from before
               pgf_max = pgf(jpgf, ipgf)%x(1)*rab2 + pgf(jpgf, ipgf)%x(2)
               IF (pgf_max + screen2(2) + log10_pmax < log10_eps_schwarz) cycle
               ! This pair passed screening. Add it to the bra.
               element_idx = element_idx + 1
               bra%pgf_idx(:, element_idx) = [ipgf, jpgf, cell_idx]
               bra%pgf_scr(1, element_idx) = pgf_max
               element_cnt = element_cnt + 1
            END DO
         END DO
 
         ! If this cell produced any pair, add this cell to the bra
         IF (element_cnt == 0) cycle
         cell_cnt = cell_cnt + 1
         bra%cell_idx(:, cell_cnt) = [cell_idx, element_cnt, element_off]
         element_off = element_off + element_cnt
 
      END DO ! cell
      bra%cell_cnt = cell_cnt
      nelements = element_idx
      bra%pgf_cnt = nelements
 
   END SUBROUTINE build_pair_list_pbc_pgf
 
! **************************************************************************************************
!> \brief ...
!> \param list1 ...
!> \param list2 ...
!> \param product_list ...
!> \param nproducts ...
!> \param log10_pmax ...
!> \param log10_eps_schwarz ...
!> \param neighbor_cells ...
!> \param cell ...
!> \param potential_parameter ...
!> \param m_max ...
!> \param do_periodic ...
! **************************************************************************************************
   SUBROUTINE build_pgf_product_list(list1, list2, product_list, nproducts, &
                                     log10_pmax, log10_eps_schwarz, neighbor_cells, &
                                     cell, potential_parameter, m_max, do_periodic)
 
      TYPE(hfx_pgf_list)                                 :: list1, list2
      TYPE(hfx_pgf_product_list), ALLOCATABLE, &
         DIMENSION(:), INTENT(INOUT)                     :: product_list
      INTEGER, INTENT(OUT)                               :: nproducts
      REAL(dp), INTENT(IN)                               :: log10_pmax, log10_eps_schwarz
      TYPE(hfx_cell_type), DIMENSION(:), POINTER         :: neighbor_cells
      TYPE(cell_type), POINTER                           :: cell
      TYPE(hfx_potential_type)                           :: potential_parameter
      INTEGER, INTENT(IN)                                :: m_max
      LOGICAL, INTENT(IN)                                :: do_periodic
 
      REAL(dp), PARAMETER                                :: r_cutoff_eps = 1.0e-12_dp
 
      INTEGER                                            :: i, j, k, k_max, l, nimages1, nimages2, &
                                                            tmp_i4
      LOGICAL                                            :: has_r_cutoff
      REAL(dp) :: c11(3), den, eta, etainv, fundamentals_work(prim_data_f_size), g(3), &
         min_img_delta, normalization, num, p(3), pgf_max_1, pgf_max_2, pq(3), pq2, prefactor, &
         q(3), r1, r2, r_cutoff_sq, r_cutoff_sq_guard, ra(3), rb(3), rc(3), rcut_eff, rd(3), rho, &
         rhoinv, rpq2, s1234, s1234a, s1234b, shift(3), ssss, t, temp(3), temp_cc(3), temp_dd(3), &
         tmp_d(3), w(3), zeta1, zeta_c, zeta_d, zetapetainv
      TYPE(hfx_pgf_product_list), ALLOCATABLE, &
         DIMENSION(:)                                    :: tmp_product_list
 
      nimages1 = list1%nimages
      nimages2 = list2%nimages
      nproducts = 0
      zeta1 = list1%zetapzetb
      eta = list2%zetapzetb
      etainv = list2%ZetaInv
      zeta_c = list2%zeta
      zeta_d = list2%zetb
      temp_cc = 0.0_dp
      temp_dd = 0.0_dp
      r_cutoff_sq = 0.0_dp
      r_cutoff_sq_guard = 0.0_dp
      SELECT CASE (potential_parameter%potential_type)
      CASE (do_potential_truncated, do_potential_short, do_potential_mix_cl_trunc)
         has_r_cutoff = .true.
         rcut_eff = potential_parameter%cutoff_radius
      CASE (do_potential_tshpsc)
         has_r_cutoff = .true.
         rcut_eff = potential_parameter%cutoff_radius*2.0_dp
      CASE DEFAULT
         has_r_cutoff = .false.
         rcut_eff = 0.0_dp
      END SELECT
      DO i = 1, nimages1
         p = list1%image_list(i)%P
         r1 = list1%image_list(i)%R
         s1234a = list1%image_list(i)%S1234
         pgf_max_1 = list1%image_list(i)%pgf_max
         ra = list1%image_list(i)%ra
         rb = list1%image_list(i)%rb
         DO j = 1, nimages2
            pgf_max_2 = list2%image_list(j)%pgf_max
            IF (pgf_max_1 + pgf_max_2 + log10_pmax < log10_eps_schwarz) cycle
            q = list2%image_list(j)%P
            r2 = list2%image_list(j)%R
            s1234b = list2%image_list(j)%S1234
            rc = list2%image_list(j)%ra
            rd = list2%image_list(j)%rb
 
            IF (has_r_cutoff) THEN
               r_cutoff_sq = (r1 + r2 + rcut_eff)**2
               r_cutoff_sq_guard = r_cutoff_sq*(1.0_dp + r_cutoff_eps)
            END IF
            k_max = SIZE(neighbor_cells)
            IF (do_periodic) THEN
               temp = p - q
               pq = pbc(temp, cell)
               IF (has_r_cutoff .AND. cell%orthorhombic .AND. all(cell%perd == 1) .AND. k_max >= 1) THEN
                  IF (all(neighbor_cells(1)%cell_r == 0.0_dp)) THEN
                     pq2 = pq(1)**2 + pq(2)**2 + pq(3)**2
                     IF (pq2 > r_cutoff_sq_guard) cycle
                     min_img_delta = min_image_gap_sq(pq, cell%hmat)
                     IF (min_img_delta >= 0.0_dp .AND. pq2 + min_img_delta > r_cutoff_sq_guard) k_max = 1
                  END IF
               END IF
               shift = -pq + temp
               temp_cc = rc + shift
               temp_dd = rd + shift
            END IF
 
            zetapetainv = zeta1 + eta
            zetapetainv = 1.0_dp/zetapetainv
            rho = zeta1*eta*zetapetainv
            rhoinv = 1.0_dp/rho
            s1234 = exp(s1234a + s1234b)
 
            DO k = 1, k_max
               IF (do_periodic) THEN
                  IF (has_r_cutoff .AND. k_max > 1) THEN
                     IF ((pq(1) - neighbor_cells(k)%cell_r(1))**2 + &
                         (pq(2) - neighbor_cells(k)%cell_r(2))**2 + &
                         (pq(3) - neighbor_cells(k)%cell_r(3))**2 > r_cutoff_sq_guard) cycle
                  END IF
                  c11 = temp_cc + neighbor_cells(k)%cell_r(:)
                  tmp_d = temp_dd + neighbor_cells(k)%cell_r(:)
               ELSE
                  c11 = rc
                  tmp_d = rd
               END IF
               q = (zeta_c*c11 + zeta_d*tmp_d)*etainv
               rpq2 = (p(1) - q(1))**2 + (p(2) - q(2))**2 + (p(3) - q(3))**2
               IF (has_r_cutoff) THEN
                  IF (rpq2 > r_cutoff_sq) cycle
               END IF
               nproducts = nproducts + 1
 
               ! allocate size as needed,
               ! updating the global size estimate to make this a rare event in longer simulations
               IF (nproducts > SIZE(product_list)) THEN
!$OMP              ATOMIC READ
                  tmp_i4 = pgf_product_list_size
                  tmp_i4 = max(pgf_product_list_size, (3*nproducts + 1)/2)
!$OMP              ATOMIC WRITE
                  pgf_product_list_size = tmp_i4
                  ALLOCATE (tmp_product_list(SIZE(product_list)))
                  tmp_product_list(:) = product_list
                  DEALLOCATE (product_list)
                  ALLOCATE (product_list(tmp_i4))
                  product_list(1:SIZE(tmp_product_list)) = tmp_product_list
                  DEALLOCATE (tmp_product_list)
               END IF
 
               t = rho*rpq2
               IF (potential_parameter%potential_type == do_potential_id) THEN
                  num = list1%zeta*list1%zetb
                  den = list1%zeta + list1%zetb
                  ssss = -num/den*sum((ra - rb)**2)
 
                  num = den*zeta_c
                  den = den + zeta_c
                  ssss = ssss - num/den*sum((p - rc)**2)
 
                  g(:) = (list1%zeta*ra(:) + list1%zetb*rb(:) + zeta_c*rc(:))/den
                  num = den*zeta_d
                  den = den + zeta_d
                  ssss = ssss - num/den*sum((g - rd)**2)
 
                  product_list(nproducts)%Fm(:) = exp(ssss)
                  ! The ID seed already contains the four-Gaussian attenuation.
                  normalization = 1.0_dp
                  IF (s1234 > epsilon(0.0_dp)) normalization = 1.0_dp/s1234
               ELSE
                  CALL compute_r12_fundamentals(product_list(nproducts)%Fm, normalization, m_max, &
                                                potential_parameter%potential_type, rho, t, fundamentals_work, &
                                                omega=potential_parameter%omega, &
                                                cutoff_radius=potential_parameter%cutoff_radius, &
                                                scale_coulomb=potential_parameter%scale_coulomb, &
                                                scale_longrange=potential_parameter%scale_longrange, &
                                                scale_gaussian=potential_parameter%scale_gaussian)
               END IF
 
               prefactor = normalization*s1234*sqrt((pi*zetapetainv)**3)
 
               DO l = 1, m_max + 1
                  product_list(nproducts)%Fm(l) = product_list(nproducts)%Fm(l)*prefactor
               END DO
 
               w = (zeta1*p + eta*q)*zetapetainv
               product_list(nproducts)%ra = ra
               product_list(nproducts)%rb = rb
               product_list(nproducts)%rc = c11
               product_list(nproducts)%rd = tmp_d
               product_list(nproducts)%ZetapEtaInv = zetapetainv
               product_list(nproducts)%Rho = rho
               product_list(nproducts)%RhoInv = rhoinv
               product_list(nproducts)%P = p
               product_list(nproducts)%Q = q
               product_list(nproducts)%W = w
               product_list(nproducts)%AB = ra - rb
               product_list(nproducts)%CD = c11 - tmp_d
            END DO
         END DO
      END DO
 
   END SUBROUTINE build_pgf_product_list
 
! **************************************************************************************************
!> \brief ...
!> \param npgfa ...
!> \param npgfb ...
!> \param list ...
!> \param zeta ...
!> \param zetb ...
!> \param screen1 ...
!> \param screen2 ...
!> \param pgf ...
!> \param R_pgf ...
!> \param log10_pmax ...
!> \param log10_eps_schwarz ...
!> \param ra ...
!> \param rb ...
!> \param nelements ...
!> \param neighbor_cells ...
!> \param nimages ...
!> \param do_periodic ...
!> \param cell ...
! **************************************************************************************************
   SUBROUTINE build_pair_list_pgf(npgfa, npgfb, list, zeta, zetb, screen1, screen2, pgf, R_pgf, &
                                  log10_pmax, log10_eps_schwarz, ra, rb, nelements, &
                                  neighbor_cells, nimages, do_periodic, cell)
      INTEGER, INTENT(IN)                                :: npgfa, npgfb
      TYPE(hfx_pgf_list), DIMENSION(npgfa*npgfb)         :: list
      REAL(dp), DIMENSION(1:npgfa), INTENT(IN)           :: zeta
      REAL(dp), DIMENSION(1:npgfb), INTENT(IN)           :: zetb
      REAL(dp), INTENT(IN)                               :: screen1(2), screen2(2)
      TYPE(hfx_screen_coeff_type), DIMENSION(:, :), &
         POINTER                                         :: pgf, r_pgf
      REAL(dp), INTENT(IN)                               :: log10_pmax, log10_eps_schwarz, ra(3), &
                                                            rb(3)
      INTEGER, INTENT(OUT)                               :: nelements
      TYPE(hfx_cell_type), DIMENSION(:), POINTER         :: neighbor_cells
      INTEGER                                            :: nimages(npgfa*npgfb)
      LOGICAL, INTENT(IN)                                :: do_periodic
      TYPE(cell_type), POINTER                           :: cell
 
      REAL(dp), PARAMETER                                :: min_img_eps = 1.0e-12_dp
 
      INTEGER                                            :: element_counter, i, i_max, ipgf, j, jpgf
      REAL(dp)                                           :: ab(3), im_b(3), min_img_delta, pgf_max, &
                                                            rab2, zeta1, zeta_a, zeta_b, zetainv
 
      i_max = SIZE(neighbor_cells)
      IF (do_periodic) THEN
         IF (cell%orthorhombic .AND. screen1(1) <= 0.0_dp .AND. i_max > 1) THEN
            IF (all(cell%perd == 1) .AND. all(neighbor_cells(1)%cell_r == 0.0_dp)) THEN
               ab = ra - rb
               rab2 = ab(1)**2 + ab(2)**2 + ab(3)**2
               min_img_delta = min_image_gap_sq(ab, cell%hmat)
               IF (min_img_delta >= 0.0_dp) THEN
                  IF (screen1(1)*rab2 + screen1(2) + screen2(2) + log10_pmax < log10_eps_schwarz) THEN
                     nelements = 0
                     RETURN
                  END IF
                  IF (screen1(1)*(rab2 + min_img_delta)*(1.0_dp - min_img_eps) + &
                      screen1(2) + screen2(2) + log10_pmax < log10_eps_schwarz) i_max = 1
               END IF
            END IF
         END IF
      ELSE
         i_max = min(i_max, 1)
      END IF
      nimages = 0
      ! ** inner loop may never be reached
      nelements = npgfa*npgfb
      DO i = 1, i_max
         IF (do_periodic) THEN
            im_b = rb + neighbor_cells(i)%cell_r(:)
         ELSE
            im_b = rb
         END IF
         ab = ra - im_b
         rab2 = ab(1)**2 + ab(2)**2 + ab(3)**2
         IF (screen1(1)*rab2 + screen1(2) + screen2(2) + log10_pmax < log10_eps_schwarz) cycle
         element_counter = 0
         DO ipgf = 1, npgfa
            DO jpgf = 1, npgfb
               element_counter = element_counter + 1
               pgf_max = pgf(jpgf, ipgf)%x(1)*rab2 + pgf(jpgf, ipgf)%x(2)
               IF (pgf_max + screen2(2) + log10_pmax < log10_eps_schwarz) THEN
                  cycle
               END IF
               nimages(element_counter) = nimages(element_counter) + 1
               list(element_counter)%image_list(nimages(element_counter))%pgf_max = pgf_max
               list(element_counter)%image_list(nimages(element_counter))%ra = ra
               list(element_counter)%image_list(nimages(element_counter))%rb = im_b
               list(element_counter)%image_list(nimages(element_counter))%rab2 = rab2
 
               zeta_a = zeta(ipgf)
               zeta_b = zetb(jpgf)
               zeta1 = zeta_a + zeta_b
               zetainv = 1.0_dp/zeta1
 
               IF (nimages(element_counter) == 1) THEN
                  list(element_counter)%ipgf = ipgf
                  list(element_counter)%jpgf = jpgf
                  list(element_counter)%zetaInv = zetainv
                  list(element_counter)%zetapzetb = zeta1
                  list(element_counter)%zeta = zeta_a
                  list(element_counter)%zetb = zeta_b
               END IF
 
               list(element_counter)%image_list(nimages(element_counter))%S1234 = (-zeta_a*zeta_b*zetainv*rab2)
               list(element_counter)%image_list(nimages(element_counter))%P = (zeta_a*ra + zeta_b*im_b)*zetainv
               list(element_counter)%image_list(nimages(element_counter))%R = &
                  max(0.0_dp, r_pgf(jpgf, ipgf)%x(1)*rab2 + r_pgf(jpgf, ipgf)%x(2))
               list(element_counter)%image_list(nimages(element_counter))%ra = ra
               list(element_counter)%image_list(nimages(element_counter))%rb = im_b
               list(element_counter)%image_list(nimages(element_counter))%rab2 = rab2
               list(element_counter)%image_list(nimages(element_counter))%bcell = neighbor_cells(i)%cell
            END DO
         END DO
         nelements = max(nelements, element_counter)
      END DO
      DO j = 1, nelements
         list(j)%nimages = nimages(j)
      END DO
      ! ** Remove unused elements
 
      element_counter = 0
      DO j = 1, nelements
         IF (list(j)%nimages == 0) cycle
         element_counter = element_counter + 1
         list(element_counter)%nimages = list(j)%nimages
         list(element_counter)%zetapzetb = list(j)%zetapzetb
         list(element_counter)%ZetaInv = list(j)%ZetaInv
         list(element_counter)%zeta = list(j)%zeta
         list(element_counter)%zetb = list(j)%zetb
         list(element_counter)%ipgf = list(j)%ipgf
         list(element_counter)%jpgf = list(j)%jpgf
         DO i = 1, list(j)%nimages
            list(element_counter)%image_list(i) = list(j)%image_list(i)
         END DO
      END DO
 
      nelements = element_counter
 
   END SUBROUTINE build_pair_list_pgf
 
! **************************************************************************************************
!> \brief ...
!> \param natom ...
!> \param list ...
!> \param set_list ...
!> \param i_start ...
!> \param i_end ...
!> \param j_start ...
!> \param j_end ...
!> \param kind_of ...
!> \param basis_parameter ...
!> \param particle_set ...
!> \param do_periodic ...
!> \param coeffs_set ...
!> \param coeffs_kind ...
!> \param coeffs_kind_max0 ...
!> \param log10_eps_schwarz ...
!> \param cell ...
!> \param pmax_blocks ...
!> \param atomic_pair_list ...
! **************************************************************************************************
   SUBROUTINE build_pair_list(natom, list, set_list, i_start, i_end, j_start, j_end, kind_of, basis_parameter, particle_set, &
                              do_periodic, coeffs_set, coeffs_kind, coeffs_kind_max0, log10_eps_schwarz, cell, &
                              pmax_blocks, atomic_pair_list)
 
      INTEGER, INTENT(IN)                                :: natom
      TYPE(pair_list_type), INTENT(INOUT)                :: list
      TYPE(pair_set_list_type), DIMENSION(:), &
         INTENT(OUT)                                     :: set_list
      INTEGER, INTENT(IN)                                :: i_start, i_end, j_start, j_end, &
                                                            kind_of(*)
      TYPE(hfx_basis_type), DIMENSION(:), INTENT(IN), &
         POINTER                                         :: basis_parameter
      TYPE(particle_type), DIMENSION(:), INTENT(IN), &
         POINTER                                         :: particle_set
      LOGICAL, INTENT(IN)                                :: do_periodic
      TYPE(hfx_screen_coeff_type), &
         DIMENSION(:, :, :, :), INTENT(IN), POINTER      :: coeffs_set
      TYPE(hfx_screen_coeff_type), DIMENSION(:, :), &
         INTENT(IN)                                      :: coeffs_kind
      REAL(kind=dp), INTENT(IN)                          :: coeffs_kind_max0, log10_eps_schwarz
      TYPE(cell_type), POINTER                           :: cell
      REAL(dp), INTENT(IN)                               :: pmax_blocks
      LOGICAL, DIMENSION(natom, natom), INTENT(IN)       :: atomic_pair_list
 
      INTEGER                                            :: iatom, ikind, iset, jatom, jkind, jset, &
                                                            n_element, nset_ij, nseta, nsetb
      REAL(kind=dp)                                      :: rab2
      REAL(kind=dp), DIMENSION(3)                        :: b11, pbc_b, ra, rb, temp
 
      n_element = 0
      nset_ij = 0
 
      DO iatom = i_start, i_end
         DO jatom = j_start, j_end
            IF (atomic_pair_list(jatom, iatom) .EQV. .false.) cycle
 
            ikind = kind_of(iatom)
            nseta = basis_parameter(ikind)%nset
            ra = particle_set(iatom)%r(:)
 
            IF (jatom < iatom) cycle
            jkind = kind_of(jatom)
            nsetb = basis_parameter(jkind)%nset
            rb = particle_set(jatom)%r(:)
 
            IF (do_periodic) THEN
               temp = rb - ra
               pbc_b = pbc(temp, cell)
               b11 = ra + pbc_b
               rab2 = (ra(1) - b11(1))**2 + (ra(2) - b11(2))**2 + (ra(3) - b11(3))**2
            ELSE
               rab2 = (ra(1) - rb(1))**2 + (ra(2) - rb(2))**2 + (ra(3) - rb(3))**2
               b11 = rb ! ra - rb
            END IF
            IF ((coeffs_kind(jkind, ikind)%x(1)*rab2 + &
                 coeffs_kind(jkind, ikind)%x(2)) + coeffs_kind_max0 + pmax_blocks < log10_eps_schwarz) cycle
 
            n_element = n_element + 1
            list%elements(n_element)%pair = [iatom, jatom]
            list%elements(n_element)%kind_pair = [ikind, jkind]
            list%elements(n_element)%r1 = ra
            list%elements(n_element)%r2 = b11
            list%elements(n_element)%dist2 = rab2
            ! build a list of guaranteed overlapping sets
            list%elements(n_element)%set_bounds(1) = nset_ij + 1
            DO iset = 1, nseta
               DO jset = 1, nsetb
                  IF (coeffs_set(jset, iset, jkind, ikind)%x(1)*rab2 + coeffs_set(jset, iset, jkind, ikind)%x(2) + &
                      coeffs_kind_max0 + pmax_blocks < log10_eps_schwarz) cycle
                  nset_ij = nset_ij + 1
                  set_list(nset_ij)%pair = [iset, jset]
               END DO
            END DO
            list%elements(n_element)%set_bounds(2) = nset_ij
         END DO
      END DO
 
      list%n_element = n_element
 
   END SUBROUTINE build_pair_list
 
! **************************************************************************************************
!> \brief ...
!> \param natom ...
!> \param atomic_pair_list ...
!> \param kind_of ...
!> \param basis_parameter ...
!> \param particle_set ...
!> \param do_periodic ...
!> \param coeffs_kind ...
!> \param coeffs_kind_max0 ...
!> \param log10_eps_schwarz ...
!> \param cell ...
!> \param blocks ...
! **************************************************************************************************
   SUBROUTINE build_atomic_pair_list(natom, atomic_pair_list, kind_of, basis_parameter, particle_set, &
                                     do_periodic, coeffs_kind, coeffs_kind_max0, log10_eps_schwarz, cell, &
                                     blocks)
      INTEGER, INTENT(IN)                                :: natom
      LOGICAL, DIMENSION(natom, natom)                   :: atomic_pair_list
      INTEGER, INTENT(IN)                                :: kind_of(*)
      TYPE(hfx_basis_type), DIMENSION(:), INTENT(IN), &
         POINTER                                         :: basis_parameter
      TYPE(particle_type), DIMENSION(:), INTENT(IN), &
         POINTER                                         :: particle_set
      LOGICAL, INTENT(IN)                                :: do_periodic
      TYPE(hfx_screen_coeff_type), DIMENSION(:, :), &
         INTENT(IN)                                      :: coeffs_kind
      REAL(kind=dp), INTENT(IN)                          :: coeffs_kind_max0, log10_eps_schwarz
      TYPE(cell_type), POINTER                           :: cell
      TYPE(hfx_block_range_type), DIMENSION(:), &
         INTENT(IN), POINTER                             :: blocks
 
      INTEGER                                            :: iatom, iatom_end, iatom_start, iblock, &
                                                            ikind, jatom, jatom_end, jatom_start, &
                                                            jblock, jkind, nseta, nsetb
      REAL(kind=dp)                                      :: rab2
      REAL(kind=dp), DIMENSION(3)                        :: b11, pbc_b, ra, rb, temp
 
      atomic_pair_list = .false.
 
      DO iblock = 1, SIZE(blocks)
         iatom_start = blocks(iblock)%istart
         iatom_end = blocks(iblock)%iend
         DO jblock = 1, SIZE(blocks)
            jatom_start = blocks(jblock)%istart
            jatom_end = blocks(jblock)%iend
 
            DO iatom = iatom_start, iatom_end
               ikind = kind_of(iatom)
               nseta = basis_parameter(ikind)%nset
               ra = particle_set(iatom)%r(:)
               DO jatom = jatom_start, jatom_end
                  IF (jatom < iatom) cycle
                  jkind = kind_of(jatom)
                  nsetb = basis_parameter(jkind)%nset
                  rb = particle_set(jatom)%r(:)
 
                  IF (do_periodic) THEN
                     temp = rb - ra
                     pbc_b = pbc(temp, cell)
                     b11 = ra + pbc_b
                     rab2 = (ra(1) - b11(1))**2 + (ra(2) - b11(2))**2 + (ra(3) - b11(3))**2
                  ELSE
                     rab2 = (ra(1) - rb(1))**2 + (ra(2) - rb(2))**2 + (ra(3) - rb(3))**2
                     b11 = rb ! ra - rb
                  END IF
                  IF ((coeffs_kind(jkind, ikind)%x(1)*rab2 + &
                       coeffs_kind(jkind, ikind)%x(2)) + coeffs_kind_max0 < log10_eps_schwarz) cycle
 
                  atomic_pair_list(jatom, iatom) = .true.
                  atomic_pair_list(iatom, jatom) = .true.
               END DO
            END DO
         END DO
      END DO
 
   END SUBROUTINE build_atomic_pair_list
 
! **************************************************************************************************
!> \brief ...
!> \param natom ...
!> \param list ...
!> \param set_list ...
!> \param i_start ...
!> \param i_end ...
!> \param j_start ...
!> \param j_end ...
!> \param kind_of ...
!> \param basis_parameter ...
!> \param particle_set ...
!> \param do_periodic ...
!> \param coeffs_set ...
!> \param coeffs_kind ...
!> \param coeffs_kind_max0 ...
!> \param log10_eps_schwarz ...
!> \param cell ...
!> \param pmax_blocks ...
!> \param atomic_pair_list ...
!> \param skip_atom_symmetry ...
! **************************************************************************************************
   SUBROUTINE build_pair_list_mp2(natom, list, set_list, i_start, i_end, j_start, j_end, kind_of, basis_parameter, particle_set, &
                                  do_periodic, coeffs_set, coeffs_kind, coeffs_kind_max0, log10_eps_schwarz, cell, &
                                  pmax_blocks, atomic_pair_list, skip_atom_symmetry)
 
      INTEGER, INTENT(IN)                                :: natom
      TYPE(pair_list_type_mp2), INTENT(INOUT)            :: list
      TYPE(pair_set_list_type), DIMENSION(:), &
         INTENT(OUT)                                     :: set_list
      INTEGER, INTENT(IN)                                :: i_start, i_end, j_start, j_end, &
                                                            kind_of(*)
      TYPE(hfx_basis_type), DIMENSION(:), INTENT(IN), &
         POINTER                                         :: basis_parameter
      TYPE(particle_type), DIMENSION(:), INTENT(IN), &
         POINTER                                         :: particle_set
      LOGICAL, INTENT(IN)                                :: do_periodic
      TYPE(hfx_screen_coeff_type), &
         DIMENSION(:, :, :, :), INTENT(IN), POINTER      :: coeffs_set
      TYPE(hfx_screen_coeff_type), DIMENSION(:, :), &
         INTENT(IN)                                      :: coeffs_kind
      REAL(kind=dp), INTENT(IN)                          :: coeffs_kind_max0, log10_eps_schwarz
      TYPE(cell_type), POINTER                           :: cell
      REAL(dp), INTENT(IN)                               :: pmax_blocks
      LOGICAL, DIMENSION(natom, natom), INTENT(IN)       :: atomic_pair_list
      LOGICAL, INTENT(IN), OPTIONAL                      :: skip_atom_symmetry
 
      INTEGER                                            :: iatom, ikind, iset, jatom, jkind, jset, &
                                                            n_element, nset_ij, nseta, nsetb
      LOGICAL                                            :: my_skip_atom_symmetry
      REAL(kind=dp)                                      :: rab2
      REAL(kind=dp), DIMENSION(3)                        :: b11, pbc_b, ra, rb, temp
 
      n_element = 0
      nset_ij = 0
 
      my_skip_atom_symmetry = .false.
      IF (PRESENT(skip_atom_symmetry)) my_skip_atom_symmetry = skip_atom_symmetry
 
      DO iatom = i_start, i_end
         DO jatom = j_start, j_end
            IF (atomic_pair_list(jatom, iatom) .EQV. .false.) cycle
 
            ikind = kind_of(iatom)
            nseta = basis_parameter(ikind)%nset
            ra = particle_set(iatom)%r(:)
 
            IF (jatom < iatom .AND. (.NOT. my_skip_atom_symmetry)) cycle
            jkind = kind_of(jatom)
            nsetb = basis_parameter(jkind)%nset
            rb = particle_set(jatom)%r(:)
 
            IF (do_periodic) THEN
               temp = rb - ra
               pbc_b = pbc(temp, cell)
               b11 = ra + pbc_b
               rab2 = (ra(1) - b11(1))**2 + (ra(2) - b11(2))**2 + (ra(3) - b11(3))**2
            ELSE
               rab2 = (ra(1) - rb(1))**2 + (ra(2) - rb(2))**2 + (ra(3) - rb(3))**2
               b11 = rb ! ra - rb
            END IF
            IF ((coeffs_kind(jkind, ikind)%x(1)*rab2 + &
                 coeffs_kind(jkind, ikind)%x(2)) + coeffs_kind_max0 + pmax_blocks < log10_eps_schwarz) cycle
 
            n_element = n_element + 1
            list%elements(n_element)%pair = [iatom, jatom]
            list%elements(n_element)%kind_pair = [ikind, jkind]
            list%elements(n_element)%r1 = ra
            list%elements(n_element)%r2 = b11
            list%elements(n_element)%dist2 = rab2
            ! build a list of guaranteed overlapping sets
            list%elements(n_element)%set_bounds(1) = nset_ij + 1
            DO iset = 1, nseta
               DO jset = 1, nsetb
                  IF (coeffs_set(jset, iset, jkind, ikind)%x(1)*rab2 + coeffs_set(jset, iset, jkind, ikind)%x(2) + &
                      coeffs_kind_max0 + pmax_blocks < log10_eps_schwarz) cycle
                  nset_ij = nset_ij + 1
                  set_list(nset_ij)%pair = [iset, jset]
               END DO
            END DO
            list%elements(n_element)%set_bounds(2) = nset_ij
         END DO
      END DO
 
      list%n_element = n_element
 
   END SUBROUTINE build_pair_list_mp2
 
! **************************************************************************************************
!> \brief Smallest increase of |r|**2 when r is moved to any other periodic image of an
!>        orthorhombic, fully periodic cell, i.e. min over n /= 0 of |r - hmat*n|**2 - |r|**2.
!>        Exact if r is a minimum-image vector (|r(i)| <= |hmat(i,i)|/2); negative otherwise.
!> \param r displacement vector
!> \param hmat cell matrix; its diagonal entries are signed cell-vector components
!> \return the squared-distance gap to the nearest other periodic image
! **************************************************************************************************
   PURE FUNCTION min_image_gap_sq(r, hmat) RESULT(gap)
      REAL(dp), DIMENSION(3), INTENT(IN)                 :: r
      REAL(dp), DIMENSION(3, 3), INTENT(IN)              :: hmat
      REAL(dp)                                           :: gap
 
      REAL(dp), DIMENSION(3)                             :: cell_len
 
      cell_len = [abs(hmat(1, 1)), abs(hmat(2, 2)), abs(hmat(3, 3))]
      gap = minval(cell_len*(cell_len - 2.0_dp*abs(r)))
 
   END FUNCTION min_image_gap_sq
 
END MODULE hfx_pair_list_methods