LCOV - code coverage report
Current view: top level - src/eri_mme - eri_mme_error_control.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 87.8 % 172 151
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 5 5

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \brief Methods aiming for error estimate and automatic cutoff calibration.
      10              : !>        integrals.
      11              : !> \par History
      12              : !>       2015 09 created
      13              : !> \author Patrick Seewald
      14              : ! **************************************************************************************************
      15              : 
      16              : MODULE eri_mme_error_control
      17              :    USE ao_util,                         ONLY: exp_radius
      18              :    USE eri_mme_gaussian,                ONLY: get_minimax_coeff_v_gspace,&
      19              :                                               hermite_gauss_norm
      20              :    USE eri_mme_lattice_summation,       ONLY: pgf_sum_2c_gspace_1d_deltal
      21              :    USE kinds,                           ONLY: dp
      22              :    USE mathconstants,                   ONLY: pi,&
      23              :                                               twopi
      24              :    USE message_passing,                 ONLY: mp_para_env_type
      25              : #include "../base/base_uses.f90"
      26              : 
      27              :    IMPLICIT NONE
      28              : 
      29              :    PRIVATE
      30              : 
      31              :    LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .FALSE.
      32              : 
      33              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'eri_mme_error_control'
      34              : 
      35              :    PUBLIC :: calibrate_cutoff, cutoff_minimax_error, minimax_error, cutoff_error
      36              : CONTAINS
      37              : 
      38              : ! **************************************************************************************************
      39              : !> \brief Find optimal cutoff minimizing errors due to minimax approximation and
      40              : !>        due to finite cutoff using bisection on the difference of the errors
      41              : !> \param hmat ...
      42              : !> \param h_inv ...
      43              : !> \param G_min ...
      44              : !> \param vol ...
      45              : !> \param zet_min   Minimum exponent
      46              : !> \param l_mm      Total ang. mom. quantum number
      47              : !> \param zet_max     Max. exponents to estimate cutoff error
      48              : !> \param l_max_zet       Max. total ang. mom. quantum numbers to estimate cutoff error
      49              : !> \param n_minimax Number of terms in minimax approximation
      50              : !> \param cutoff_l  Initial guess of lower bound for cutoff
      51              : !> \param cutoff_r  Initial guess of upper bound for cutoff
      52              : !> \param tol       Tolerance (cutoff precision)
      53              : !> \param delta     to modify initial guess interval
      54              : !> \param cutoff    Best cutoff
      55              : !> \param err_mm    Minimax error
      56              : !> \param err_c     Cutoff error
      57              : !> \param C_mm      Scaling constant to generalize AM-GM upper bound estimate to
      58              : !>                  minimax approx.
      59              : !> \param para_env ...
      60              : !> \param print_calib ...
      61              : !> \param unit_nr ...
      62              : ! **************************************************************************************************
      63           84 :    SUBROUTINE calibrate_cutoff(hmat, h_inv, G_min, vol, zet_min, l_mm, zet_max, l_max_zet, &
      64              :                                n_minimax, cutoff_l, cutoff_r, tol, delta, &
      65              :                                cutoff, err_mm, err_c, C_mm, para_env, print_calib, unit_nr)
      66              :       REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN)         :: hmat, h_inv
      67              :       REAL(KIND=dp), INTENT(IN)                          :: G_min
      68              :       REAL(KIND=dp)                                      :: vol
      69              :       REAL(KIND=dp), INTENT(IN)                          :: zet_min
      70              :       INTEGER, INTENT(IN)                                :: l_mm
      71              :       REAL(KIND=dp), INTENT(IN)                          :: zet_max
      72              :       INTEGER, INTENT(IN)                                :: l_max_zet, n_minimax
      73              :       REAL(KIND=dp), INTENT(IN)                          :: cutoff_l, cutoff_r, tol, delta
      74              :       REAL(KIND=dp), INTENT(OUT)                         :: cutoff, err_mm, err_c, C_mm
      75              :       TYPE(mp_para_env_type), INTENT(IN)                 :: para_env
      76              :       LOGICAL, INTENT(IN)                                :: print_calib
      77              :       INTEGER, INTENT(IN)                                :: unit_nr
      78              : 
      79              :       INTEGER                                            :: i, iter1, iter2, max_iter
      80              :       LOGICAL                                            :: do_print, valid_initial
      81              :       REAL(KIND=dp)                                      :: cutoff_mid, delta_c_mid, delta_mm_mid
      82           84 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: minimax_aw
      83              :       REAL(KIND=dp), DIMENSION(2)                        :: cutoff_lr, delta_c, delta_mm
      84              : 
      85           84 :       do_print = unit_nr > 0 .AND. print_calib
      86              :       IF (do_print) THEN
      87            0 :          WRITE (unit_nr, '(/T2, A)') "ERI_MME| Basis set parameters for estimating minimax error"
      88            0 :          WRITE (unit_nr, '(T2, A, T67, ES12.2, 1X, I1)') "ERI_MME|   exp, l:", zet_min, l_mm
      89            0 :          WRITE (unit_nr, '(T2, A)') "ERI_MME| Basis set parameters for estimating cutoff error"
      90            0 :          WRITE (unit_nr, '(T2, A, T67, ES12.2, 1X, I1)') "ERI_MME|   exp, l:", zet_max, l_max_zet
      91              :       END IF
      92              : 
      93           84 :       max_iter = 100
      94              : 
      95           84 :       IF ((cutoff_r - cutoff_l)/(0.5_dp*(cutoff_r + cutoff_l)) .LE. tol) &
      96              :          CALL cp_abort(__LOCATION__, "difference of boundaries for cutoff "// &
      97            0 :                        "(MAX - MIN) must be greater than cutoff precision.")
      98              : 
      99           84 :       IF ((delta .GE. 1.0_dp) .OR. (delta .LE. 0.0_dp)) &
     100              :          CALL cp_abort(__LOCATION__, &
     101            0 :                        "relative delta to modify initial cutoff interval (DELTA) must be in (0, 1)")
     102              : 
     103           84 :       cutoff_lr(1) = cutoff_l
     104           84 :       cutoff_lr(2) = cutoff_r
     105              : 
     106          252 :       ALLOCATE (minimax_aw(2*n_minimax))
     107              : 
     108           84 :       IF (do_print) THEN
     109            0 :          WRITE (unit_nr, '(/T2, A)') "ERI_MME| Calibrating cutoff by bisecting error(minimax) - error(cutoff)"
     110            0 :          WRITE (unit_nr, '(T2, A, T72, ES9.2)') "ERI_MME| Rel. cutoff precision", tol
     111            0 :          WRITE (unit_nr, '(T2, A, T77, F4.1)') "ERI_MME| Rel. cutoff delta to modify initial interval", delta
     112              :       END IF
     113              : 
     114              :       ! 1) find valid initial values for bisection
     115           84 :       DO iter1 = 1, max_iter + 1
     116           84 :          IF (iter1 .GT. max_iter) &
     117              :             CALL cp_abort(__LOCATION__, &
     118              :                           "Maximum number of iterations in bisection to determine initial "// &
     119            0 :                           "cutoff interval has been exceeded.")
     120              : 
     121           84 :          cutoff_lr(1) = MAX(cutoff_lr(1), 0.5_dp*G_min**2)
     122              :          ! approx.) is hit
     123              : 
     124          252 :          DO i = 1, 2
     125              :             CALL cutoff_minimax_error(cutoff_lr(i), hmat, h_inv, vol, G_min, zet_min, l_mm, zet_max, l_max_zet, &
     126          252 :                                       n_minimax, minimax_aw, delta_mm(i), delta_c(i), C_mm, para_env)
     127              :          END DO
     128              : 
     129           84 :          valid_initial = .TRUE.
     130           84 :          IF ((delta_mm(1) - delta_c(1)) .GT. 0) THEN
     131            0 :             cutoff_lr(1) = cutoff_lr(1)*(1.0_dp - ABS(delta))
     132            0 :             valid_initial = .FALSE.
     133              :          END IF
     134           84 :          IF ((delta_mm(2) - delta_c(2)) .LT. 0) THEN
     135            0 :             cutoff_lr(2) = cutoff_lr(2)*(1.0_dp + ABS(delta))
     136              :             valid_initial = .FALSE.
     137              :          END IF
     138              : 
     139           84 :          IF (valid_initial) EXIT
     140              :       END DO
     141              : 
     142              :       ! 2) bisection to find cutoff s.t. err_minimax(cutoff) - err_cutoff(cutoff) = 0
     143           84 :       IF (do_print) WRITE (unit_nr, '(/T2, A)') &
     144            0 :          "ERI_MME| Step, cutoff (min, max, mid), err(minimax), err(cutoff), err diff"
     145              : 
     146         1190 :       DO iter2 = 1, max_iter + 1
     147         1190 :          IF (iter2 .GT. max_iter) &
     148              :             CALL cp_abort(__LOCATION__, &
     149            0 :                           "Maximum number of iterations in bisection to determine cutoff has been exceeded")
     150              : 
     151         1190 :          cutoff_mid = 0.5_dp*(cutoff_lr(1) + cutoff_lr(2))
     152              :          CALL cutoff_minimax_error(cutoff_mid, hmat, h_inv, vol, G_min, zet_min, l_mm, zet_max, l_max_zet, &
     153         1190 :                                    n_minimax, minimax_aw, delta_mm_mid, delta_c_mid, C_mm, para_env)
     154         1190 :          IF (do_print) WRITE (unit_nr, '(T11, I2, F11.1, F11.1, F11.1, 3X, ES9.2, 3X, ES9.2, 3X, ES9.2)') &
     155            0 :             iter2, cutoff_lr(1), cutoff_lr(2), cutoff_mid, &
     156            0 :             delta_mm_mid, delta_c_mid, delta_mm_mid - delta_c_mid
     157              : 
     158         1190 :          IF ((cutoff_lr(2) - cutoff_lr(1))/cutoff_mid .LT. tol) EXIT
     159         2380 :          IF (delta_mm_mid - delta_c_mid .GT. 0) THEN
     160          776 :             cutoff_lr(2) = cutoff_mid
     161          776 :             delta_mm(2) = delta_mm_mid
     162          776 :             delta_c(2) = delta_c_mid
     163              :          ELSE
     164          330 :             cutoff_lr(1) = cutoff_mid
     165          330 :             delta_mm(1) = delta_mm_mid
     166          330 :             delta_c(1) = delta_c_mid
     167              :          END IF
     168              :       END DO
     169           84 :       err_mm = delta_mm_mid
     170           84 :       err_c = delta_c_mid
     171           84 :       cutoff = cutoff_mid
     172              : 
     173           84 :       IF (do_print) THEN
     174            0 :          WRITE (unit_nr, '(/T2, A)') "ERI_MME| Cutoff calibration number of steps:"
     175            0 :          WRITE (unit_nr, '(T2, A, T79, I2)') "ERI_MME|   Steps for initial interval", iter1 - 1
     176            0 :          WRITE (unit_nr, '(T2, A, T79, I2/)') "ERI_MME|   Bisection iteration steps", iter2 - 1
     177              :       END IF
     178              : 
     179           84 :    END SUBROUTINE calibrate_cutoff
     180              : 
     181              : ! **************************************************************************************************
     182              : !> \brief Compute upper bounds for the errors of 2-center ERI's (P|P) due
     183              : !>        to minimax approximation and due to finite cutoff, where P is a
     184              : !>        normalized Hermite Gaussian.
     185              : !> \param cutoff ...
     186              : !> \param hmat ...
     187              : !> \param h_inv ...
     188              : !> \param vol ...
     189              : !> \param G_min ...
     190              : !> \param zet_min     Exponent of P to estimate minimax error
     191              : !> \param l_mm       total ang. mom. quantum number of P to estimate minimax error
     192              : !> \param zet_max   Max. exponents of P to estimate cutoff error
     193              : !> \param l_max_zet     Max. total ang. mom. quantum numbers of P to estimate cutoff error
     194              : !> \param n_minimax  Number of terms in minimax approximation
     195              : !> \param minimax_aw Minimax coefficients
     196              : !> \param err_mm     Minimax error
     197              : !> \param err_ctff   Cutoff error
     198              : !> \param C_mm       Scaling constant to generalize AM-GM upper bound estimate to
     199              : !>                   minimax approx.
     200              : !> \param para_env ...
     201              : ! **************************************************************************************************
     202         1396 :    SUBROUTINE cutoff_minimax_error(cutoff, hmat, h_inv, vol, G_min, zet_min, l_mm, zet_max, l_max_zet, &
     203         1396 :                                    n_minimax, minimax_aw, err_mm, err_ctff, C_mm, para_env)
     204              :       REAL(KIND=dp), INTENT(IN)                          :: cutoff
     205              :       REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN)         :: hmat, h_inv
     206              :       REAL(KIND=dp), INTENT(IN)                          :: vol, G_min, zet_min
     207              :       INTEGER, INTENT(IN)                                :: l_mm
     208              :       REAL(KIND=dp), INTENT(IN)                          :: zet_max
     209              :       INTEGER, INTENT(IN)                                :: l_max_zet, n_minimax
     210              :       REAL(KIND=dp), DIMENSION(:), INTENT(OUT)           :: minimax_aw
     211              :       REAL(KIND=dp), INTENT(OUT)                         :: err_mm, err_ctff, C_mm
     212              :       TYPE(mp_para_env_type), INTENT(IN)                 :: para_env
     213              : 
     214              :       REAL(KIND=dp)                                      :: delta_mm
     215              : 
     216              :       CALL minimax_error(cutoff, hmat, vol, G_min, zet_min, l_mm, &
     217         1396 :                          n_minimax, minimax_aw, err_mm, delta_mm)
     218              :       CALL cutoff_error(cutoff, h_inv, G_min, zet_max, l_max_zet, &
     219         1396 :                         n_minimax, minimax_aw, err_ctff, C_mm, para_env)
     220              : 
     221         1396 :    END SUBROUTINE cutoff_minimax_error
     222              : 
     223              : ! **************************************************************************************************
     224              : !> \brief   Minimax error, simple analytical formula
     225              : !>          Note minimax error may blow up for small exponents. This is also observed numerically,
     226              : !>          but in this case, error estimate is no upper bound.
     227              : !> \param cutoff ...
     228              : !> \param hmat ...
     229              : !> \param vol ...
     230              : !> \param G_min ...
     231              : !> \param zet_min    Exponent of P to estimate minimax error
     232              : !> \param l_mm       total ang. mom. quantum number of P to estimate minimax error
     233              : !> \param n_minimax  Number of terms in minimax approximation
     234              : !> \param minimax_aw Minimax coefficients
     235              : !> \param err_mm     Minimax error
     236              : !> \param delta_mm ...
     237              : !> \param potential ...
     238              : !> \param pot_par ...
     239              : ! **************************************************************************************************
     240        86888 :    SUBROUTINE minimax_error(cutoff, hmat, vol, G_min, zet_min, l_mm, &
     241        86888 :                             n_minimax, minimax_aw, err_mm, delta_mm, potential, pot_par)
     242              :       REAL(KIND=dp), INTENT(IN)                          :: cutoff
     243              :       REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN)         :: hmat
     244              :       REAL(KIND=dp), INTENT(IN)                          :: vol, G_min, zet_min
     245              :       INTEGER, INTENT(IN)                                :: l_mm, n_minimax
     246              :       REAL(KIND=dp), DIMENSION(:), INTENT(OUT)           :: minimax_aw
     247              :       REAL(KIND=dp), INTENT(OUT)                         :: err_mm, delta_mm
     248              :       INTEGER, INTENT(IN), OPTIONAL                      :: potential
     249              :       REAL(KIND=dp), INTENT(IN), OPTIONAL                :: pot_par
     250              : 
     251              :       INTEGER                                            :: i_xyz
     252              :       REAL(KIND=dp)                                      :: prod_mm_k
     253              : 
     254              :       CALL get_minimax_coeff_v_gspace(n_minimax, cutoff, G_min, minimax_aw(:), &
     255        86888 :                                       potential=potential, pot_par=pot_par, err_minimax=delta_mm)
     256              : 
     257        86888 :       prod_mm_k = 1.0_dp
     258       347552 :       DO i_xyz = 1, 3
     259              :          prod_mm_k = prod_mm_k*(ABS(hmat(i_xyz, i_xyz))/twopi + &
     260       348032 :                                 MERGE(SQRT(2.0_dp/(zet_min*pi))*EXP(-1.0_dp), 0.0_dp, l_mm .GT. 0))
     261              :       END DO
     262        86888 :       err_mm = 32*pi**4/vol*delta_mm*prod_mm_k
     263              : 
     264        86888 :    END SUBROUTINE
     265              : 
     266              : ! **************************************************************************************************
     267              : !> \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
     268              : !>        upper bound for 1/G^2 (AM-GM inequality) and its minimax approximation (factor C).
     269              : !>
     270              : !> Note: usually, minimax approx. falls off faster than 1/G**2, so C should be approximately 1.
     271              : !> The error is calculated for all l up to l_max and golden section search algorithm is
     272              : !> applied to find the exponent that maximizes cutoff error.
     273              : !> \param cutoff ...
     274              : !> \param h_inv ...
     275              : !> \param G_min ...
     276              : !> \param zet_max   Max. exponents of P to estimate cutoff error
     277              : !> \param l_max_zet     Max. total ang. mom. quantum numbers of P to estimate cutoff error
     278              : !> \param n_minimax  Number of terms in minimax approximation
     279              : !> \param minimax_aw Minimax coefficients
     280              : !> \param err_ctff   Cutoff error
     281              : !> \param C_mm       Scaling constant to generalize AM-GM upper bound estimate to
     282              : !>                   minimax approx.
     283              : !> \param para_env ...
     284              : ! **************************************************************************************************
     285         1396 :    SUBROUTINE cutoff_error(cutoff, h_inv, G_min, zet_max, l_max_zet, &
     286         1396 :                            n_minimax, minimax_aw, err_ctff, C_mm, para_env)
     287              :       REAL(KIND=dp), INTENT(IN)                          :: cutoff
     288              :       REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN)         :: h_inv
     289              :       REAL(KIND=dp), INTENT(IN)                          :: G_min, zet_max
     290              :       INTEGER, INTENT(IN)                                :: l_max_zet, n_minimax
     291              :       REAL(KIND=dp), DIMENSION(:), INTENT(INOUT)         :: minimax_aw
     292              :       REAL(KIND=dp), INTENT(OUT)                         :: err_ctff, C_mm
     293              :       TYPE(mp_para_env_type), INTENT(IN)                 :: para_env
     294              : 
     295              :       INTEGER                                            :: i_aw, iG, iter, max_iter, nG
     296              :       REAL(KIND=dp) :: C, dG, eps_zet, err0, err1, err_c, err_ctff_curr, err_ctff_prev, err_d, G, &
     297              :          G_1, G_c, gr, zet_a, zet_b, zet_c, zet_d, zet_div, zet_max_tmp
     298              : 
     299              :       ! parameters for finding exponent maximizing cutoff error
     300              : 
     301         1396 :       eps_zet = 1.0E-05_dp ! tolerance for exponent
     302         1396 :       zet_div = 2.0_dp ! sampling constant for finding initial values of exponents
     303         1396 :       max_iter = 100 ! maximum number of iterations in golden section search
     304         1396 :       G_c = SQRT(2.0*cutoff)
     305              : 
     306         1396 :       zet_max_tmp = zet_max
     307              : 
     308              :       ! 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
     309              :       !                  upper bound for 1/G^2 (AM-GM inequality) and its minimax approximation (factor C).
     310              :       !                  Note: usually, minimax approx. falls off faster than 1/G**2, so C should be approximately 1.
     311              :       !                  The error is calculated for all l up to l_max and golden section search algorithm is
     312              :       !                  applied to find the exponent that maximizes cutoff error.
     313        25586 :       G_1 = SQRT(1.0_dp/(3.0_dp*MINVAL(minimax_aw(1:n_minimax))))
     314              : 
     315         1396 :       C_mm = 0.0_dp
     316         1396 :       IF (G_1 .GT. G_c) THEN
     317          762 :          nG = 1000
     318          762 :          dG = (G_1 - G_c)/nG
     319          762 :          G = G_c
     320       762762 :          DO iG = 1, nG
     321       762000 :             G = MIN(G, G_c)
     322       762000 :             C = 0.0_dp
     323     20478000 :             DO i_aw = 1, n_minimax
     324     20478000 :                C = C + 3.0_dp*minimax_aw(n_minimax + i_aw)*EXP(-3.0_dp*minimax_aw(i_aw)*G**2)*G**2
     325              :             END DO
     326       762000 :             C_mm = MAX(C, C_mm)
     327       762762 :             G = G + dG
     328              :          END DO
     329              :       ELSE
     330         3712 :          DO i_aw = 1, n_minimax
     331         3712 :             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
     332              :          END DO
     333              :       END IF
     334         1396 :       C = MAX(1.0_dp, C_mm)
     335              : 
     336         1396 :       err_ctff_prev = 0.0_dp
     337         1396 :       gr = 0.5_dp*(SQRT(5.0_dp) - 1.0_dp) ! golden ratio
     338              :       ! Find valid starting values for golden section search
     339         2754 :       DO iter = 1, max_iter + 1
     340         2754 :          IF (iter .GT. max_iter) &
     341              :             CALL cp_abort(__LOCATION__, "Maximum number of iterations for finding "// &
     342            0 :                           "exponent maximizing cutoff error has been exceeded.")
     343              : 
     344         2754 :          CALL cutoff_error_fixed_exp(cutoff, h_inv, G_min, l_max_zet, zet_max_tmp, C, err_ctff_curr, para_env)
     345         2754 :          IF (err_ctff_prev .GE. err_ctff_curr) THEN
     346         1396 :             zet_a = zet_max_tmp
     347         1396 :             zet_b = MIN(zet_max_tmp*zet_div**2, zet_max)
     348         1396 :             EXIT
     349              :          ELSE
     350         1358 :             err_ctff_prev = err_ctff_curr
     351              :          END IF
     352         4112 :          zet_max_tmp = zet_max_tmp/zet_div
     353              :       END DO
     354              : 
     355              :       ! Golden section search
     356         1396 :       zet_c = zet_b - gr*(zet_b - zet_a)
     357         1396 :       zet_d = zet_a + gr*(zet_b - zet_a)
     358        23210 :       DO iter = 1, max_iter + 1
     359        23210 :          IF (ABS(zet_c - zet_d) .LT. eps_zet*(zet_a + zet_b)) THEN
     360         1396 :             CALL cutoff_error_fixed_exp(cutoff, h_inv, G_min, l_max_zet, zet_a, C, err0, para_env)
     361         1396 :             CALL cutoff_error_fixed_exp(cutoff, h_inv, G_min, l_max_zet, zet_b, C, err1, para_env)
     362         1396 :             err_ctff_curr = MAX(err0, err1)
     363         1396 :             EXIT
     364              :          END IF
     365        21814 :          CALL cutoff_error_fixed_exp(cutoff, h_inv, G_min, l_max_zet, zet_c, C, err_c, para_env)
     366        21814 :          CALL cutoff_error_fixed_exp(cutoff, h_inv, G_min, l_max_zet, zet_d, C, err_d, para_env)
     367        21814 :          IF (err_c .GT. err_d) THEN
     368         2132 :             zet_b = zet_d
     369         2132 :             zet_d = zet_c
     370         2132 :             zet_c = zet_b - gr*(zet_b - zet_a)
     371              :          ELSE
     372        19682 :             zet_a = zet_c
     373        19682 :             zet_c = zet_d
     374        19682 :             zet_d = zet_a + gr*(zet_b - zet_a)
     375              :          END IF
     376              :       END DO
     377         1396 :       err_ctff = err_ctff_curr
     378              : 
     379         1396 :    END SUBROUTINE
     380              : 
     381              : ! **************************************************************************************************
     382              : !> \brief Calculate cutoff error estimate by using C_mm/3 * 1/(Gx^2*Gy^2*Gz^2)^1/3
     383              : !>        as an upper bound for 1/G^2 (and its minimax approximation) for |G| > G_c.
     384              : !>        Error is referring to a basis function P with fixed exponent zet_max and
     385              : !>        max. angular momentum l_max_zet.
     386              : !> \param cutoff ...
     387              : !> \param h_inv ...
     388              : !> \param G_min ...
     389              : !> \param l_max_zet ...
     390              : !> \param zet_max ...
     391              : !> \param C_mm ...
     392              : !> \param err_c ...
     393              : !> \param para_env ...
     394              : ! **************************************************************************************************
     395        49174 :    SUBROUTINE cutoff_error_fixed_exp(cutoff, h_inv, G_min, l_max_zet, zet_max, C_mm, err_c, para_env)
     396              :       REAL(KIND=dp), INTENT(IN)                          :: cutoff
     397              :       REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN)         :: h_inv
     398              :       REAL(KIND=dp), INTENT(IN)                          :: G_min
     399              :       INTEGER, INTENT(IN)                                :: l_max_zet
     400              :       REAL(KIND=dp), INTENT(IN)                          :: zet_max, C_mm
     401              :       REAL(KIND=dp), INTENT(OUT)                         :: err_c
     402              :       TYPE(mp_para_env_type), INTENT(IN)                 :: para_env
     403              : 
     404              :       INTEGER                                            :: ax, ay, az, G_l, G_u, Gl_first, Gl_last, &
     405              :                                                             Gu_first, Gu_last, i_xyz, l, my_p, &
     406              :                                                             n_Gl, n_Gl_left, n_Gl_p, n_Gu, &
     407              :                                                             n_Gu_left, n_Gu_p, n_p
     408              :       REAL(KIND=dp)                                      :: alpha_G, eps_G, err_c_l, G_c, G_rad, &
     409              :                                                             G_res, inv_lgth, prefactor, sum_G_diff
     410              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: S_G_l, S_G_u
     411              : 
     412        49174 :       G_c = SQRT(2.0_dp*cutoff)
     413        49174 :       eps_G = TINY(eps_G) ! sum up to machine precision
     414        49174 :       G_res = 0.5_dp*G_min ! resolution for screening
     415              : 
     416        49174 :       err_c = 0.0_dp
     417        49174 :       alpha_G = 1.0_dp/(2.0_dp*zet_max)
     418        49174 :       prefactor = 1.0_dp/zet_max
     419              : 
     420       196696 :       ALLOCATE (S_G_l(0:2*l_max_zet, 3))
     421        98348 :       ALLOCATE (S_G_u(0:2*l_max_zet, 3))
     422              : 
     423        49174 :       G_rad = exp_radius(2*l_max_zet, alpha_G, eps_G, prefactor, epsabs=G_res)
     424              : 
     425              :       ! Parallelization of sum over G vectors
     426        49174 :       my_p = para_env%mepos ! mpi rank
     427        49174 :       n_p = para_env%num_pe ! total number of processes
     428              : 
     429       196696 :       DO i_xyz = 1, 3
     430       147522 :          inv_lgth = ABS(h_inv(i_xyz, i_xyz))
     431              : 
     432       147522 :          G_l = FLOOR(G_c/(inv_lgth*twopi))
     433       147522 :          G_u = FLOOR(G_rad/(inv_lgth*twopi))
     434              : 
     435       147522 :          IF (G_u .LT. G_l) G_u = G_l
     436              : 
     437              :          ! Serial code:
     438              :          ! !Sum |G| <= G_c
     439              :          ! CALL pgf_sum_2c_gspace_1d_deltal(S_G_l(:, i_xyz), alpha_G, inv_lgth, -G_l, G_l, &
     440              :          !                               2.0_dp/3.0_dp, prefactor)
     441              :          ! !Sum |G| > G_c
     442              :          ! CALL pgf_sum_2c_gspace_1d_deltal(S_G_u(:, i_xyz), alpha_G, inv_lgth, G_l + 1, G_u, &
     443              :          !                               2.0_dp/3.0_dp, prefactor)
     444              : 
     445              :          ! Parallel code:
     446       147522 :          n_Gu = MAX((G_u - G_l), 0)
     447       147522 :          n_Gl = 2*G_l + 1
     448       147522 :          n_Gu_p = n_Gu/n_p
     449       147522 :          n_Gl_p = n_Gl/n_p
     450       147522 :          n_Gu_left = MOD(n_Gu, n_p)
     451       147522 :          n_Gl_left = MOD(n_Gl, n_p)
     452              : 
     453       147522 :          IF (my_p .LT. n_Gu_left) THEN
     454        34831 :             Gu_first = G_l + 1 + (n_Gu_p + 1)*my_p
     455        34831 :             Gu_last = G_l + 1 + (n_Gu_p + 1)*(my_p + 1) - 1
     456              :          ELSE
     457       112691 :             Gu_first = G_l + 1 + n_Gu_left + n_Gu_p*my_p
     458       112691 :             Gu_last = G_l + 1 + n_Gu_left + n_Gu_p*(my_p + 1) - 1
     459              :          END IF
     460              : 
     461       147522 :          IF (my_p .LT. n_Gl_left) THEN
     462        73761 :             Gl_first = -G_l + (n_Gl_p + 1)*my_p
     463        73761 :             Gl_last = -G_l + (n_Gl_p + 1)*(my_p + 1) - 1
     464              :          ELSE
     465        73761 :             Gl_first = -G_l + n_Gl_left + n_Gl_p*my_p
     466        73761 :             Gl_last = -G_l + n_Gl_left + n_Gl_p*(my_p + 1) - 1
     467              :          END IF
     468              : 
     469              :          ! Sum |G| <= G_c
     470              :          CALL pgf_sum_2c_gspace_1d_deltal(S_G_l(:, i_xyz), alpha_G, inv_lgth, Gl_first, Gl_last, &
     471       147522 :                                           2.0_dp/3.0_dp, prefactor)
     472              : 
     473              :          ! Sum |G| > G_c
     474              :          CALL pgf_sum_2c_gspace_1d_deltal(S_G_u(:, i_xyz), alpha_G, inv_lgth, Gu_first, Gu_last, &
     475       196696 :                                           2.0_dp/3.0_dp, prefactor)
     476              :       END DO
     477              : 
     478        49174 :       CALL para_env%sum(S_G_l)
     479        49174 :       CALL para_env%sum(S_G_u)
     480              : 
     481       879226 :       S_G_u = S_G_u*2.0_dp ! to include negative values of G
     482              : 
     483       187516 :       DO l = 0, l_max_zet
     484       482956 :       DO ax = 0, l
     485       994330 :       DO ay = 0, l - ax
     486       560548 :          az = l - ax - ay
     487              : 
     488              :          ! 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}
     489              :          ! Note: term by term multiplication to avoid subtraction for numerical stability
     490              :          sum_G_diff = S_G_u(2*ax, 1)*S_G_u(2*ay, 2)*S_G_u(2*az, 3) + &
     491              :                       S_G_u(2*ax, 1)*S_G_u(2*ay, 2)*S_G_l(2*az, 3) + &
     492              :                       S_G_u(2*ax, 1)*S_G_l(2*ay, 2)*S_G_u(2*az, 3) + &
     493              :                       S_G_l(2*ax, 1)*S_G_u(2*ay, 2)*S_G_u(2*az, 3) + &
     494              :                       S_G_u(2*ax, 1)*S_G_l(2*ay, 2)*S_G_l(2*az, 3) + &
     495              :                       S_G_l(2*ax, 1)*S_G_u(2*ay, 2)*S_G_l(2*az, 3) + &
     496       560548 :                       S_G_l(2*ax, 1)*S_G_l(2*ay, 2)*S_G_u(2*az, 3)
     497              : 
     498              :          err_c_l = 4.0_dp*pi**4*hermite_gauss_norm(zet_max, [ax, ay, az])**2* &
     499      2242192 :                    C_mm/3.0_dp*sum_G_diff
     500              : 
     501       855988 :          err_c = MAX(err_c, err_c_l)
     502              :       END DO
     503              :       END DO
     504              :       END DO
     505              : 
     506        49174 :       DEALLOCATE (S_G_u, S_G_l)
     507              : 
     508        49174 :    END SUBROUTINE cutoff_error_fixed_exp
     509              : 
     510              : END MODULE eri_mme_error_control
        

Generated by: LCOV version 2.0-1