LCOV - code coverage report
Current view: top level - src/eri_mme - eri_mme_integrate.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 78.7 % 286 225
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 6 6

            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 Minimax-Ewald (MME) method for calculating 2-center and 3-center
      10              : !>        electron repulsion integrals (ERI) of periodic systems using a
      11              : !>        Hermite Gaussian basis.
      12              : !>        The method relies on analytical Fourier transforms of Cartesian and
      13              : !>        Hermite Gaussian functions and Poisson summation formula to represent
      14              : !>        ERIs as a discrete sum over direct lattice vectors or reciprocal
      15              : !>        lattice vectors. The reciprocal space potential 1/G^2 is approximated
      16              : !>        by a linear combination of Gaussians employing minimax approximation.
      17              : !>        Not yet implemented: 3c ERIs for nonorthogonal cells.
      18              : !> \par History
      19              : !>       2015 09 created
      20              : !> \author Patrick Seewald
      21              : ! **************************************************************************************************
      22              : 
      23              : MODULE eri_mme_integrate
      24              :    USE ao_util,                         ONLY: exp_radius
      25              :    USE eri_mme_gaussian,                ONLY: hermite_gauss_norm
      26              :    USE eri_mme_lattice_summation,       ONLY: &
      27              :         ellipsoid_bounds, eri_mme_2c_get_bounds, eri_mme_2c_get_rads, eri_mme_3c_get_bounds, &
      28              :         eri_mme_3c_get_rads, get_l, pgf_sum_2c_gspace_1d, pgf_sum_2c_gspace_3d, &
      29              :         pgf_sum_2c_rspace_1d, pgf_sum_2c_rspace_3d, pgf_sum_3c_1d, pgf_sum_3c_3d
      30              :    USE eri_mme_types,                   ONLY: eri_mme_param,&
      31              :                                               get_minimax_from_cutoff
      32              :    USE kinds,                           ONLY: dp,&
      33              :                                               int_8
      34              :    USE mathconstants,                   ONLY: pi
      35              :    USE orbital_pointers,                ONLY: coset,&
      36              :                                               ncoset
      37              : #include "../base/base_uses.f90"
      38              : 
      39              :    IMPLICIT NONE
      40              : 
      41              :    PRIVATE
      42              : 
      43              :    LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .FALSE.
      44              : 
      45              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'eri_mme_integrate'
      46              : 
      47              :    PUBLIC :: eri_mme_2c_integrate, eri_mme_3c_integrate
      48              : 
      49              : CONTAINS
      50              : 
      51              : ! **************************************************************************************************
      52              : !> \brief Low-level integration routine for 2-center ERIs.
      53              : !> \param param ...
      54              : !> \param la_min ...
      55              : !> \param la_max ...
      56              : !> \param lb_min ...
      57              : !> \param lb_max ...
      58              : !> \param zeta ...
      59              : !> \param zetb ...
      60              : !> \param rab ...
      61              : !> \param hab ...
      62              : !> \param o1 ...
      63              : !> \param o2 ...
      64              : !> \param G_count ...
      65              : !> \param R_count ...
      66              : !> \param normalize calculate integrals w.r.t. normalized Hermite-Gaussians
      67              : !> \param potential use exact potential instead of minimax approx. (for testing only)
      68              : !> \param pot_par potential parameter
      69              : ! **************************************************************************************************
      70        55934 :    SUBROUTINE eri_mme_2c_integrate(param, la_min, la_max, lb_min, lb_max, zeta, zetb, rab, &
      71        55934 :                                    hab, o1, o2, G_count, R_count, normalize, potential, pot_par)
      72              :       TYPE(eri_mme_param), INTENT(IN)                    :: param
      73              :       INTEGER, INTENT(IN)                                :: la_min, la_max, lb_min, lb_max
      74              :       REAL(KIND=dp), INTENT(IN)                          :: zeta, zetb
      75              :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rab
      76              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT)      :: hab
      77              :       INTEGER, INTENT(IN)                                :: o1, o2
      78              :       INTEGER, INTENT(INOUT), OPTIONAL                   :: G_count, R_count
      79              :       LOGICAL, INTENT(IN), OPTIONAL                      :: normalize
      80              :       INTEGER, INTENT(IN), OPTIONAL                      :: potential
      81              :       REAL(KIND=dp), INTENT(IN), OPTIONAL                :: pot_par
      82              : 
      83              :       INTEGER                                            :: ax, ay, az, bx, by, bz, grid, i_aw, &
      84              :                                                             i_xyz, ico, jco, l_max, la, lb, n_aw
      85              :       INTEGER(KIND=int_8), DIMENSION(2)                  :: n_sum_3d
      86              :       INTEGER(KIND=int_8), DIMENSION(2, 3)               :: n_sum_1d
      87              :       INTEGER, DIMENSION(3)                              :: la_xyz, lb_xyz
      88              :       LOGICAL                                            :: do_g_sum, exact, is_ortho, norm
      89              :       REAL(KIND=dp)                                      :: alpha_G, alpha_R, cutoff, G_rad, G_res, &
      90              :                                                             Imm, inv_lgth, Ixyz, lgth, max_error, &
      91              :                                                             prefac, R_rad, R_res, vol
      92        55934 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: S_G_1, S_G_2, S_G_no, S_G_no_H
      93        55934 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: S_G
      94              :       REAL(KIND=dp), DIMENSION(3)                        :: G_bounds, R_bounds
      95              :       REAL(KIND=dp), DIMENSION(3, 3)                     :: h_inv, hmat
      96        55934 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: aw
      97              : 
      98            0 :       CPASSERT(param%is_valid)
      99              : 
     100              :       grid = 0
     101              :       CALL eri_mme_2c_get_rads(la_max, lb_max, zeta, zetb, 0.0_dp, param%G_min, param%R_min, param%sum_precision, &
     102        55934 :                                G_rad=G_rad)
     103        55934 :       cutoff = G_rad**2/2
     104        55934 :       CALL get_minimax_from_cutoff(param%minimax_grid, cutoff, n_aw, aw, grid)
     105              : 
     106        55934 :       CPASSERT(grid .GT. 0)
     107              : 
     108              :       ! cell info
     109       727142 :       h_inv = param%h_inv
     110       727142 :       hmat = param%hmat
     111        55934 :       vol = param%vol
     112              : 
     113        55934 :       IF (PRESENT(normalize)) THEN
     114         2304 :          norm = normalize
     115              :       ELSE
     116              :          norm = .FALSE.
     117              :       END IF
     118              : 
     119        55934 :       l_max = la_max + lb_max
     120              : 
     121        55934 :       IF (PRESENT(potential)) THEN
     122              :          exact = .TRUE.
     123              :       ELSE
     124              :          exact = .FALSE.
     125              :       END IF
     126              : 
     127              :       IF (exact) THEN
     128          192 :          is_ortho = .FALSE.
     129              :       ELSE
     130        55742 :          is_ortho = param%is_ortho
     131              :       END IF
     132              : 
     133        55934 :       IF (is_ortho) THEN
     134       225320 :          ALLOCATE (S_G(0:l_max, 3, n_aw))
     135     10584840 :          S_G = 0.0_dp
     136              : 
     137        45064 :          IF (param%debug) THEN
     138            0 :             ALLOCATE (S_G_1(0:l_max))
     139            0 :             ALLOCATE (S_G_2(0:l_max))
     140              :          END IF
     141              :       ELSE
     142        32610 :          ALLOCATE (S_G_no(ncoset(l_max)))
     143       242259 :          S_G_no(:) = 0.0_dp
     144        21740 :          ALLOCATE (S_G_no_H(ncoset(l_max)))
     145              :       END IF
     146              : 
     147        55934 :       IF (exact) THEN
     148          192 :          alpha_G = 0.25_dp/zeta + 0.25_dp/zetb
     149              :          ! resolution for Gaussian width
     150          192 :          G_res = 0.5_dp*param%G_min
     151          192 :          R_res = 0.5_dp*param%R_min
     152              : 
     153          192 :          G_rad = exp_radius(l_max, alpha_G, 0.01*param%sum_precision, 1.0_dp, epsabs=G_res)
     154         2496 :          G_bounds(:) = ellipsoid_bounds(G_rad, TRANSPOSE(hmat)/(2.0_dp*pi))
     155          768 :          CALL pgf_sum_2c_gspace_3d(S_G_no, l_max, -rab, alpha_G, h_inv, G_bounds, G_rad, vol, potential, pot_par)
     156              :       ELSE
     157              : 
     158       947846 :          DO i_aw = 1, n_aw
     159              : 
     160              :             CALL eri_mme_2c_get_bounds(hmat, h_inv, vol, is_ortho, param%G_min, param%R_min, la_max, lb_max, &
     161              :                                        zeta, zetb, aw(i_aw), param%sum_precision, n_sum_1d, n_sum_3d, &
     162       892104 :                                        G_bounds, G_rad, R_bounds, R_rad)
     163       892104 :             alpha_G = aw(i_aw) + 0.25_dp/zeta + 0.25_dp/zetb
     164       892104 :             alpha_R = 0.25_dp/alpha_G
     165       947846 :             IF (is_ortho) THEN ! orthorhombic cell
     166              : 
     167              :                ! 1) precompute Ewald-like sum
     168              : 
     169      2981984 :                DO i_xyz = 1, 3
     170      2236488 :                   lgth = ABS(hmat(i_xyz, i_xyz))
     171      2236488 :                   inv_lgth = ABS(h_inv(i_xyz, i_xyz))
     172              : 
     173              :                   ! perform sum in R or G space. Choose the space in which less summands are required for convergence
     174      2236488 :                   do_g_sum = n_sum_1d(1, i_xyz) < n_sum_1d(2, i_xyz) !G_bounds < R_bounds
     175              : 
     176      2236488 :                   IF (do_g_sum) THEN
     177       202238 :                      CALL pgf_sum_2c_gspace_1d(S_G(:, i_xyz, i_aw), -rab(i_xyz), alpha_G, inv_lgth, G_bounds(i_xyz))
     178       202238 :                      IF (PRESENT(G_count)) G_count = G_count + 1
     179              :                   ELSE
     180      2034250 :                      CALL pgf_sum_2c_rspace_1d(S_G(:, i_xyz, i_aw), -rab(i_xyz), alpha_R, lgth, R_bounds(i_xyz))
     181      2034250 :                      IF (PRESENT(R_count)) R_count = R_count + 1
     182              :                   END IF
     183              : 
     184      2981984 :                   IF (param%debug) THEN
     185              :                      ! check consistency of summation methods
     186            0 :                      CALL pgf_sum_2c_gspace_1d(S_G_1, -rab(i_xyz), alpha_G, inv_lgth, G_bounds(i_xyz))
     187            0 :                      CALL pgf_sum_2c_rspace_1d(S_G_2, -rab(i_xyz), alpha_R, lgth, R_bounds(i_xyz))
     188            0 :                      max_error = MAXVAL(ABS(S_G_1 - S_G_2)/(0.5_dp*(ABS(S_G_1) + ABS(S_G_2)) + 1.0_dp))
     189              : 
     190            0 :                      CPASSERT(max_error .LE. param%debug_delta)
     191              :                   END IF
     192              :                END DO
     193              : 
     194              :             ELSE ! general cell
     195              : 
     196       146608 :                do_g_sum = n_sum_3d(1) < n_sum_3d(2) !PRODUCT(2*R_bounds+1) .GT. PRODUCT(2*G_bounds+1)
     197              : 
     198       146608 :                IF (do_g_sum) THEN
     199       109556 :                   CALL pgf_sum_2c_gspace_3d(S_G_no_H, l_max, -rab, alpha_G, h_inv, G_bounds, G_rad, vol)
     200        27389 :                   IF (PRESENT(G_count)) G_count = G_count + 1
     201              :                ELSE
     202       476876 :                   CALL pgf_sum_2c_rspace_3d(S_G_no_H, l_max, -rab, alpha_R, hmat, h_inv, R_bounds, R_rad)
     203       119219 :                   IF (PRESENT(R_count)) R_count = R_count + 1
     204              :                END IF
     205      2579036 :                S_G_no(:) = S_G_no(:) + aw(n_aw + i_aw)*S_G_no_H
     206              :             END IF
     207              :          END DO
     208              :       END IF
     209              : 
     210              :       ! prefactor for integral values (unnormalized Hermite Gaussians)
     211        55934 :       prefac = SQRT(1.0_dp/(zeta*zetb))
     212              : 
     213              :       ! 2) Assemble integral values from Ewald sums
     214       294428 :       DO jco = ncoset(lb_min - 1) + 1, ncoset(lb_max)
     215       238494 :          CALL get_l(jco, lb, bx, by, bz)
     216       953976 :          lb_xyz = [bx, by, bz]
     217      3577034 :          DO ico = ncoset(la_min - 1) + 1, ncoset(la_max)
     218      3282606 :             CALL get_l(ico, la, ax, ay, az)
     219     13130424 :             la_xyz = [ax, ay, az]
     220      3282606 :             IF (is_ortho) THEN
     221      2481319 :                Imm = 0.0_dp
     222     41625051 :                DO i_aw = 1, n_aw
     223              :                   Ixyz = 1.0_dp
     224    156574928 :                   DO i_xyz = 1, 3
     225    156574928 :                      Ixyz = Ixyz*S_G(la_xyz(i_xyz) + lb_xyz(i_xyz), i_xyz, i_aw)*prefac
     226              :                   END DO
     227     41625051 :                   Imm = Imm + aw(n_aw + i_aw)*Ixyz
     228              :                END DO
     229              :             ELSE
     230       801287 :                Imm = S_G_no(coset(ax + bx, ay + by, az + bz))*prefac**3
     231              :             END IF
     232      3282606 :             IF (la + lb .EQ. 0 .AND. .NOT. exact) THEN
     233       226050 :                Imm = Imm - SUM(aw(n_aw + 1:2*n_aw))*prefac**3/vol ! subtracting G = 0 term
     234              :             END IF
     235      3521100 :             IF (.NOT. norm) THEN
     236              :                ! rescaling needed due to Hermite Gaussians (such that they can be contracted same way as Cartesian Gaussians)
     237              :                ! and factor of 4 pi**4 (-1)**lb
     238       963438 :                hab(o1 + ico, o2 + jco) = Imm*4.0_dp*pi**4/((2.0_dp*zeta)**la*(-2.0_dp*zetb)**lb)
     239              :             ELSE
     240              :                ! same thing for normalized Hermite Gaussians
     241              :                hab(o1 + ico, o2 + jco) = Imm*4.0_dp*pi**4*(-1.0_dp)**lb*hermite_gauss_norm(zeta, la_xyz)* &
     242      2319168 :                                          hermite_gauss_norm(zetb, lb_xyz)
     243              :             END IF
     244              :          END DO ! la
     245              :       END DO ! lb
     246              : 
     247       111868 :    END SUBROUTINE eri_mme_2c_integrate
     248              : 
     249              : ! **************************************************************************************************
     250              : !> \brief Low-level integration routine for 3-center ERIs
     251              : !> \param param ...
     252              : !> \param la_min ...
     253              : !> \param la_max ...
     254              : !> \param lb_min ...
     255              : !> \param lb_max ...
     256              : !> \param lc_min ...
     257              : !> \param lc_max ...
     258              : !> \param zeta ...
     259              : !> \param zetb ...
     260              : !> \param zetc ...
     261              : !> \param RA ...
     262              : !> \param RB ...
     263              : !> \param RC ...
     264              : !> \param habc ...
     265              : !> \param o1 ...
     266              : !> \param o2 ...
     267              : !> \param o3 ...
     268              : !> \param GG_count ...
     269              : !> \param GR_count ...
     270              : !> \param RR_count ...
     271              : ! **************************************************************************************************
     272       154554 :    SUBROUTINE eri_mme_3c_integrate(param, la_min, la_max, lb_min, lb_max, lc_min, lc_max, zeta, zetb, zetc, RA, RB, RC, &
     273       154554 :                                    habc, o1, o2, o3, GG_count, GR_count, RR_count)
     274              :       TYPE(eri_mme_param), INTENT(IN)                    :: param
     275              :       INTEGER, INTENT(IN)                                :: la_min, la_max, lb_min, lb_max, lc_min, &
     276              :                                                             lc_max
     277              :       REAL(KIND=dp), INTENT(IN)                          :: zeta, zetb, zetc
     278              :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: RA, RB, RC
     279              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT)   :: habc
     280              :       INTEGER, INTENT(IN)                                :: o1, o2, o3
     281              :       INTEGER, INTENT(INOUT), OPTIONAL                   :: GG_count, GR_count, RR_count
     282              : 
     283       154554 :       CPASSERT(param%is_valid)
     284       154554 :       IF (param%is_ortho) THEN
     285              :          CALL eri_mme_3c_integrate_ortho(param, la_min, la_max, lb_min, lb_max, lc_min, lc_max, zeta, zetb, zetc, RA, RB, RC, &
     286       134154 :                                          habc, o1, o2, o3, RR_count)
     287              : 
     288              :       ELSE
     289              :          CALL eri_mme_3c_integrate_nonortho(param, la_min, la_max, lb_min, lb_max, lc_min, lc_max, zeta, zetb, zetc, RA, RB, RC, &
     290        20400 :                                             habc, o1, o2, o3, GG_count, GR_count, RR_count)
     291              : 
     292              :       END IF
     293       154554 :    END SUBROUTINE eri_mme_3c_integrate
     294              : 
     295              : ! **************************************************************************************************
     296              : !> \brief ...
     297              : !> \param param ...
     298              : !> \param la_min ...
     299              : !> \param la_max ...
     300              : !> \param lb_min ...
     301              : !> \param lb_max ...
     302              : !> \param lc_min ...
     303              : !> \param lc_max ...
     304              : !> \param zeta ...
     305              : !> \param zetb ...
     306              : !> \param zetc ...
     307              : !> \param RA ...
     308              : !> \param RB ...
     309              : !> \param RC ...
     310              : !> \param habc ...
     311              : !> \param o1 ...
     312              : !> \param o2 ...
     313              : !> \param o3 ...
     314              : !> \param RR_count ...
     315              : ! **************************************************************************************************
     316       134154 :    SUBROUTINE eri_mme_3c_integrate_ortho(param, la_min, la_max, lb_min, lb_max, lc_min, lc_max, zeta, zetb, zetc, RA, RB, RC, &
     317       134154 :                                          habc, o1, o2, o3, RR_count)
     318              :       TYPE(eri_mme_param), INTENT(IN)                    :: param
     319              :       INTEGER, INTENT(IN)                                :: la_min, la_max, lb_min, lb_max, lc_min, &
     320              :                                                             lc_max
     321              :       REAL(KIND=dp), INTENT(IN)                          :: zeta, zetb, zetc
     322              :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: RA, RB, RC
     323              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT)   :: habc
     324              :       INTEGER, INTENT(IN)                                :: o1, o2, o3
     325              :       INTEGER, INTENT(INOUT), OPTIONAL                   :: RR_count
     326              : 
     327              :       INTEGER                                            :: grid, i_aw, lmax_0, n_aw
     328              :       INTEGER(KIND=int_8), DIMENSION(3)                  :: n_sum_3d
     329              :       INTEGER(KIND=int_8), DIMENSION(3, 3)               :: n_sum_1d
     330              :       REAL(KIND=dp)                                      :: alpha_R_0, cutoff, G_res, lgth, prefac, &
     331              :                                                             R_rad_0, R_res, vol
     332       134154 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: S_G_0
     333              :       REAL(KIND=dp), ALLOCATABLE, &
     334       134154 :          DIMENSION(:, :, :, :, :)                        :: S_G
     335              :       REAL(KIND=dp), DIMENSION(2)                        :: R_rads_3
     336              :       REAL(KIND=dp), DIMENSION(2, 3)                     :: R_bounds_3
     337              :       REAL(KIND=dp), DIMENSION(3)                        :: G_rads_1, R_rads_2
     338              :       REAL(KIND=dp), DIMENSION(3, 3)                     :: G_bounds_1, h_inv, hmat, R_bounds_2
     339       134154 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: aw
     340              : 
     341              :       grid = 0
     342              : 
     343              :       CALL eri_mme_3c_get_rads(la_max, lb_max, lc_max, zeta, zetb, zetc, 0.0_dp, param%G_min, param%R_min, &
     344       134154 :                                param%sum_precision, G_rads_1=G_rads_1)
     345              : 
     346       134154 :       cutoff = (MIN(G_rads_1(1), G_rads_1(2) + G_rads_1(3)))**2/2
     347              : 
     348       134154 :       CALL get_minimax_from_cutoff(param%minimax_grid, cutoff, n_aw, aw, grid)
     349              : 
     350       134154 :       CPASSERT(grid .GT. 0)
     351              : 
     352              :       ! minimax coeffs
     353       134154 :       n_aw = param%minimax_grid(grid)%n_minimax
     354       134154 :       aw => param%minimax_grid(grid)%minimax_aw
     355              : 
     356              :       ! cell info
     357      1744002 :       h_inv = param%h_inv
     358      1744002 :       hmat = param%hmat
     359       134154 :       vol = param%vol
     360              : 
     361              :       ! prefactor for integral values (unnormalized Hermite Gaussians)
     362       134154 :       prefac = (zeta*zetb*zetc)**(-1.5_dp)*pi**(11.0_dp/2.0_dp)*4.0_dp
     363              : 
     364              :       ! Preparations for G=0 component
     365       134154 :       G_res = 0.5_dp*param%G_min
     366       134154 :       R_res = 0.5_dp*param%R_min
     367              : 
     368       804924 :       ALLOCATE (S_G(n_aw, 3, 0:la_max, 0:lb_max, 0:lc_max))
     369              : 
     370              :       ! G= 0 components
     371       134154 :       IF (lc_min == 0) THEN
     372       237912 :          ALLOCATE (S_G_0(0:la_max + lb_max, 3))
     373              : 
     374        59478 :          alpha_R_0 = zeta*zetb/(zeta + zetb)
     375        59478 :          lmax_0 = la_max + lb_max
     376        59478 :          R_rad_0 = exp_radius(lmax_0, alpha_R_0, param%sum_precision, 1.0_dp, epsabs=R_res)
     377              : 
     378        59478 :          lgth = ABS(hmat(1, 1))
     379        59478 :          CALL pgf_sum_2c_rspace_1d(S_G_0(:, 1), RB(1) - RA(1), alpha_R_0, lgth, R_rad_0/lgth)
     380        59478 :          lgth = ABS(hmat(2, 2))
     381        59478 :          CALL pgf_sum_2c_rspace_1d(S_G_0(:, 2), RB(2) - RA(2), alpha_R_0, lgth, R_rad_0/lgth)
     382        59478 :          lgth = ABS(hmat(3, 3))
     383        59478 :          CALL pgf_sum_2c_rspace_1d(S_G_0(:, 3), RB(3) - RA(3), alpha_R_0, lgth, R_rad_0/lgth)
     384              :       END IF
     385              : 
     386      2242194 :       DO i_aw = 1, n_aw
     387              :          CALL eri_mme_3c_get_bounds(hmat, h_inv, vol, .TRUE., param%G_min, param%R_min, la_max, lb_max, lc_max, &
     388              :                                     zeta, zetb, zetc, aw(i_aw), param%sum_precision, n_sum_1d, n_sum_3d, &
     389      2108040 :                                     G_bounds_1, G_rads_1, R_bounds_2, R_rads_2, R_bounds_3, R_rads_3)
     390              :          CALL pgf_sum_3c_1d(S_G(i_aw, 1, :, :, :), RA(1), RB(1), RC(1), &
     391              :                             zeta, zetb, zetc, aw(i_aw), ABS(hmat(1, 1)), &
     392      2108040 :                             R_bounds_3(:, 1))
     393              : 
     394              :          CALL pgf_sum_3c_1d(S_G(i_aw, 2, :, :, :), RA(2), RB(2), RC(2), &
     395              :                             zeta, zetb, zetc, aw(i_aw), ABS(hmat(2, 2)), &
     396      2108040 :                             R_bounds_3(:, 2))
     397              : 
     398              :          CALL pgf_sum_3c_1d(S_G(i_aw, 3, :, :, :), RA(3), RB(3), RC(3), &
     399              :                             zeta, zetb, zetc, aw(i_aw), ABS(hmat(3, 3)), &
     400      2108040 :                             R_bounds_3(:, 3))
     401              : 
     402      2242194 :          IF (PRESENT(RR_count)) RR_count = RR_count + 3
     403              :       END DO
     404              : 
     405              :       CALL eri_mme_3c_collect_ortho(vol, la_min, la_max, lb_min, lb_max, lc_min, lc_max, &
     406       134154 :                                     zeta, zetb, zetc, n_aw, aw, S_G, S_G_0, prefac, habc, o1, o2, o3)
     407              : 
     408       268308 :    END SUBROUTINE
     409              : 
     410              : ! **************************************************************************************************
     411              : !> \brief ...
     412              : !> \param param ...
     413              : !> \param la_min ...
     414              : !> \param la_max ...
     415              : !> \param lb_min ...
     416              : !> \param lb_max ...
     417              : !> \param lc_min ...
     418              : !> \param lc_max ...
     419              : !> \param zeta ...
     420              : !> \param zetb ...
     421              : !> \param zetc ...
     422              : !> \param RA ...
     423              : !> \param RB ...
     424              : !> \param RC ...
     425              : !> \param habc ...
     426              : !> \param o1 ...
     427              : !> \param o2 ...
     428              : !> \param o3 ...
     429              : !> \param GG_count ...
     430              : !> \param GR_count ...
     431              : !> \param RR_count ...
     432              : ! **************************************************************************************************
     433        20400 :    SUBROUTINE eri_mme_3c_integrate_nonortho(param, la_min, la_max, lb_min, lb_max, lc_min, lc_max, zeta, zetb, zetc, RA, RB, RC, &
     434        20400 :                                             habc, o1, o2, o3, GG_count, GR_count, RR_count)
     435              : 
     436              :       TYPE(eri_mme_param), INTENT(IN)                    :: param
     437              :       INTEGER, INTENT(IN)                                :: la_min, la_max, lb_min, lb_max, lc_min, &
     438              :                                                             lc_max
     439              :       REAL(KIND=dp), INTENT(IN)                          :: zeta, zetb, zetc
     440              :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: RA, RB, RC
     441              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT)   :: habc
     442              :       INTEGER, INTENT(IN)                                :: o1, o2, o3
     443              :       INTEGER, INTENT(INOUT), OPTIONAL                   :: GG_count, GR_count, RR_count
     444              : 
     445              :       INTEGER                                            :: ax, ay, az, bx, by, bz, cx, cy, cz, &
     446              :                                                             grid, i_aw, ico, ir, jco, kco, la, lb, &
     447              :                                                             lc, lmax_0, method, n_aw, nresults, &
     448              :                                                             sum_method
     449              :       INTEGER(KIND=int_8), DIMENSION(3)                  :: n_sum_3d
     450              :       INTEGER(KIND=int_8), DIMENSION(3, 3)               :: n_sum_1d
     451              :       LOGICAL                                            :: db_sum1, db_sum2, db_sum3, do_g_sum_0
     452              :       REAL(KIND=dp)                                      :: alpha_G_0, alpha_R_0, cutoff, G_rad_0, &
     453              :                                                             G_res, max_error, max_result, &
     454              :                                                             min_result, prefac, R_rad_0, R_res, vol
     455        20400 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: results_no, S_G_0_no
     456        20400 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: S_G_no, S_G_no_1, S_G_no_2, S_G_no_3, &
     457        20400 :                                                             S_G_no_H
     458              :       REAL(KIND=dp), DIMENSION(2)                        :: R_rads_3
     459              :       REAL(KIND=dp), DIMENSION(2, 3)                     :: R_bounds_3
     460              :       REAL(KIND=dp), DIMENSION(3)                        :: G_bound_0, G_rads_1, R_0, R_bound_0, &
     461              :                                                             R_rads_2
     462              :       REAL(KIND=dp), DIMENSION(3, 3)                     :: G_bounds_1, h_inv, hmat, R_bounds_2
     463        20400 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: aw
     464              : 
     465            0 :       CPASSERT(param%is_valid)
     466              : 
     467              :       grid = 0
     468              : 
     469              :       CALL eri_mme_3c_get_rads(la_max, lb_max, lc_max, zeta, zetb, zetc, 0.0_dp, param%G_min, param%R_min, &
     470        20400 :                                param%sum_precision, G_rads_1=G_rads_1)
     471              : 
     472        20400 :       cutoff = (MIN(G_rads_1(1), G_rads_1(2) + G_rads_1(3)))**2/2
     473              : 
     474        20400 :       CALL get_minimax_from_cutoff(param%minimax_grid, cutoff, n_aw, aw, grid)
     475              : 
     476        20400 :       CPASSERT(grid .GT. 0)
     477              : 
     478              :       ! minimax coeffs
     479        20400 :       n_aw = param%minimax_grid(grid)%n_minimax
     480        20400 :       aw => param%minimax_grid(grid)%minimax_aw
     481              : 
     482              :       ! cell info
     483       265200 :       h_inv = param%h_inv
     484       265200 :       hmat = param%hmat
     485        20400 :       vol = param%vol
     486              : 
     487              :       ! prefactor for integral values (unnormalized Hermite Gaussians)
     488        20400 :       prefac = (zeta*zetb*zetc)**(-1.5_dp)*pi**(11.0_dp/2.0_dp)*4.0_dp
     489              : 
     490        20400 :       IF (param%debug) THEN
     491            0 :          ALLOCATE (S_G_no_1(ncoset(la_max), ncoset(lb_max), ncoset(lc_max)))
     492            0 :          ALLOCATE (S_G_no_2(ncoset(la_max), ncoset(lb_max), ncoset(lc_max)))
     493            0 :          ALLOCATE (S_G_no_3(ncoset(la_max), ncoset(lb_max), ncoset(lc_max)))
     494              :       END IF
     495              : 
     496              :       ! Preparations for G=0 component
     497        20400 :       G_res = 0.5_dp*param%G_min
     498        20400 :       R_res = 0.5_dp*param%R_min
     499              : 
     500       102000 :       ALLOCATE (S_G_no(ncoset(la_max), ncoset(lb_max), ncoset(lc_max)))
     501              : 
     502      6330636 :       S_G_no(:, :, :) = 0.0_dp
     503        20400 :       IF (param%debug) THEN
     504            0 :          S_G_no_1(:, :, :) = -1.0_dp
     505            0 :          S_G_no_2(:, :, :) = -1.0_dp
     506            0 :          S_G_no_3(:, :, :) = -1.0_dp
     507              :       END IF
     508       102000 :       ALLOCATE (S_G_no_H(ncoset(la_max), ncoset(lb_max), ncoset(lc_max)))
     509              : 
     510              :       ! G= 0 components
     511        20400 :       IF (lc_min == 0) THEN
     512        27000 :          ALLOCATE (S_G_0_no(ncoset(la_max + lb_max)))
     513         9000 :          alpha_G_0 = 0.25_dp/zetb + 0.25_dp/zeta
     514         9000 :          alpha_R_0 = 0.25_dp/alpha_G_0
     515         9000 :          lmax_0 = la_max + lb_max
     516        36000 :          R_0 = RB - RA
     517         9000 :          G_rad_0 = exp_radius(lmax_0, alpha_G_0, param%sum_precision, 1.0_dp, epsabs=G_res)
     518         9000 :          R_rad_0 = exp_radius(lmax_0, alpha_R_0, param%sum_precision, 1.0_dp, epsabs=R_res)
     519       117000 :          G_bound_0 = ellipsoid_bounds(G_rad_0, TRANSPOSE(hmat)/(2.0_dp*pi))
     520         9000 :          R_bound_0 = ellipsoid_bounds(R_rad_0, h_inv)
     521        63000 :          do_g_sum_0 = PRODUCT(2*R_bound_0 + 1) .GT. PRODUCT(2*G_bound_0 + 1)
     522         9000 :          IF (do_g_sum_0) THEN
     523            0 :             CALL pgf_sum_2c_gspace_3d(S_G_0_no, lmax_0, R_0, alpha_G_0, h_inv, G_bound_0, G_rad_0, vol)
     524              :          ELSE
     525         9000 :             CALL pgf_sum_2c_rspace_3d(S_G_0_no, lmax_0, R_0, alpha_R_0, hmat, h_inv, R_bound_0, R_rad_0)
     526              :          END IF
     527              :       END IF
     528              : 
     529       168900 :       DO i_aw = 1, n_aw
     530              :          CALL eri_mme_3c_get_bounds(hmat, h_inv, vol, .FALSE., param%G_min, param%R_min, la_max, lb_max, lc_max, &
     531              :                                     zeta, zetb, zetc, aw(i_aw), param%sum_precision, n_sum_1d, n_sum_3d, &
     532       148500 :                                     G_bounds_1, G_rads_1, R_bounds_2, R_rads_2, R_bounds_3, R_rads_3)
     533       594000 :          sum_method = MINLOC(n_sum_3d, DIM=1)
     534              : 
     535              :          CALL pgf_sum_3c_3d(S_G_no_H, la_max, lb_max, lc_max, RA, RB, RC, &
     536              :                             zeta, zetb, zetc, aw(i_aw), hmat, h_inv, vol, &
     537              :                             G_bounds_1, R_bounds_2, R_bounds_3, &
     538              :                             G_rads_1, R_rads_2, R_rads_3, &
     539       148500 :                             method=sum_method, method_out=method)
     540     24569280 :          S_G_no(:, :, :) = S_G_no(:, :, :) + aw(n_aw + i_aw)*S_G_no_H(:, :, :)
     541              : 
     542            0 :          SELECT CASE (method)
     543              :          CASE (1)
     544            0 :             IF (PRESENT(GG_count)) GG_count = GG_count + 1
     545              :          CASE (2)
     546         5846 :             IF (PRESENT(GR_count)) GR_count = GR_count + 1
     547              :          CASE (3)
     548       142654 :             IF (PRESENT(RR_count)) RR_count = RR_count + 1
     549              :          CASE DEFAULT
     550       148500 :             CPABORT("")
     551              :          END SELECT
     552              : 
     553       317400 :          IF (param%debug) THEN
     554            0 :             nresults = 0
     555              :             ! check consistency of summation methods
     556              : 
     557            0 :             db_sum1 = (n_sum_3d(1)) .LT. INT(param%debug_nsum, KIND=int_8)
     558            0 :             db_sum2 = (n_sum_3d(2)) .LT. INT(param%debug_nsum, KIND=int_8)
     559            0 :             db_sum3 = (n_sum_3d(3)) .LT. INT(param%debug_nsum, KIND=int_8)
     560              : 
     561            0 :             IF (param%unit_nr > 0) THEN
     562            0 :                WRITE (param%unit_nr, *) "ERI_MME DEBUG | number of summands (GG / GR / RR)", n_sum_3d
     563            0 :                WRITE (param%unit_nr, *) "ERI_MME DEBUG | sum methods to be compared (GG / GR / RR)", db_sum1, db_sum2, db_sum3
     564              :             END IF
     565              : 
     566            0 :             S_G_no_1(:, :, :) = 0.0_dp
     567            0 :             S_G_no_2(:, :, :) = 0.0_dp
     568            0 :             S_G_no_3(:, :, :) = 0.0_dp
     569              : 
     570            0 :             IF (db_sum1) THEN
     571              :                CALL pgf_sum_3c_3d(S_G_no_1, la_max, lb_max, lc_max, RA, RB, RC, &
     572              :                                   zeta, zetb, zetc, aw(i_aw), hmat, h_inv, vol, &
     573              :                                   G_bounds_1, R_bounds_2, R_bounds_3, &
     574              :                                   G_rads_1, R_rads_2, R_rads_3, &
     575            0 :                                   method=1)
     576            0 :                nresults = nresults + 1
     577              :             END IF
     578              : 
     579            0 :             IF (db_sum2) THEN
     580              :                CALL pgf_sum_3c_3d(S_G_no_2, la_max, lb_max, lc_max, RA, RB, RC, &
     581              :                                   zeta, zetb, zetc, aw(i_aw), hmat, h_inv, vol, &
     582              :                                   G_bounds_1, R_bounds_2, R_bounds_3, &
     583              :                                   G_rads_1, R_rads_2, R_rads_3, &
     584            0 :                                   method=2)
     585            0 :                nresults = nresults + 1
     586              :             END IF
     587              : 
     588            0 :             IF (db_sum3) THEN
     589              :                CALL pgf_sum_3c_3d(S_G_no_3, la_max, lb_max, lc_max, RA, RB, RC, &
     590              :                                   zeta, zetb, zetc, aw(i_aw), hmat, h_inv, vol, &
     591              :                                   G_bounds_1, R_bounds_2, R_bounds_3, &
     592              :                                   G_rads_1, R_rads_2, R_rads_3, &
     593            0 :                                   method=3)
     594            0 :                nresults = nresults + 1
     595              :             END IF
     596              : 
     597            0 :             max_error = 0.0_dp
     598            0 :             ALLOCATE (results_no(nresults))
     599              : 
     600            0 :             DO kco = ncoset(lc_min - 1) + 1, ncoset(lc_max)
     601            0 :                CALL get_l(kco, lc, cx, cy, cz)
     602            0 :                DO jco = ncoset(lb_min - 1) + 1, ncoset(lb_max)
     603            0 :                   CALL get_l(jco, lb, bx, by, bz)
     604            0 :                   DO ico = ncoset(la_min - 1) + 1, ncoset(la_max)
     605            0 :                      CALL get_l(ico, la, ax, ay, az)
     606              : 
     607            0 :                      max_error = 0.0_dp
     608            0 :                      ir = 0
     609            0 :                      IF (db_sum1) THEN
     610            0 :                         ir = ir + 1
     611            0 :                         results_no(ir) = S_G_no_1(ico, jco, kco)
     612              :                      END IF
     613              : 
     614            0 :                      IF (db_sum2) THEN
     615            0 :                         ir = ir + 1
     616            0 :                         results_no(ir) = S_G_no_2(ico, jco, kco)
     617              :                      END IF
     618              : 
     619            0 :                      IF (db_sum3) THEN
     620            0 :                         ir = ir + 1
     621            0 :                         results_no(ir) = S_G_no_3(ico, jco, kco)
     622              :                      END IF
     623              : 
     624            0 :                      max_result = MAXVAL(results_no)
     625            0 :                      min_result = MINVAL(results_no)
     626            0 :                      IF (nresults > 0) max_error = MAX(max_error, &
     627            0 :                                                     (max_result - min_result)/(0.5_dp*(ABS(max_result) + ABS(min_result)) + 1.0_dp))
     628              :                   END DO
     629              :                END DO
     630              :             END DO
     631              : 
     632            0 :             CPASSERT(max_error .LE. param%debug_delta)
     633            0 :             DEALLOCATE (results_no)
     634              :          END IF
     635              :       END DO
     636              : 
     637              :       ! assemble integral values
     638              : 
     639              :       CALL eri_mme_3c_collect_nonortho(vol, la_min, la_max, lb_min, lb_max, lc_min, lc_max, &
     640        20400 :                                        zeta, zetb, zetc, n_aw, aw, S_G_no, S_G_0_no, prefac, habc, o1, o2, o3)
     641              : 
     642        61200 :    END SUBROUTINE
     643              : 
     644              : ! **************************************************************************************************
     645              : !> \brief ...
     646              : !> \param vol ...
     647              : !> \param la_min ...
     648              : !> \param la_max ...
     649              : !> \param lb_min ...
     650              : !> \param lb_max ...
     651              : !> \param lc_min ...
     652              : !> \param lc_max ...
     653              : !> \param zeta ...
     654              : !> \param zetb ...
     655              : !> \param zetc ...
     656              : !> \param n_aw ...
     657              : !> \param aw ...
     658              : !> \param S_G ...
     659              : !> \param S_G_0 ...
     660              : !> \param prefac ...
     661              : !> \param habc ...
     662              : !> \param o1 ...
     663              : !> \param o2 ...
     664              : !> \param o3 ...
     665              : ! **************************************************************************************************
     666       134154 :    SUBROUTINE eri_mme_3c_collect_ortho(vol, la_min, la_max, lb_min, lb_max, lc_min, lc_max, &
     667       134154 :                                        zeta, zetb, zetc, n_aw, aw, S_G, S_G_0, prefac, habc, o1, o2, o3)
     668              :       REAL(KIND=dp), INTENT(IN)                          :: vol
     669              :       INTEGER, INTENT(IN)                                :: la_min, la_max, lb_min, lb_max, lc_min, &
     670              :                                                             lc_max
     671              :       REAL(KIND=dp), INTENT(IN)                          :: zeta, zetb, zetc
     672              :       INTEGER, INTENT(IN)                                :: n_aw
     673              :       REAL(KIND=dp), DIMENSION(2*n_aw), INTENT(IN)       :: aw
     674              :       REAL(KIND=dp), DIMENSION(:, :, 0:, 0:, 0:), &
     675              :          INTENT(IN)                                      :: S_G
     676              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
     677              :          INTENT(IN)                                      :: S_G_0
     678              :       REAL(KIND=dp), INTENT(IN)                          :: prefac
     679              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT)   :: habc
     680              :       INTEGER, INTENT(IN)                                :: o1, o2, o3
     681              : 
     682              :       INTEGER                                            :: ax, ay, az, bx, by, bz, cx, cy, cz, ico, &
     683              :                                                             jco, kco, la, la_prev, lb, lb_prev, &
     684              :                                                             lc, lc_prev
     685              :       REAL(KIND=dp)                                      :: Imm, Ixyz_0, mone_lb, resc_a, &
     686              :                                                             resc_a_init, resc_b, resc_b_init, &
     687              :                                                             resc_c, resc_c_init
     688              : 
     689              :       ! Initialization of rescaling factors due to Hermite Gaussians
     690       134154 :       resc_a_init = (2.0_dp*zeta)**la_min
     691       134154 :       resc_b_init = (2.0_dp*zetb)**lb_min
     692       134154 :       resc_c_init = (2.0_dp*zetc)**lc_min
     693              : 
     694       134154 :       resc_c = resc_c_init
     695       134154 :       lc_prev = lc_min
     696       903872 :       DO kco = ncoset(lc_min - 1) + 1, ncoset(lc_max)
     697       769718 :          CALL get_l(kco, lc, cx, cy, cz)
     698       769718 :          IF (lc_prev < lc) resc_c = resc_c*(2.0_dp*zetc)
     699              : 
     700       769718 :          resc_b = resc_b_init
     701       769718 :          lb_prev = lb_min
     702      6336708 :          DO jco = ncoset(lb_min - 1) + 1, ncoset(lb_max)
     703      5566990 :             CALL get_l(jco, lb, bx, by, bz)
     704      5566990 :             mone_lb = (-1.0_dp)**lb
     705      5566990 :             IF (lb_prev < lb) resc_b = resc_b*(2.0_dp*zetb)
     706              : 
     707      5566990 :             resc_a = resc_a_init
     708      5566990 :             la_prev = la_min
     709     47679670 :             DO ico = ncoset(la_min - 1) + 1, ncoset(la_max)
     710     42112680 :                CALL get_l(ico, la, ax, ay, az)
     711              : 
     712     42112680 :                IF (la_prev < la) resc_a = resc_a*(2.0_dp*zeta)
     713     42112680 :                Ixyz_0 = 0.0_dp
     714     42112680 :                IF (lc == 0) THEN
     715              :                   Ixyz_0 = S_G_0(ax + bx, 1)* &
     716              :                            S_G_0(ay + by, 2)* &
     717              :                            S_G_0(az + bz, 3) &
     718      2781458 :                            /vol*mone_lb
     719              :                END IF
     720              :                Imm = SUM(aw(n_aw + 1:2*n_aw)*( &
     721              :                          S_G(1:n_aw, 1, ax, bx, cx)* &
     722              :                          S_G(1:n_aw, 2, ay, by, cy)* &
     723    741910920 :                          S_G(1:n_aw, 3, az, bz, cz)) - Ixyz_0)
     724              : 
     725              :                ! rescaling needed due to Hermite Gaussians
     726     42112680 :                habc(o1 + ico, o2 + jco, o3 + kco) = Imm*prefac/(resc_a*resc_b*resc_c)
     727     47679670 :                la_prev = la
     728              :             END DO ! la
     729      6336708 :             lb_prev = lb
     730              :          END DO ! lb
     731       903872 :          lc_prev = lc
     732              :       END DO ! lc
     733              : 
     734       134154 :    END SUBROUTINE
     735              : 
     736              : ! **************************************************************************************************
     737              : !> \brief ...
     738              : !> \param vol ...
     739              : !> \param la_min ...
     740              : !> \param la_max ...
     741              : !> \param lb_min ...
     742              : !> \param lb_max ...
     743              : !> \param lc_min ...
     744              : !> \param lc_max ...
     745              : !> \param zeta ...
     746              : !> \param zetb ...
     747              : !> \param zetc ...
     748              : !> \param n_aw ...
     749              : !> \param aw ...
     750              : !> \param S_G ...
     751              : !> \param S_G_0 ...
     752              : !> \param prefac ...
     753              : !> \param habc ...
     754              : !> \param o1 ...
     755              : !> \param o2 ...
     756              : !> \param o3 ...
     757              : ! **************************************************************************************************
     758        20400 :    PURE SUBROUTINE eri_mme_3c_collect_nonortho(vol, la_min, la_max, lb_min, lb_max, lc_min, lc_max, &
     759        20400 :                                                zeta, zetb, zetc, n_aw, aw, S_G, S_G_0, prefac, habc, o1, o2, o3)
     760              :       REAL(KIND=dp), INTENT(IN)                          :: vol
     761              :       INTEGER, INTENT(IN)                                :: la_min, la_max, lb_min, lb_max, lc_min, &
     762              :                                                             lc_max
     763              :       REAL(KIND=dp), INTENT(IN)                          :: zeta, zetb, zetc
     764              :       INTEGER, INTENT(IN)                                :: n_aw
     765              :       REAL(KIND=dp), DIMENSION(2*n_aw), INTENT(IN)       :: aw
     766              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN)      :: S_G
     767              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
     768              :          INTENT(IN)                                      :: S_G_0
     769              :       REAL(KIND=dp), INTENT(IN)                          :: prefac
     770              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT)   :: habc
     771              :       INTEGER, INTENT(IN)                                :: o1, o2, o3
     772              : 
     773              :       INTEGER                                            :: ax, ay, az, bx, by, bz, cx, cy, cz, ico, &
     774              :                                                             ijco, jco, kco, la, la_prev, lb, &
     775              :                                                             lb_prev, lc, lc_prev
     776              :       REAL(KIND=dp)                                      :: Imm, mone_lb, resc_a, resc_a_init, &
     777              :                                                             resc_b, resc_b_init, resc_c, &
     778              :                                                             resc_c_init
     779              : 
     780              :       ! Initialization of rescaling factors due to Hermite Gaussians
     781        20400 :       resc_a_init = (2.0_dp*zeta)**la_min
     782        20400 :       resc_b_init = (2.0_dp*zetb)**lb_min
     783        20400 :       resc_c_init = (2.0_dp*zetc)**lc_min
     784              : 
     785        20400 :       resc_c = resc_c_init
     786        20400 :       lc_prev = lc_min
     787       100350 :       DO kco = ncoset(lc_min - 1) + 1, ncoset(lc_max)
     788        79950 :          CALL get_l(kco, lc, cx, cy, cz)
     789        79950 :          IF (lc_prev < lc) resc_c = resc_c*(2.0_dp*zetc)
     790              : 
     791        79950 :          resc_b = resc_b_init
     792        79950 :          lb_prev = lb_min
     793       505890 :          DO jco = ncoset(lb_min - 1) + 1, ncoset(lb_max)
     794       425940 :             CALL get_l(jco, lb, bx, by, bz)
     795       425940 :             mone_lb = (-1.0_dp)**lb
     796       425940 :             IF (lb_prev < lb) resc_b = resc_b*(2.0_dp*zetb)
     797              : 
     798       425940 :             resc_a = resc_a_init
     799       425940 :             la_prev = la_min
     800      4209543 :             DO ico = ncoset(la_min - 1) + 1, ncoset(la_max)
     801      3783603 :                CALL get_l(ico, la, ax, ay, az)
     802              : 
     803      3783603 :                IF (la_prev < la) resc_a = resc_a*(2.0_dp*zeta)
     804      3783603 :                IF (lc .GT. 0) THEN
     805      3631431 :                   Imm = S_G(ico, jco, kco)
     806              :                ELSE
     807       152172 :                   ijco = coset(ax + bx, ay + by, az + bz)
     808       810828 :                   Imm = S_G(ico, jco, kco) - SUM(aw(n_aw + 1:2*n_aw))*S_G_0(ijco)/vol*mone_lb
     809              :                END IF
     810              : 
     811              :                ! rescaling needed due to Hermite Gaussians
     812      3783603 :                habc(o1 + ico, o2 + jco, o3 + kco) = Imm*prefac/(resc_a*resc_b*resc_c)
     813      4209543 :                la_prev = la
     814              :             END DO ! la
     815       505890 :             lb_prev = lb
     816              :          END DO ! lb
     817       100350 :          lc_prev = lc
     818              :       END DO ! lc
     819              : 
     820        20400 :    END SUBROUTINE
     821              : 
     822              : END MODULE eri_mme_integrate
        

Generated by: LCOV version 2.0-1