eri_mme_error_control.F

Go to the documentation of this file.

!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!
 
! **************************************************************************************************
!> \brief Methods aiming for error estimate and automatic cutoff calibration.
!>        integrals.
!> \par History
!>       2015 09 created
!> \author Patrick Seewald
! **************************************************************************************************
 
MODULE eri_mme_error_control
   USE ao_util,                         ONLY: exp_radius
   USE eri_mme_gaussian,                ONLY: get_minimax_coeff_v_gspace,&
                                              hermite_gauss_norm
   USE eri_mme_lattice_summation,       ONLY: pgf_sum_2c_gspace_1d_deltal
   USE kinds,                           ONLY: dp
   USE mathconstants,                   ONLY: pi,&
                                              twopi
   USE message_passing,                 ONLY: mp_para_env_type
#include "../base/base_uses.f90"
 
   IMPLICIT NONE
 
   PRIVATE
 
   LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .false.
 
   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'eri_mme_error_control'
 
   PUBLIC :: calibrate_cutoff, cutoff_minimax_error, minimax_error, cutoff_error
CONTAINS
 
! **************************************************************************************************
!> \brief Find optimal cutoff minimizing errors due to minimax approximation and
!>        due to finite cutoff using bisection on the difference of the errors
!> \param hmat ...
!> \param h_inv ...
!> \param G_min ...
!> \param vol ...
!> \param zet_min   Minimum exponent
!> \param l_mm      Total ang. mom. quantum number
!> \param zet_max     Max. exponents to estimate cutoff error
!> \param l_max_zet       Max. total ang. mom. quantum numbers to estimate cutoff error
!> \param n_minimax Number of terms in minimax approximation
!> \param cutoff_l  Initial guess of lower bound for cutoff
!> \param cutoff_r  Initial guess of upper bound for cutoff
!> \param tol       Tolerance (cutoff precision)
!> \param delta     to modify initial guess interval
!> \param cutoff    Best cutoff
!> \param err_mm    Minimax error
!> \param err_c     Cutoff error
!> \param C_mm      Scaling constant to generalize AM-GM upper bound estimate to
!>                  minimax approx.
!> \param para_env ...
!> \param print_calib ...
!> \param unit_nr ...
! **************************************************************************************************
   SUBROUTINE calibrate_cutoff(hmat, h_inv, G_min, vol, zet_min, l_mm, zet_max, l_max_zet, &
                               n_minimax, cutoff_l, cutoff_r, tol, delta, &
                               cutoff, err_mm, err_c, C_mm, para_env, print_calib, unit_nr)
      REAL(kind=dp), DIMENSION(3, 3), INTENT(IN)         :: hmat, h_inv
      REAL(kind=dp), INTENT(IN)                          :: g_min
      REAL(kind=dp)                                      :: vol
      REAL(kind=dp), INTENT(IN)                          :: zet_min
      INTEGER, INTENT(IN)                                :: l_mm
      REAL(kind=dp), INTENT(IN)                          :: zet_max
      INTEGER, INTENT(IN)                                :: l_max_zet, n_minimax
      REAL(kind=dp), INTENT(IN)                          :: cutoff_l, cutoff_r, tol, delta
      REAL(kind=dp), INTENT(OUT)                         :: cutoff, err_mm, err_c, c_mm
      TYPE(mp_para_env_type), INTENT(IN)                 :: para_env
      LOGICAL, INTENT(IN)                                :: print_calib
      INTEGER, INTENT(IN)                                :: unit_nr
 
      INTEGER                                            :: i, iter1, iter2, max_iter
      LOGICAL                                            :: do_print, valid_initial
      REAL(kind=dp)                                      :: cutoff_mid, delta_c_mid, delta_mm_mid
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: minimax_aw
      REAL(kind=dp), DIMENSION(2)                        :: cutoff_lr, delta_c, delta_mm
 
      do_print = unit_nr > 0 .AND. print_calib
      IF (do_print) THEN
         WRITE (unit_nr, '(/T2, A)') "ERI_MME| Basis set parameters for estimating minimax error"
         WRITE (unit_nr, '(T2, A, T67, ES12.2, 1X, I1)') "ERI_MME|   exp, l:", zet_min, l_mm
         WRITE (unit_nr, '(T2, A)') "ERI_MME| Basis set parameters for estimating cutoff error"
         WRITE (unit_nr, '(T2, A, T67, ES12.2, 1X, I1)') "ERI_MME|   exp, l:", zet_max, l_max_zet
      END IF
 
      max_iter = 100
 
      IF ((cutoff_r - cutoff_l)/(0.5_dp*(cutoff_r + cutoff_l)) .LE. tol) &
         CALL cp_abort(__location__, "difference of boundaries for cutoff "// &
                       "(MAX - MIN) must be greater than cutoff precision.")
 
      IF ((delta .GE. 1.0_dp) .OR. (delta .LE. 0.0_dp)) &
         CALL cp_abort(__location__, &
                       "relative delta to modify initial cutoff interval (DELTA) must be in (0, 1)")
 
      cutoff_lr(1) = cutoff_l
      cutoff_lr(2) = cutoff_r
 
      ALLOCATE (minimax_aw(2*n_minimax))
 
      IF (do_print) THEN
         WRITE (unit_nr, '(/T2, A)') "ERI_MME| Calibrating cutoff by bisecting error(minimax) - error(cutoff)"
         WRITE (unit_nr, '(T2, A, T72, ES9.2)') "ERI_MME| Rel. cutoff precision", tol
         WRITE (unit_nr, '(T2, A, T77, F4.1)') "ERI_MME| Rel. cutoff delta to modify initial interval", delta
      END IF
 
      ! 1) find valid initial values for bisection
      DO iter1 = 1, max_iter + 1
         IF (iter1 .GT. max_iter) &
            CALL cp_abort(__location__, &
                          "Maximum number of iterations in bisection to determine initial "// &
                          "cutoff interval has been exceeded.")
 
         cutoff_lr(1) = max(cutoff_lr(1), 0.5_dp*g_min**2)
         ! approx.) is hit
 
         DO i = 1, 2
            CALL cutoff_minimax_error(cutoff_lr(i), hmat, h_inv, vol, g_min, zet_min, l_mm, zet_max, l_max_zet, &
                                      n_minimax, minimax_aw, delta_mm(i), delta_c(i), c_mm, para_env)
         END DO
 
         valid_initial = .true.
         IF ((delta_mm(1) - delta_c(1)) .GT. 0) THEN
            cutoff_lr(1) = cutoff_lr(1)*(1.0_dp - abs(delta))
            valid_initial = .false.
         END IF
         IF ((delta_mm(2) - delta_c(2)) .LT. 0) THEN
            cutoff_lr(2) = cutoff_lr(2)*(1.0_dp + abs(delta))
            valid_initial = .false.
         END IF
 
         IF (valid_initial) EXIT
      END DO
 
      ! 2) bisection to find cutoff s.t. err_minimax(cutoff) - err_cutoff(cutoff) = 0
      IF (do_print) WRITE (unit_nr, '(/T2, A)') &
         "ERI_MME| Step, cutoff (min, max, mid), err(minimax), err(cutoff), err diff"
 
      DO iter2 = 1, max_iter + 1
         IF (iter2 .GT. max_iter) &
            CALL cp_abort(__location__, &
                          "Maximum number of iterations in bisection to determine cutoff has been exceeded")
 
         cutoff_mid = 0.5_dp*(cutoff_lr(1) + cutoff_lr(2))
         CALL cutoff_minimax_error(cutoff_mid, hmat, h_inv, vol, g_min, zet_min, l_mm, zet_max, l_max_zet, &
                                   n_minimax, minimax_aw, delta_mm_mid, delta_c_mid, c_mm, para_env)
         IF (do_print) WRITE (unit_nr, '(T11, I2, F11.1, F11.1, F11.1, 3X, ES9.2, 3X, ES9.2, 3X, ES9.2)') &
            iter2, cutoff_lr(1), cutoff_lr(2), cutoff_mid, &
            delta_mm_mid, delta_c_mid, delta_mm_mid - delta_c_mid
 
         IF ((cutoff_lr(2) - cutoff_lr(1))/cutoff_mid .LT. tol) EXIT
         IF (delta_mm_mid - delta_c_mid .GT. 0) THEN
            cutoff_lr(2) = cutoff_mid
            delta_mm(2) = delta_mm_mid
            delta_c(2) = delta_c_mid
         ELSE
            cutoff_lr(1) = cutoff_mid
            delta_mm(1) = delta_mm_mid
            delta_c(1) = delta_c_mid
         END IF
      END DO
      err_mm = delta_mm_mid
      err_c = delta_c_mid
      cutoff = cutoff_mid
 
      IF (do_print) THEN
         WRITE (unit_nr, '(/T2, A)') "ERI_MME| Cutoff calibration number of steps:"
         WRITE (unit_nr, '(T2, A, T79, I2)') "ERI_MME|   Steps for initial interval", iter1 - 1
         WRITE (unit_nr, '(T2, A, T79, I2/)') "ERI_MME|   Bisection iteration steps", iter2 - 1
      END IF
 
   END SUBROUTINE calibrate_cutoff
 
! **************************************************************************************************
!> \brief Compute upper bounds for the errors of 2-center ERI's (P|P) due
!>        to minimax approximation and due to finite cutoff, where P is a
!>        normalized Hermite Gaussian.
!> \param cutoff ...
!> \param hmat ...
!> \param h_inv ...
!> \param vol ...
!> \param G_min ...
!> \param zet_min     Exponent of P to estimate minimax error
!> \param l_mm       total ang. mom. quantum number of P to estimate minimax error
!> \param zet_max   Max. exponents of P to estimate cutoff error
!> \param l_max_zet     Max. total ang. mom. quantum numbers of P to estimate cutoff error
!> \param n_minimax  Number of terms in minimax approximation
!> \param minimax_aw Minimax coefficients
!> \param err_mm     Minimax error
!> \param err_ctff   Cutoff error
!> \param C_mm       Scaling constant to generalize AM-GM upper bound estimate to
!>                   minimax approx.
!> \param para_env ...
! **************************************************************************************************
   SUBROUTINE cutoff_minimax_error(cutoff, hmat, h_inv, vol, G_min, zet_min, l_mm, zet_max, l_max_zet, &
                                   n_minimax, minimax_aw, err_mm, err_ctff, C_mm, para_env)
      REAL(kind=dp), INTENT(IN)                          :: cutoff
      REAL(kind=dp), DIMENSION(3, 3), INTENT(IN)         :: hmat, h_inv
      REAL(kind=dp), INTENT(IN)                          :: vol, g_min, zet_min
      INTEGER, INTENT(IN)                                :: l_mm
      REAL(kind=dp), INTENT(IN)                          :: zet_max
      INTEGER, INTENT(IN)                                :: l_max_zet, n_minimax
      REAL(kind=dp), DIMENSION(:), INTENT(OUT)           :: minimax_aw
      REAL(kind=dp), INTENT(OUT)                         :: err_mm, err_ctff, c_mm
      TYPE(mp_para_env_type), INTENT(IN)                 :: para_env
 
      REAL(kind=dp)                                      :: delta_mm
 
      CALL minimax_error(cutoff, hmat, vol, g_min, zet_min, l_mm, &
                         n_minimax, minimax_aw, err_mm, delta_mm)
      CALL cutoff_error(cutoff, h_inv, g_min, zet_max, l_max_zet, &
                        n_minimax, minimax_aw, err_ctff, c_mm, para_env)
 
   END SUBROUTINE cutoff_minimax_error
 
! **************************************************************************************************
!> \brief   Minimax error, simple analytical formula
!>          Note minimax error may blow up for small exponents. This is also observed numerically,
!>          but in this case, error estimate is no upper bound.
!> \param cutoff ...
!> \param hmat ...
!> \param vol ...
!> \param G_min ...
!> \param zet_min    Exponent of P to estimate minimax error
!> \param l_mm       total ang. mom. quantum number of P to estimate minimax error
!> \param n_minimax  Number of terms in minimax approximation
!> \param minimax_aw Minimax coefficients
!> \param err_mm     Minimax error
!> \param delta_mm ...
!> \param potential ...
!> \param pot_par ...
! **************************************************************************************************
   SUBROUTINE minimax_error(cutoff, hmat, vol, G_min, zet_min, l_mm, &
                            n_minimax, minimax_aw, err_mm, delta_mm, potential, pot_par)
      REAL(kind=dp), INTENT(IN)                          :: cutoff
      REAL(kind=dp), DIMENSION(3, 3), INTENT(IN)         :: hmat
      REAL(kind=dp), INTENT(IN)                          :: vol, g_min, zet_min
      INTEGER, INTENT(IN)                                :: l_mm, n_minimax
      REAL(kind=dp), DIMENSION(:), INTENT(OUT)           :: minimax_aw
      REAL(kind=dp), INTENT(OUT)                         :: err_mm, delta_mm
      INTEGER, INTENT(IN), OPTIONAL                      :: potential
      REAL(kind=dp), INTENT(IN), OPTIONAL                :: pot_par
 
      INTEGER                                            :: i_xyz
      REAL(kind=dp)                                      :: prod_mm_k
 
      CALL get_minimax_coeff_v_gspace(n_minimax, cutoff, g_min, minimax_aw(:), &
                                      potential=potential, pot_par=pot_par, err_minimax=delta_mm)
 
      prod_mm_k = 1.0_dp
      DO i_xyz = 1, 3
         prod_mm_k = prod_mm_k*(abs(hmat(i_xyz, i_xyz))/twopi + &
                                merge(sqrt(2.0_dp/(zet_min*pi))*exp(-1.0_dp), 0.0_dp, l_mm .GT. 0))
      END DO
      err_mm = 32*pi**4/vol*delta_mm*prod_mm_k
 
   END SUBROUTINE
 
! **************************************************************************************************
!> \brief Cutoff error, estimating G > G_c part of Ewald sum by using C/3 * 1/(Gx^2*Gy^2*Gz^2)^1/3 as an
!>        upper bound for 1/G^2 (AM-GM inequality) and its minimax approximation (factor C).
!>
!> Note: usually, minimax approx. falls off faster than 1/G**2, so C should be approximately 1.
!> The error is calculated for all l up to l_max and golden section search algorithm is
!> applied to find the exponent that maximizes cutoff error.
!> \param cutoff ...
!> \param h_inv ...
!> \param G_min ...
!> \param zet_max   Max. exponents of P to estimate cutoff error
!> \param l_max_zet     Max. total ang. mom. quantum numbers of P to estimate cutoff error
!> \param n_minimax  Number of terms in minimax approximation
!> \param minimax_aw Minimax coefficients
!> \param err_ctff   Cutoff error
!> \param C_mm       Scaling constant to generalize AM-GM upper bound estimate to
!>                   minimax approx.
!> \param para_env ...
! **************************************************************************************************
   SUBROUTINE cutoff_error(cutoff, h_inv, G_min, zet_max, l_max_zet, &
                           n_minimax, minimax_aw, err_ctff, C_mm, para_env)
      REAL(kind=dp), INTENT(IN)                          :: cutoff
      REAL(kind=dp), DIMENSION(3, 3), INTENT(IN)         :: h_inv
      REAL(kind=dp), INTENT(IN)                          :: g_min, zet_max
      INTEGER, INTENT(IN)                                :: l_max_zet, n_minimax
      REAL(kind=dp), DIMENSION(:), INTENT(INOUT)         :: minimax_aw
      REAL(kind=dp), INTENT(OUT)                         :: err_ctff, c_mm
      TYPE(mp_para_env_type), INTENT(IN)                 :: para_env
 
      INTEGER                                            :: i_aw, ig, iter, max_iter, ng
      REAL(kind=dp) :: c, dg, eps_zet, err0, err1, err_c, err_ctff_curr, err_ctff_prev, err_d, g, &
         g_1, g_c, gr, zet_a, zet_b, zet_c, zet_d, zet_div, zet_max_tmp
 
      ! parameters for finding exponent maximizing cutoff error
 
      eps_zet = 1.0e-05_dp ! tolerance for exponent
      zet_div = 2.0_dp ! sampling constant for finding initial values of exponents
      max_iter = 100 ! maximum number of iterations in golden section search
      g_c = sqrt(2.0*cutoff)
 
      zet_max_tmp = zet_max
 
      ! 2) Cutoff error, estimating G > G_c part of Ewald sum by using C/3 * 1/(Gx^2*Gy^2*Gz^2)^1/3 as an
      !                  upper bound for 1/G^2 (AM-GM inequality) and its minimax approximation (factor C).
      !                  Note: usually, minimax approx. falls off faster than 1/G**2, so C should be approximately 1.
      !                  The error is calculated for all l up to l_max and golden section search algorithm is
      !                  applied to find the exponent that maximizes cutoff error.
      g_1 = sqrt(1.0_dp/(3.0_dp*minval(minimax_aw(1:n_minimax))))
 
      c_mm = 0.0_dp
      IF (g_1 .GT. g_c) THEN
         ng = 1000
         dg = (g_1 - g_c)/ng
         g = g_c
         DO ig = 1, ng
            g = min(g, g_c)
            c = 0.0_dp
            DO i_aw = 1, n_minimax
               c = c + 3.0_dp*minimax_aw(n_minimax + i_aw)*exp(-3.0_dp*minimax_aw(i_aw)*g**2)*g**2
            END DO
            c_mm = max(c, c_mm)
            g = g + dg
         END DO
      ELSE
         DO i_aw = 1, n_minimax
            c_mm = c_mm + 3.0_dp*minimax_aw(n_minimax + i_aw)*exp(-3.0_dp*minimax_aw(i_aw)*g_c**2)*g_c**2
         END DO
      END IF
      c = max(1.0_dp, c_mm)
 
      err_ctff_prev = 0.0_dp
      gr = 0.5_dp*(sqrt(5.0_dp) - 1.0_dp) ! golden ratio
      ! Find valid starting values for golden section search
      DO iter = 1, max_iter + 1
         IF (iter .GT. max_iter) &
            CALL cp_abort(__location__, "Maximum number of iterations for finding "// &
                          "exponent maximizing cutoff error has been exceeded.")
 
         CALL cutoff_error_fixed_exp(cutoff, h_inv, g_min, l_max_zet, zet_max_tmp, c, err_ctff_curr, para_env)
         IF (err_ctff_prev .GE. err_ctff_curr) THEN
            zet_a = zet_max_tmp
            zet_b = min(zet_max_tmp*zet_div**2, zet_max)
            EXIT
         ELSE
            err_ctff_prev = err_ctff_curr
         END IF
         zet_max_tmp = zet_max_tmp/zet_div
      END DO
 
      ! Golden section search
      zet_c = zet_b - gr*(zet_b - zet_a)
      zet_d = zet_a + gr*(zet_b - zet_a)
      DO iter = 1, max_iter + 1
         IF (abs(zet_c - zet_d) .LT. eps_zet*(zet_a + zet_b)) THEN
            CALL cutoff_error_fixed_exp(cutoff, h_inv, g_min, l_max_zet, zet_a, c, err0, para_env)
            CALL cutoff_error_fixed_exp(cutoff, h_inv, g_min, l_max_zet, zet_b, c, err1, para_env)
            err_ctff_curr = max(err0, err1)
            EXIT
         END IF
         CALL cutoff_error_fixed_exp(cutoff, h_inv, g_min, l_max_zet, zet_c, c, err_c, para_env)
         CALL cutoff_error_fixed_exp(cutoff, h_inv, g_min, l_max_zet, zet_d, c, err_d, para_env)
         IF (err_c .GT. err_d) THEN
            zet_b = zet_d
            zet_d = zet_c
            zet_c = zet_b - gr*(zet_b - zet_a)
         ELSE
            zet_a = zet_c
            zet_c = zet_d
            zet_d = zet_a + gr*(zet_b - zet_a)
         END IF
      END DO
      err_ctff = err_ctff_curr
 
   END SUBROUTINE
 
! **************************************************************************************************
!> \brief Calculate cutoff error estimate by using C_mm/3 * 1/(Gx^2*Gy^2*Gz^2)^1/3
!>        as an upper bound for 1/G^2 (and its minimax approximation) for |G| > G_c.
!>        Error is referring to a basis function P with fixed exponent zet_max and
!>        max. angular momentum l_max_zet.
!> \param cutoff ...
!> \param h_inv ...
!> \param G_min ...
!> \param l_max_zet ...
!> \param zet_max ...
!> \param C_mm ...
!> \param err_c ...
!> \param para_env ...
! **************************************************************************************************
   SUBROUTINE cutoff_error_fixed_exp(cutoff, h_inv, G_min, l_max_zet, zet_max, C_mm, err_c, para_env)
      REAL(kind=dp), INTENT(IN)                          :: cutoff
      REAL(kind=dp), DIMENSION(3, 3), INTENT(IN)         :: h_inv
      REAL(kind=dp), INTENT(IN)                          :: g_min
      INTEGER, INTENT(IN)                                :: l_max_zet
      REAL(kind=dp), INTENT(IN)                          :: zet_max, c_mm
      REAL(kind=dp), INTENT(OUT)                         :: err_c
      TYPE(mp_para_env_type), INTENT(IN)                 :: para_env
 
      INTEGER                                            :: ax, ay, az, g_l, g_u, gl_first, gl_last, &
                                                            gu_first, gu_last, i_xyz, l, my_p, &
                                                            n_gl, n_gl_left, n_gl_p, n_gu, &
                                                            n_gu_left, n_gu_p, n_p
      REAL(kind=dp)                                      :: alpha_g, eps_g, err_c_l, g_c, g_rad, &
                                                            g_res, inv_lgth, prefactor, sum_g_diff
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :)        :: s_g_l, s_g_u
 
      g_c = sqrt(2.0_dp*cutoff)
      eps_g = tiny(eps_g) ! sum up to machine precision
      g_res = 0.5_dp*g_min ! resolution for screening
 
      err_c = 0.0_dp
      alpha_g = 1.0_dp/(2.0_dp*zet_max)
      prefactor = 1.0_dp/zet_max
 
      ALLOCATE (s_g_l(0:2*l_max_zet, 3))
      ALLOCATE (s_g_u(0:2*l_max_zet, 3))
 
      g_rad = exp_radius(2*l_max_zet, alpha_g, eps_g, prefactor, epsabs=g_res)
 
      ! Parallelization of sum over G vectors
      my_p = para_env%mepos ! mpi rank
      n_p = para_env%num_pe ! total number of processes
 
      DO i_xyz = 1, 3
         inv_lgth = abs(h_inv(i_xyz, i_xyz))
 
         g_l = floor(g_c/(inv_lgth*twopi))
         g_u = floor(g_rad/(inv_lgth*twopi))
 
         IF (g_u .LT. g_l) g_u = g_l
 
         ! Serial code:
         ! !Sum |G| <= G_c
         ! CALL pgf_sum_2c_gspace_1d_deltal(S_G_l(:, i_xyz), alpha_G, inv_lgth, -G_l, G_l, &
         !                               2.0_dp/3.0_dp, prefactor)
         ! !Sum |G| > G_c
         ! CALL pgf_sum_2c_gspace_1d_deltal(S_G_u(:, i_xyz), alpha_G, inv_lgth, G_l + 1, G_u, &
         !                               2.0_dp/3.0_dp, prefactor)
 
         ! Parallel code:
         n_gu = max((g_u - g_l), 0)
         n_gl = 2*g_l + 1
         n_gu_p = n_gu/n_p
         n_gl_p = n_gl/n_p
         n_gu_left = mod(n_gu, n_p)
         n_gl_left = mod(n_gl, n_p)
 
         IF (my_p .LT. n_gu_left) THEN
            gu_first = g_l + 1 + (n_gu_p + 1)*my_p
            gu_last = g_l + 1 + (n_gu_p + 1)*(my_p + 1) - 1
         ELSE
            gu_first = g_l + 1 + n_gu_left + n_gu_p*my_p
            gu_last = g_l + 1 + n_gu_left + n_gu_p*(my_p + 1) - 1
         END IF
 
         IF (my_p .LT. n_gl_left) THEN
            gl_first = -g_l + (n_gl_p + 1)*my_p
            gl_last = -g_l + (n_gl_p + 1)*(my_p + 1) - 1
         ELSE
            gl_first = -g_l + n_gl_left + n_gl_p*my_p
            gl_last = -g_l + n_gl_left + n_gl_p*(my_p + 1) - 1
         END IF
 
         ! Sum |G| <= G_c
         CALL pgf_sum_2c_gspace_1d_deltal(s_g_l(:, i_xyz), alpha_g, inv_lgth, gl_first, gl_last, &
                                          2.0_dp/3.0_dp, prefactor)
 
         ! Sum |G| > G_c
         CALL pgf_sum_2c_gspace_1d_deltal(s_g_u(:, i_xyz), alpha_g, inv_lgth, gu_first, gu_last, &
                                          2.0_dp/3.0_dp, prefactor)
      END DO
 
      CALL para_env%sum(s_g_l)
      CALL para_env%sum(s_g_u)
 
      s_g_u = s_g_u*2.0_dp ! to include negative values of G
 
      DO l = 0, l_max_zet
      DO ax = 0, l
      DO ay = 0, l - ax
         az = l - ax - ay
 
         ! Compute prod_k (S_G_l(l_k,k) + S_G_u(l_k,k)) - prod_k (S_G_l(l_k,k)) with k in {x, y, z}
         ! Note: term by term multiplication to avoid subtraction for numerical stability
         sum_g_diff = s_g_u(2*ax, 1)*s_g_u(2*ay, 2)*s_g_u(2*az, 3) + &
                      s_g_u(2*ax, 1)*s_g_u(2*ay, 2)*s_g_l(2*az, 3) + &
                      s_g_u(2*ax, 1)*s_g_l(2*ay, 2)*s_g_u(2*az, 3) + &
                      s_g_l(2*ax, 1)*s_g_u(2*ay, 2)*s_g_u(2*az, 3) + &
                      s_g_u(2*ax, 1)*s_g_l(2*ay, 2)*s_g_l(2*az, 3) + &
                      s_g_l(2*ax, 1)*s_g_u(2*ay, 2)*s_g_l(2*az, 3) + &
                      s_g_l(2*ax, 1)*s_g_l(2*ay, 2)*s_g_u(2*az, 3)
 
         err_c_l = 4.0_dp*pi**4*hermite_gauss_norm(zet_max, [ax, ay, az])**2* &
                   c_mm/3.0_dp*sum_g_diff
 
         err_c = max(err_c, err_c_l)
      END DO
      END DO
      END DO
 
      DEALLOCATE (s_g_u, s_g_l)
 
   END SUBROUTINE cutoff_error_fixed_exp
 
END MODULE eri_mme_error_control