LCOV - code coverage report
Current view: top level - src/eri_mme - eri_mme_integrate.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:d0bd076) Lines: 225 286 78.7 %
Date: 2021-09-15 13:52:28 Functions: 6 6 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 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       35426 :    SUBROUTINE eri_mme_2c_integrate(param, la_min, la_max, lb_min, lb_max, zeta, zetb, rab, &
      71       35426 :                                    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       35426 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: S_G_1, S_G_2, S_G_no, S_G_no_H
      93       35426 :       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       35426 :       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       35426 :                                G_rad=G_rad)
     103       35426 :       cutoff = G_rad**2/2
     104       35426 :       CALL get_minimax_from_cutoff(param%minimax_grid, cutoff, n_aw, aw, grid)
     105             : 
     106       35426 :       CPASSERT(grid .GT. 0)
     107             : 
     108             :       ! cell info
     109      460538 :       h_inv = param%h_inv
     110      460538 :       hmat = param%hmat
     111       35426 :       vol = param%vol
     112             : 
     113       35426 :       IF (PRESENT(normalize)) THEN
     114        2304 :          norm = normalize
     115             :       ELSE
     116             :          norm = .FALSE.
     117             :       ENDIF
     118             : 
     119       35426 :       l_max = la_max + lb_max
     120             : 
     121       35426 :       IF (PRESENT(potential)) THEN
     122             :          exact = .TRUE.
     123             :       ELSE
     124             :          exact = .FALSE.
     125             :       ENDIF
     126             : 
     127             :       IF (exact) THEN
     128         192 :          is_ortho = .FALSE.
     129             :       ELSE
     130       35234 :          is_ortho = param%is_ortho
     131             :       ENDIF
     132             : 
     133       35426 :       IF (is_ortho) THEN
     134      151620 :          ALLOCATE (S_G(0:l_max, 3, n_aw))
     135     6499020 :          S_G = 0.0_dp
     136             : 
     137       30324 :          IF (param%debug) THEN
     138           0 :             ALLOCATE (S_G_1(0:l_max))
     139           0 :             ALLOCATE (S_G_2(0:l_max))
     140             :          ENDIF
     141             :       ELSE
     142       15306 :          ALLOCATE (S_G_no(ncoset(l_max)))
     143      136955 :          S_G_no(:) = 0.0_dp
     144       15306 :          ALLOCATE (S_G_no_H(ncoset(l_max)))
     145             :       ENDIF
     146             : 
     147       35426 :       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      517178 :          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      481944 :                                        G_bounds, G_rad, R_bounds, R_rad)
     163      481944 :             alpha_G = aw(i_aw) + 0.25_dp/zeta + 0.25_dp/zetb
     164      481944 :             alpha_R = 0.25_dp/alpha_G
     165      517178 :             IF (is_ortho) THEN ! orthorhombic cell
     166             : 
     167             :                ! 1) precompute Ewald-like sum
     168             : 
     169     1802784 :                DO i_xyz = 1, 3
     170     1352088 :                   lgth = ABS(hmat(i_xyz, i_xyz))
     171     1352088 :                   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     1352088 :                   do_g_sum = n_sum_1d(1, i_xyz) < n_sum_1d(2, i_xyz) !G_bounds < R_bounds
     175             : 
     176     1352088 :                   IF (do_g_sum) THEN
     177      123584 :                      CALL pgf_sum_2c_gspace_1d(S_G(:, i_xyz, i_aw), -rab(i_xyz), alpha_G, inv_lgth, G_bounds(i_xyz))
     178      123584 :                      IF (PRESENT(G_count)) G_count = G_count + 1
     179             :                   ELSE
     180     1228504 :                      CALL pgf_sum_2c_rspace_1d(S_G(:, i_xyz, i_aw), -rab(i_xyz), alpha_R, lgth, R_bounds(i_xyz))
     181     1228504 :                      IF (PRESENT(R_count)) R_count = R_count + 1
     182             :                   ENDIF
     183             : 
     184     1802784 :                   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             :                   ENDIF
     192             :                ENDDO
     193             : 
     194             :             ELSE ! general cell
     195             : 
     196       31248 :                do_g_sum = n_sum_3d(1) < n_sum_3d(2) !PRODUCT(2*R_bounds+1) .GT. PRODUCT(2*G_bounds+1)
     197             : 
     198       31248 :                IF (do_g_sum) THEN
     199       20340 :                   CALL pgf_sum_2c_gspace_3d(S_G_no_H, l_max, -rab, alpha_G, h_inv, G_bounds, G_rad, vol)
     200        5085 :                   IF (PRESENT(G_count)) G_count = G_count + 1
     201             :                ELSE
     202      104652 :                   CALL pgf_sum_2c_rspace_3d(S_G_no_H, l_max, -rab, alpha_R, hmat, h_inv, R_bounds, R_rad)
     203       26163 :                   IF (PRESENT(R_count)) R_count = R_count + 1
     204             :                ENDIF
     205      472956 :                S_G_no(:) = S_G_no(:) + aw(n_aw + i_aw)*S_G_no_H
     206             :             ENDIF
     207             :          ENDDO
     208             :       ENDIF
     209             : 
     210             :       ! prefactor for integral values (unnormalized Hermite Gaussians)
     211       35426 :       prefac = SQRT(1.0_dp/(zeta*zetb))
     212             : 
     213             :       ! 2) Assemble integral values from Ewald sums
     214      208200 :       DO jco = ncoset(lb_min - 1) + 1, ncoset(lb_max)
     215      172774 :          CALL get_l(jco, lb, bx, by, bz)
     216      691096 :          lb_xyz = [bx, by, bz]
     217     3084426 :          DO ico = ncoset(la_min - 1) + 1, ncoset(la_max)
     218     2876226 :             CALL get_l(ico, la, ax, ay, az)
     219    11504904 :             la_xyz = [ax, ay, az]
     220     2876226 :             IF (is_ortho) THEN
     221     2188859 :                Imm = 0.0_dp
     222    35483391 :                DO i_aw = 1, n_aw
     223             :                   Ixyz = 1.0_dp
     224   133178128 :                   DO i_xyz = 1, 3
     225   133178128 :                      Ixyz = Ixyz*S_G(la_xyz(i_xyz) + lb_xyz(i_xyz), i_xyz, i_aw)*prefac
     226             :                   ENDDO
     227    35483391 :                   Imm = Imm + aw(n_aw + i_aw)*Ixyz
     228             :                ENDDO
     229             :             ELSE
     230      687367 :                Imm = S_G_no(coset(ax + bx, ay + by, az + bz))*prefac**3
     231             :             ENDIF
     232     2876226 :             IF (la + lb .EQ. 0 .AND. .NOT. exact) THEN
     233      149316 :                Imm = Imm - SUM(aw(n_aw + 1:2*n_aw))*prefac**3/vol ! subtracting G = 0 term
     234             :             ENDIF
     235     3049000 :             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      557058 :                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             :             ENDIF
     244             :          ENDDO ! la
     245             :       ENDDO ! lb
     246             : 
     247       70852 :    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      112434 :    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      112434 :                                    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      112434 :       CPASSERT(param%is_valid)
     284      112434 :       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       92034 :                                          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             :       ENDIF
     293      112434 :    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       92034 :    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       92034 :                                          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       92034 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: S_G_0
     333             :       REAL(KIND=dp), ALLOCATABLE, &
     334       92034 :          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       92034 :       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       92034 :                                param%sum_precision, G_rads_1=G_rads_1)
     345             : 
     346       92034 :       cutoff = (MIN(G_rads_1(1), G_rads_1(2) + G_rads_1(3)))**2/2
     347             : 
     348       92034 :       CALL get_minimax_from_cutoff(param%minimax_grid, cutoff, n_aw, aw, grid)
     349             : 
     350       92034 :       CPASSERT(grid .GT. 0)
     351             : 
     352             :       ! minimax coeffs
     353       92034 :       n_aw = param%minimax_grid(grid)%n_minimax
     354       92034 :       aw => param%minimax_grid(grid)%minimax_aw
     355             : 
     356             :       ! cell info
     357     1196442 :       h_inv = param%h_inv
     358     1196442 :       hmat = param%hmat
     359       92034 :       vol = param%vol
     360             : 
     361             :       ! prefactor for integral values (unnormalized Hermite Gaussians)
     362       92034 :       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       92034 :       G_res = 0.5_dp*param%G_min
     366       92034 :       R_res = 0.5_dp*param%R_min
     367             : 
     368      552204 :       ALLOCATE (S_G(n_aw, 3, 0:la_max, 0:lb_max, 0:lc_max))
     369             : 
     370             :       ! G= 0 components
     371       92034 :       IF (lc_min == 0) THEN
     372      165480 :          ALLOCATE (S_G_0(0:la_max + lb_max, 3))
     373             : 
     374       41370 :          alpha_R_0 = zeta*zetb/(zeta + zetb)
     375       41370 :          lmax_0 = la_max + lb_max
     376       41370 :          R_rad_0 = exp_radius(lmax_0, alpha_R_0, param%sum_precision, 1.0_dp, epsabs=R_res)
     377             : 
     378       41370 :          lgth = ABS(hmat(1, 1))
     379       41370 :          CALL pgf_sum_2c_rspace_1d(S_G_0(:, 1), RB(1) - RA(1), alpha_R_0, lgth, R_rad_0/lgth)
     380       41370 :          lgth = ABS(hmat(2, 2))
     381       41370 :          CALL pgf_sum_2c_rspace_1d(S_G_0(:, 2), RB(2) - RA(2), alpha_R_0, lgth, R_rad_0/lgth)
     382       41370 :          lgth = ABS(hmat(3, 3))
     383       41370 :          CALL pgf_sum_2c_rspace_1d(S_G_0(:, 3), RB(3) - RA(3), alpha_R_0, lgth, R_rad_0/lgth)
     384             :       ENDIF
     385             : 
     386     1357674 :       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     1265640 :                                     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     1265640 :                             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     1265640 :                             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     1265640 :                             R_bounds_3(:, 3))
     401             : 
     402     1357674 :          IF (PRESENT(RR_count)) RR_count = RR_count + 3
     403             :       ENDDO
     404             : 
     405             :       CALL eri_mme_3c_collect_ortho(vol, la_min, la_max, lb_min, lb_max, lc_min, lc_max, &
     406       92034 :                                     zeta, zetb, zetc, n_aw, aw, S_G, S_G_0, prefac, habc, o1, o2, o3)
     407             : 
     408      184068 :    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             :       ENDIF
     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             :       ENDIF
     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             :          ENDIF
     527             :       ENDIF
     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             :             ENDIF
     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             :             ENDIF
     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             :             ENDIF
     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             :             ENDIF
     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             :                      ENDIF
     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             :                      ENDIF
     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             :                      ENDIF
     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             :                   ENDDO
     629             :                ENDDO
     630             :             ENDDO
     631             : 
     632           0 :             CPASSERT(max_error .LE. param%debug_delta)
     633           0 :             DEALLOCATE (results_no)
     634             :          ENDIF
     635             :       ENDDO
     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       92034 :    SUBROUTINE eri_mme_3c_collect_ortho(vol, la_min, la_max, lb_min, lb_max, lc_min, lc_max, &
     667       92034 :                                        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       92034 :       resc_a_init = (2.0_dp*zeta)**la_min
     691       92034 :       resc_b_init = (2.0_dp*zetb)**lb_min
     692       92034 :       resc_c_init = (2.0_dp*zetc)**lc_min
     693             : 
     694       92034 :       resc_c = resc_c_init
     695       92034 :       lc_prev = lc_min
     696      625142 :       DO kco = ncoset(lc_min - 1) + 1, ncoset(lc_max)
     697      533108 :          CALL get_l(kco, lc, cx, cy, cz)
     698      533108 :          IF (lc_prev < lc) resc_c = resc_c*(2.0_dp*zetc)
     699             : 
     700      533108 :          resc_b = resc_b_init
     701      533108 :          lb_prev = lb_min
     702     4924248 :          DO jco = ncoset(lb_min - 1) + 1, ncoset(lb_max)
     703     4391140 :             CALL get_l(jco, lb, bx, by, bz)
     704     4391140 :             mone_lb = (-1.0_dp)**lb
     705     4391140 :             IF (lb_prev < lb) resc_b = resc_b*(2.0_dp*zetb)
     706             : 
     707     4391140 :             resc_a = resc_a_init
     708     4391140 :             la_prev = la_min
     709    39103168 :             DO ico = ncoset(la_min - 1) + 1, ncoset(la_max)
     710    34712028 :                CALL get_l(ico, la, ax, ay, az)
     711             : 
     712    34712028 :                IF (la_prev < la) resc_a = resc_a*(2.0_dp*zeta)
     713    34712028 :                Ixyz_0 = 0.0_dp
     714    34712028 :                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     2450822 :                            /vol*mone_lb
     719             :                ENDIF
     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   586497228 :                          S_G(1:n_aw, 3, az, bz, cz)) - Ixyz_0)
     724             : 
     725             :                ! rescaling needed due to Hermite Gaussians
     726    34712028 :                habc(o1 + ico, o2 + jco, o3 + kco) = Imm*prefac/(resc_a*resc_b*resc_c)
     727    39103168 :                la_prev = la
     728             :             ENDDO ! la
     729     4924248 :             lb_prev = lb
     730             :          ENDDO ! lb
     731      625142 :          lc_prev = lc
     732             :       ENDDO ! lc
     733             : 
     734       92034 :    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             :                ENDIF
     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             :             ENDDO ! la
     815      505890 :             lb_prev = lb
     816             :          ENDDO ! lb
     817      100350 :          lc_prev = lc
     818             :       ENDDO ! lc
     819             : 
     820       20400 :    END SUBROUTINE
     821             : 
     822             : END MODULE eri_mme_integrate

Generated by: LCOV version 1.15