LCOV - code coverage report
Current view: top level - src/eri_mme - eri_mme_error_control.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:d0bd076) Lines: 150 171 87.7 %
Date: 2021-09-15 13:52:28 Functions: 5 5 100.0 %

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

Generated by: LCOV version 1.15