LCOV - code coverage report
Current view: top level - src/eri_mme - eri_mme_lattice_summation.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:d0bd076) Lines: 543 681 79.7 %
Date: 2021-09-15 13:52:28 Functions: 152 281 54.1 %

          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 Ewald sums to represent integrals in direct and reciprocal lattice.
      10             : !> \par History
      11             : !>       2015 09 created
      12             : !> \author Patrick Seewald
      13             : ! **************************************************************************************************
      14             : 
      15             : MODULE eri_mme_lattice_summation
      16             : #:include "eri_mme_unroll.fypp"
      17             : 
      18             :    USE ao_util, ONLY: exp_radius
      19             :    USE eri_mme_gaussian, ONLY: create_gaussian_overlap_dist_to_hermite, &
      20             :                                create_hermite_to_cartesian, &
      21             :                                eri_mme_coulomb, eri_mme_yukawa, &
      22             :                                eri_mme_longrange
      23             :    USE kinds, ONLY: dp, &
      24             :                     int_8
      25             :    USE mathconstants, ONLY: gaussi, &
      26             :                             pi, &
      27             :                             twopi
      28             :    USE orbital_pointers, ONLY: coset, &
      29             :                                indco, &
      30             :                                ncoset
      31             : #include "../base/base_uses.f90"
      32             : 
      33             :    IMPLICIT NONE
      34             : 
      35             :    PRIVATE
      36             : 
      37             :    LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .FALSE.
      38             : 
      39             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'eri_mme_lattice_summation'
      40             : 
      41             :    PUBLIC :: &
      42             :       ellipsoid_bounds, &
      43             :       eri_mme_2c_get_bounds, &
      44             :       eri_mme_2c_get_rads, &
      45             :       eri_mme_3c_get_bounds, &
      46             :       eri_mme_3c_get_rads, &
      47             :       get_l, &
      48             :       pgf_sum_2c_gspace_1d, &
      49             :       pgf_sum_2c_gspace_1d_deltal, &
      50             :       pgf_sum_2c_gspace_3d, &
      51             :       pgf_sum_2c_rspace_1d, &
      52             :       pgf_sum_2c_rspace_3d, &
      53             :       pgf_sum_3c_1d, &
      54             :       pgf_sum_3c_3d
      55             : 
      56             :    INTEGER, PARAMETER, PRIVATE :: exp_w = 50, div_w = 10
      57             : 
      58             : CONTAINS
      59             : 
      60             : ! **************************************************************************************************
      61             : !> \brief Get summation radii for 2c integrals
      62             : !> \param la_max ...
      63             : !> \param lb_max ...
      64             : !> \param zeta ...
      65             : !> \param zetb ...
      66             : !> \param a_mm ...
      67             : !> \param G_min ...
      68             : !> \param R_min ...
      69             : !> \param sum_precision ...
      70             : !> \param G_rad ...
      71             : !> \param R_rad ...
      72             : ! **************************************************************************************************
      73      517370 :    SUBROUTINE eri_mme_2c_get_rads(la_max, lb_max, zeta, zetb, a_mm, G_min, R_min, sum_precision, G_rad, R_rad)
      74             :       INTEGER, INTENT(IN)                                :: la_max, lb_max
      75             :       REAL(KIND=dp), INTENT(IN)                          :: zeta, zetb, a_mm, G_min, R_min, &
      76             :                                                             sum_precision
      77             :       REAL(KIND=dp), INTENT(OUT), OPTIONAL               :: G_rad, R_rad
      78             : 
      79             :       INTEGER                                            :: l_max
      80             :       REAL(KIND=dp)                                      :: alpha_G, alpha_R, G_res, R_res
      81             : 
      82      517370 :       l_max = la_max + lb_max
      83      517370 :       alpha_G = a_mm + 0.25_dp/zeta + 0.25_dp/zetb
      84      517370 :       alpha_R = 0.25_dp/alpha_G
      85             : 
      86      517370 :       G_res = 0.5_dp*G_min
      87      517370 :       R_res = 0.5_dp*R_min
      88             : 
      89      517370 :       IF (PRESENT(G_rad)) G_rad = exp_radius(l_max, alpha_G, sum_precision, 1.0_dp, epsabs=G_res)
      90      517370 :       IF (PRESENT(R_rad)) R_rad = exp_radius(l_max, alpha_R, sum_precision, 1.0_dp, epsabs=R_res)
      91             : 
      92      517370 :    END SUBROUTINE
      93             : 
      94             : ! **************************************************************************************************
      95             : !> \brief Get summation radii for 3c integrals
      96             : !> \param la_max ...
      97             : !> \param lb_max ...
      98             : !> \param lc_max ...
      99             : !> \param zeta ...
     100             : !> \param zetb ...
     101             : !> \param zetc ...
     102             : !> \param a_mm ...
     103             : !> \param G_min ...
     104             : !> \param R_min ...
     105             : !> \param sum_precision ...
     106             : !> \param G_rads_1 ...
     107             : !> \param R_rads_2 ...
     108             : !> \param R_rads_3 ...
     109             : ! **************************************************************************************************
     110     1526574 :    SUBROUTINE eri_mme_3c_get_rads(la_max, lb_max, lc_max, zeta, zetb, zetc, a_mm, G_min, R_min, &
     111             :                                   sum_precision, G_rads_1, R_rads_2, R_rads_3)
     112             :       INTEGER, INTENT(IN)                                :: la_max, lb_max, lc_max
     113             :       REAL(KIND=dp), INTENT(IN)                          :: zeta, zetb, zetc, a_mm
     114             :       REAL(KIND=dp), INTENT(IN), OPTIONAL                :: G_min, R_min
     115             :       REAL(KIND=dp), INTENT(IN)                          :: sum_precision
     116             :       REAL(KIND=dp), DIMENSION(3), INTENT(OUT), OPTIONAL :: G_rads_1, R_rads_2
     117             :       REAL(KIND=dp), DIMENSION(2), INTENT(OUT), OPTIONAL :: R_rads_3
     118             : 
     119             :       REAL(KIND=dp)                                      :: alpha, alpha_R, beta, G_res, gamma, R_res
     120             : 
     121             :       ! exponents in G space
     122     1526574 :       alpha = 0.25_dp/zeta
     123     1526574 :       beta = 0.25_dp/zetb
     124     1526574 :       gamma = 0.25_dp/zetc + a_mm
     125             : 
     126             :       ! Summation radii and number of summands for all lattice summation methods
     127             :       ! sum method 1
     128     1526574 :       IF (PRESENT(G_rads_1)) THEN
     129     1526574 :          G_res = 0.5_dp*G_min
     130     1526574 :          G_rads_1(1) = exp_radius(la_max, alpha, sum_precision, 1.0_dp, G_res)
     131     1526574 :          G_rads_1(2) = exp_radius(lb_max, beta, sum_precision, 1.0_dp, G_res)
     132     1526574 :          G_rads_1(3) = exp_radius(lc_max, gamma, sum_precision, 1.0_dp, G_res)
     133             :       ENDIF
     134             : 
     135             :       ! sum method 2
     136     1526574 :       IF (PRESENT(R_rads_2)) THEN
     137     1414140 :          R_res = 0.5_dp*R_min
     138     1414140 :          R_rads_2(1) = exp_radius(lb_max + lc_max, 0.25_dp/(beta + gamma), sum_precision, 1.0_dp, epsabs=R_res)
     139     1414140 :          R_rads_2(2) = exp_radius(lc_max + la_max, 0.25_dp/(alpha + gamma), sum_precision, 1.0_dp, epsabs=R_res)
     140     1414140 :          R_rads_2(3) = exp_radius(lb_max + la_max, 0.25_dp/(alpha + beta), sum_precision, 1.0_dp, epsabs=R_res)
     141             :       ENDIF
     142             : 
     143             :       ! sum method 3
     144             : 
     145     1526574 :       IF (PRESENT(R_rads_3)) THEN
     146     1414140 :          R_res = 0.5_dp*R_min
     147     1414140 :          alpha_R = 1.0_dp/((zeta + zetb + zetc)/((zeta + zetb)*zetc) + 4.0_dp*a_mm)
     148     1414140 :          R_rads_3(1) = exp_radius(la_max + lb_max, zeta*zetb/(zeta + zetb), sum_precision, 1.0_dp, R_res)
     149     1414140 :          R_rads_3(2) = exp_radius(la_max + lb_max + lc_max, alpha_R, sum_precision, 1.0_dp, R_res)
     150             :       ENDIF
     151             : 
     152     1526574 :    END SUBROUTINE
     153             : 
     154             : ! **************************************************************************************************
     155             : !> \brief Get summation bounds for 2c integrals
     156             : !> \param hmat ...
     157             : !> \param h_inv ...
     158             : !> \param vol ...
     159             : !> \param is_ortho ...
     160             : !> \param G_min ...
     161             : !> \param R_min ...
     162             : !> \param la_max ...
     163             : !> \param lb_max ...
     164             : !> \param zeta ...
     165             : !> \param zetb ...
     166             : !> \param a_mm ...
     167             : !> \param sum_precision ...
     168             : !> \param n_sum_1d ...
     169             : !> \param n_sum_3d ...
     170             : !> \param G_bounds ...
     171             : !> \param G_rad ...
     172             : !> \param R_bounds ...
     173             : !> \param R_rad ...
     174             : ! **************************************************************************************************
     175      481944 :    SUBROUTINE eri_mme_2c_get_bounds(hmat, h_inv, vol, is_ortho, G_min, R_min, la_max, lb_max, &
     176             :                                     zeta, zetb, a_mm, sum_precision, n_sum_1d, n_sum_3d, &
     177             :                                     G_bounds, G_rad, R_bounds, R_rad)
     178             :       REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN)         :: hmat, h_inv
     179             :       REAL(KIND=dp), INTENT(IN)                          :: vol
     180             :       LOGICAL, INTENT(IN)                                :: is_ortho
     181             :       REAL(KIND=dp), INTENT(IN)                          :: G_min, R_min
     182             :       INTEGER, INTENT(IN)                                :: la_max, lb_max
     183             :       REAL(KIND=dp), INTENT(IN)                          :: zeta, zetb, a_mm, sum_precision
     184             :       INTEGER(KIND=int_8), DIMENSION(2, 3), INTENT(OUT)  :: n_sum_1d
     185             :       INTEGER(KIND=int_8), DIMENSION(2), INTENT(OUT)     :: n_sum_3d
     186             :       REAL(KIND=dp), DIMENSION(3), INTENT(OUT)           :: G_bounds
     187             :       REAL(KIND=dp), INTENT(OUT)                         :: G_rad
     188             :       REAL(KIND=dp), DIMENSION(3), INTENT(OUT)           :: R_bounds
     189             :       REAL(KIND=dp), INTENT(OUT)                         :: R_rad
     190             : 
     191             :       INTEGER                                            :: i_xyz
     192             :       REAL(KIND=dp)                                      :: ns_G, ns_R, vol_G
     193             : 
     194      481944 :       CALL eri_mme_2c_get_rads(la_max, lb_max, zeta, zetb, a_mm, G_min, R_min, sum_precision, G_rad, R_rad)
     195             : 
     196     6265272 :       G_bounds = ellipsoid_bounds(G_rad, TRANSPOSE(hmat)/(2.0_dp*pi))
     197      481944 :       R_bounds = ellipsoid_bounds(R_rad, h_inv)
     198             : 
     199      481944 :       vol_G = twopi**3/vol
     200             : 
     201      481944 :       IF (is_ortho) THEN
     202     1802784 :          DO i_xyz = 1, 3
     203     1352088 :             ns_G = 2.0_dp*G_bounds(i_xyz)
     204     1352088 :             ns_R = 2.0_dp*R_bounds(i_xyz)
     205     1352088 :             n_sum_1d(1, i_xyz) = nsum_2c_gspace_1d(ns_G, la_max, lb_max)
     206     1802784 :             n_sum_1d(2, i_xyz) = nsum_2c_rspace_1d(ns_R, la_max, lb_max)
     207             :          ENDDO
     208             :       ELSE
     209       31248 :          ns_G = 4._dp/3._dp*pi*G_rad**3/vol_G
     210       31248 :          ns_R = 4._dp/3._dp*pi*R_rad**3/vol
     211       31248 :          n_sum_3d(1) = nsum_2c_gspace_3d(ns_G, la_max, lb_max)
     212       31248 :          n_sum_3d(2) = nsum_2c_rspace_3d(ns_R, la_max, lb_max)
     213             :       ENDIF
     214             : 
     215      481944 :    END SUBROUTINE
     216             : 
     217             : ! **************************************************************************************************
     218             : !> \brief Get summation bounds for 3c integrals
     219             : !> \param hmat ...
     220             : !> \param h_inv ...
     221             : !> \param vol ...
     222             : !> \param is_ortho ...
     223             : !> \param G_min ...
     224             : !> \param R_min ...
     225             : !> \param la_max ...
     226             : !> \param lb_max ...
     227             : !> \param lc_max ...
     228             : !> \param zeta ...
     229             : !> \param zetb ...
     230             : !> \param zetc ...
     231             : !> \param a_mm ...
     232             : !> \param sum_precision ...
     233             : !> \param n_sum_1d ...
     234             : !> \param n_sum_3d ...
     235             : !> \param G_bounds_1 ...
     236             : !> \param G_rads_1 ...
     237             : !> \param R_bounds_2 ...
     238             : !> \param R_rads_2 ...
     239             : !> \param R_bounds_3 ...
     240             : !> \param R_rads_3 ...
     241             : ! **************************************************************************************************
     242     1414140 :    SUBROUTINE eri_mme_3c_get_bounds(hmat, h_inv, vol, is_ortho, G_min, R_min, la_max, lb_max, lc_max, &
     243             :                                     zeta, zetb, zetc, a_mm, sum_precision, n_sum_1d, n_sum_3d, &
     244             :                                     G_bounds_1, G_rads_1, R_bounds_2, R_rads_2, R_bounds_3, R_rads_3)
     245             :       REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN)         :: hmat, h_inv
     246             :       REAL(KIND=dp), INTENT(IN)                          :: vol
     247             :       LOGICAL, INTENT(IN)                                :: is_ortho
     248             :       REAL(KIND=dp), INTENT(IN)                          :: G_min, R_min
     249             :       INTEGER, INTENT(IN)                                :: la_max, lb_max, lc_max
     250             :       REAL(KIND=dp), INTENT(IN)                          :: zeta, zetb, zetc, a_mm, sum_precision
     251             :       INTEGER(KIND=int_8), DIMENSION(3, 3), INTENT(OUT)  :: n_sum_1d
     252             :       INTEGER(KIND=int_8), DIMENSION(3), INTENT(OUT)     :: n_sum_3d
     253             :       REAL(KIND=dp), DIMENSION(3, 3)                     :: G_bounds_1
     254             :       REAL(KIND=dp), DIMENSION(3), INTENT(OUT)           :: G_rads_1
     255             :       REAL(KIND=dp), DIMENSION(3, 3)                     :: R_bounds_2
     256             :       REAL(KIND=dp), DIMENSION(3), INTENT(OUT)           :: R_rads_2
     257             :       REAL(KIND=dp), DIMENSION(2, 3)                     :: R_bounds_3
     258             :       REAL(KIND=dp), DIMENSION(2), INTENT(OUT)           :: R_rads_3
     259             : 
     260             :       INTEGER                                            :: i, i_xyz, order_1, order_2
     261             :       REAL(KIND=dp)                                      :: ns1_G1, ns1_G2, ns2_G, ns2_R, ns3_R1, &
     262             :                                                             ns3_R2, vol_G
     263             :       CALL eri_mme_3c_get_rads(la_max, lb_max, lc_max, zeta, zetb, zetc, a_mm, G_min, R_min, sum_precision, &
     264     1414140 :                                G_rads_1=G_rads_1, R_rads_2=R_rads_2, R_rads_3=R_rads_3)
     265             : 
     266     1414140 :       vol_G = twopi**3/vol
     267             : 
     268     5656560 :       order_1 = MAXLOC(G_rads_1, DIM=1)
     269     5656560 :       order_2 = MINLOC(G_rads_1, DIM=1)
     270             : 
     271     5656560 :       DO i = 1, 3
     272    69292860 :          G_bounds_1(i, :) = ellipsoid_bounds(G_rads_1(i), TRANSPOSE(hmat)/(2.0_dp*pi))
     273             :       ENDDO
     274             : 
     275     5656560 :       DO i = 1, 3
     276    18383820 :          R_bounds_2(i, :) = ellipsoid_bounds(R_rads_2(i), h_inv)
     277             :       ENDDO
     278             : 
     279     4242420 :       DO i = 1, 2
     280    12727260 :          R_bounds_3(i, :) = ellipsoid_bounds(R_rads_3(i), h_inv)
     281             :       ENDDO
     282             : 
     283     1414140 :       IF (is_ortho) THEN
     284     5062560 :          DO i_xyz = 1, 3
     285             : 
     286     3796920 :             ns3_R1 = 2.0_dp*R_bounds_3(1, i_xyz)
     287     3796920 :             ns3_R2 = 2.0_dp*R_bounds_3(2, i_xyz)
     288             : 
     289     5062560 :             n_sum_1d(3, i_xyz) = nsum_3c_rspace_1d(ns3_R1, ns3_R2)
     290             :          ENDDO
     291             : 
     292             :       ELSE
     293             : 
     294      594000 :          order_1 = MAXLOC(G_rads_1, DIM=1)
     295      594000 :          order_2 = MINLOC(G_rads_1, DIM=1)
     296             : 
     297       71062 :          SELECT CASE (order_1)
     298             :          CASE (1)
     299       71062 :             ns1_G1 = 4._dp/3._dp*pi*G_rads_1(2)**3/vol_G
     300       71062 :             ns1_G2 = 4._dp/3._dp*pi*G_rads_1(3)**3/vol_G
     301             :          CASE (2)
     302       53482 :             ns1_G1 = 4._dp/3._dp*pi*G_rads_1(1)**3/vol_G
     303       53482 :             ns1_G2 = 4._dp/3._dp*pi*G_rads_1(3)**3/vol_G
     304             :          CASE (3)
     305       23956 :             ns1_G1 = 4._dp/3._dp*pi*G_rads_1(1)**3/vol_G
     306      148500 :             ns1_G2 = 4._dp/3._dp*pi*G_rads_1(2)**3/vol_G
     307             :          END SELECT
     308             : 
     309      148500 :          n_sum_3d(1) = nsum_3c_gspace_3d(ns1_G1, ns1_G2, la_max, lb_max, lc_max)
     310             : 
     311      148500 :          ns2_G = 4._dp/3._dp*pi*G_rads_1(order_2)**3/vol_G
     312      148500 :          ns2_R = 4._dp/3._dp*pi*R_rads_2(order_2)**3/vol
     313             : 
     314      148500 :          ns3_R1 = 4._dp/3._dp*pi*R_rads_3(1)**3/vol
     315      148500 :          ns3_R2 = 4._dp/3._dp*pi*R_rads_3(2)**3/vol
     316             : 
     317       30492 :          SELECT CASE (order_2)
     318             :          CASE (1)
     319       30492 :             n_sum_3d(2) = nsum_product_3c_gspace_3d(ns2_G, ns2_R, la_max, lb_max, lc_max)
     320             :          CASE (2)
     321       24372 :             n_sum_3d(2) = nsum_product_3c_gspace_3d(ns2_G, ns2_R, lb_max, la_max, lc_max)
     322             :          CASE (3)
     323      148500 :             n_sum_3d(2) = nsum_product_3c_gspace_3d(ns2_G, ns2_R, lc_max, lb_max, la_max)
     324             :          END SELECT
     325             : 
     326      148500 :          n_sum_3d(3) = nsum_3c_rspace_3d(ns3_R1, ns3_R2, la_max, lb_max, lc_max)
     327             :       ENDIF
     328             : 
     329     1414140 :    END SUBROUTINE
     330             : 
     331             : ! **************************************************************************************************
     332             : !> \brief Roughly estimated number of floating point operations
     333             : !> \param ns_G ...
     334             : !> \param l ...
     335             : !> \param m ...
     336             : !> \return ...
     337             : ! **************************************************************************************************
     338     1352088 :    PURE FUNCTION nsum_2c_gspace_1d(ns_G, l, m)
     339             :       REAL(KIND=dp), INTENT(IN)                          :: ns_G
     340             :       INTEGER, INTENT(IN)                                :: l, m
     341             :       INTEGER(KIND=int_8)                                :: nsum_2c_gspace_1d
     342             : 
     343     1352088 :       nsum_2c_gspace_1d = NINT(ns_G*(2*exp_w + (l + m + 1)*5), KIND=int_8)
     344     1352088 :    END FUNCTION
     345             : 
     346             : ! **************************************************************************************************
     347             : !> \brief Compute Ewald-like sum for 2-center ERIs in G space in 1 dimension
     348             : !>        S_G(l, alpha) = (-i)^l*inv_lgth*sum_G( C(l, alpha, G) exp(iGR) ), with
     349             : !>                        C(l, alpha, r) = r^l exp(-alpha*r^2),
     350             : !>        dG = inv_lgth*twopi and G = -G_bound*dG, (-G_bound + 1)*dG, ..., G_bound*dG
     351             : !>             for all l < = l_max.
     352             : !> \param S_G ...
     353             : !> \param R ...
     354             : !> \param alpha ...
     355             : !> \param inv_lgth ...
     356             : !> \param G_c ...
     357             : !> \note  S_G is real.
     358             : ! **************************************************************************************************
     359      123584 :    PURE SUBROUTINE pgf_sum_2c_gspace_1d(S_G, R, alpha, inv_lgth, G_c)
     360             :       REAL(KIND=dp), DIMENSION(0:), INTENT(OUT)          :: S_G
     361             :       REAL(KIND=dp), INTENT(IN)                          :: R, alpha, inv_lgth, G_c
     362             : 
     363             :       COMPLEX(KIND=dp)                                   :: exp_tot
     364      123584 :       COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:)        :: S_G_c
     365             :       INTEGER                                            :: gg, l, l_max
     366             :       REAL(KIND=dp)                                      :: dG, G, G_pow_l
     367             : 
     368      123584 :       dG = inv_lgth*twopi
     369      123584 :       l_max = UBOUND(S_G, 1)
     370             : 
     371      370752 :       ALLOCATE (S_G_c(0:l_max))
     372      673681 :       S_G_c(:) = 0.0_dp
     373      623238 :       DO gg = -FLOOR(G_c), FLOOR(G_c)
     374      499654 :          G = gg*dG
     375      499654 :          exp_tot = EXP(-alpha*G**2)*EXP(gaussi*G*R) ! cost: 2 exp_w flops
     376      499654 :          G_pow_l = 1.0_dp
     377     2701067 :          DO l = 0, l_max
     378     2077829 :             S_G_c(l) = S_G_c(l) + G_pow_l*(-1.0_dp)**l*exp_tot ! cost: 4 flops
     379     2577483 :             G_pow_l = G_pow_l*G ! cost: 1 flop
     380             :          ENDDO
     381             :       ENDDO
     382             : 
     383     1223778 :       S_G(:) = REAL(S_G_c(0:l_max)*i_pow((/(l, l=0, l_max)/)))*inv_lgth
     384      123584 :    END SUBROUTINE pgf_sum_2c_gspace_1d
     385             : 
     386             : ! **************************************************************************************************
     387             : !> \brief Roughly estimated number of floating point operations
     388             : !> \param ns_G ...
     389             : !> \param l ...
     390             : !> \param m ...
     391             : !> \return ...
     392             : ! **************************************************************************************************
     393       31248 :    PURE FUNCTION nsum_2c_gspace_3d(ns_G, l, m)
     394             :       REAL(KIND=dp), INTENT(IN)                          :: ns_G
     395             :       INTEGER, INTENT(IN)                                :: l, m
     396             :       INTEGER(KIND=int_8)                                :: nsum_2c_gspace_3d
     397             : 
     398       31248 :       nsum_2c_gspace_3d = NINT(ns_G*(2*exp_w + ncoset(l + m)*7), KIND=int_8)
     399       31248 :    END FUNCTION
     400             : 
     401             : ! **************************************************************************************************
     402             : !> \brief As pgf_sum_2c_gspace_1d but 3d sum required for non-orthorhombic cells
     403             : !> \param S_G ...
     404             : !> \param l_max ...
     405             : !> \param R ...
     406             : !> \param alpha ...
     407             : !> \param h_inv ...
     408             : !> \param G_c ...
     409             : !> \param G_rad ...
     410             : !> \param vol ...
     411             : !> \param coulomb ...
     412             : !> \note  MMME Method is not very efficient for non-orthorhombic cells
     413             : ! **************************************************************************************************
     414        5277 :    PURE SUBROUTINE pgf_sum_2c_gspace_3d(S_G, l_max, R, alpha, h_inv, G_c, G_rad, vol, potential, pot_par)
     415             :       REAL(KIND=dp), DIMENSION(:), INTENT(OUT)           :: S_G
     416             :       INTEGER, INTENT(IN)                                :: l_max
     417             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: R
     418             :       REAL(KIND=dp), INTENT(IN)                          :: alpha
     419             :       REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN)         :: h_inv
     420             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: G_c
     421             :       REAL(KIND=dp), INTENT(IN)                          :: G_rad, vol
     422             :       INTEGER, INTENT(IN), OPTIONAL                      :: potential
     423             :       REAL(KIND=dp), INTENT(IN), OPTIONAL                :: pot_par
     424             : 
     425             :       COMPLEX(KIND=dp)                                   :: exp_tot
     426             :       INTEGER, DIMENSION(3)                              :: l_xyz
     427             :       INTEGER                                            :: gx, gy, gz, k, l, lco, lx, ly, lz
     428       10554 :       COMPLEX(KIND=dp), DIMENSION(ncoset(l_max))         :: Ig
     429             :       REAL(KIND=dp)                                      :: G_rads_sq, G_sq
     430             :       REAL(KIND=dp), DIMENSION(3)                        :: G, G_x, G_y, G_z
     431        5277 :       REAL(KIND=dp), DIMENSION(3, 0:l_max)               :: G_pow_l
     432             :       REAL(KIND=dp), DIMENSION(3, 3)                     :: ht
     433             : 
     434       68601 :       ht = twopi*TRANSPOSE(h_inv)
     435      166900 :       Ig(:) = 0.0_dp
     436             : 
     437        5277 :       G_rads_sq = G_rad**2
     438             : 
     439       37652 :       DO gx = -FLOOR(G_c(1)), FLOOR(G_c(1))
     440      129500 :          G_x = ht(:, 1)*gx
     441      277609 :          DO gy = -FLOOR(G_c(2)), FLOOR(G_c(2))
     442      959828 :             G_y = ht(:, 2)*gy
     443     3495819 :             DO gz = -FLOOR(G_c(3)), FLOOR(G_c(3))
     444    12893948 :                G_z = ht(:, 3)*gz
     445             : 
     446    12893948 :                G = G_x + G_y + G_z
     447     3223487 :                G_sq = G(1)**2 + G(2)**2 + G(3)**2
     448     3223487 :                IF (G_sq .GT. G_rads_sq) CYCLE
     449             : 
     450     1738767 :                IF (PRESENT(potential)) THEN
     451     1033920 :                   IF (gx .EQ. 0 .AND. gy .EQ. 0 .AND. gz .EQ. 0) CYCLE
     452             :                ENDIF
     453             : 
     454     6954300 :                exp_tot = EXP(-alpha*G_sq)*EXP(gaussi*DOT_PRODUCT(G, R)) ! cost: 2 exp_w flops
     455     1738575 :                IF (PRESENT(potential)) THEN
     456     1378304 :                   SELECT CASE (potential)
     457             :                   CASE (eri_mme_coulomb)
     458      344576 :                      exp_tot = exp_tot/G_sq
     459             :                   CASE (eri_mme_yukawa)
     460      344576 :                      exp_tot = exp_tot/(G_sq + pot_par**2)
     461             :                      !exp_tot = exp_tot/G_sq
     462             :                   CASE (eri_mme_longrange)
     463     1033728 :                      exp_tot = exp_tot/G_sq*EXP(-G_sq/pot_par**2)
     464             :                   END SELECT
     465             :                ENDIF
     466     6954300 :                DO k = 1, 3
     467     5215725 :                   G_pow_l(k, 0) = 1.0_dp
     468    44295588 :                   DO l = 1, l_max
     469    42557013 :                      G_pow_l(k, l) = G_pow_l(k, l - 1)*G(k)
     470             :                   ENDDO
     471             :                ENDDO
     472   314858909 :                DO lco = 1, ncoset(l_max)
     473   312880377 :                   CALL get_l(lco, l, lx, ly, lz)
     474  1251521508 :                   l_xyz = [lx, ly, lz]
     475             :                   Ig(coset(lx, ly, lz)) = Ig(coset(lx, ly, lz)) + & ! cost: 7 flops
     476             :                                           G_pow_l(1, lx)*G_pow_l(2, ly)*G_pow_l(3, lz)* &
     477   316103864 :                                           exp_tot*(-1.0_dp)**l*i_pow(l)
     478             :                ENDDO
     479             :             ENDDO
     480             :          ENDDO
     481             :       ENDDO
     482      166900 :       S_G(:) = REAL(Ig(:), KIND=dp)/vol
     483        5277 :    END SUBROUTINE pgf_sum_2c_gspace_3d
     484             : 
     485             : ! **************************************************************************************************
     486             : !> \brief Roughly estimated number of floating point operations
     487             : !> \param ns_R ...
     488             : !> \param l ...
     489             : !> \param m ...
     490             : !> \return ...
     491             : ! **************************************************************************************************
     492     1352088 :    PURE FUNCTION nsum_2c_rspace_1d(ns_R, l, m)
     493             :       REAL(KIND=dp), INTENT(IN)                          :: ns_R
     494             :       INTEGER, INTENT(IN)                                :: l, m
     495             :       INTEGER(KIND=int_8)                                :: nsum_2c_rspace_1d
     496             : 
     497     1352088 :       nsum_2c_rspace_1d = NINT(ns_R*(exp_w + (l + m + 1)*3), KIND=int_8)
     498     1352088 :    END FUNCTION
     499             : 
     500             : ! **************************************************************************************************
     501             : !> \brief Compute Ewald-like sum for 2-center ERIs in R space in 1 dimension
     502             : !>        S_R(l, alpha) = SQRT(alpha/pi) sum_R'( H(l, alpha, R-R') ),
     503             : !>        with H(l, alpha, R) = (-d/dR)^l exp(-alpha*R^2),
     504             : !>        dR = lgth and R' = -R_min*dR, (-R_min + 1)*dR, ..., R_max*dR,
     505             : !>        for all l < = l_max.
     506             : !> \param S_R ...
     507             : !> \param R ...
     508             : !> \param alpha ...
     509             : !> \param lgth ...
     510             : !> \param R_c ...
     511             : !> \note  result is equivalent to pgf_sum_2c_gspace_1d with
     512             : !>              S_R(l, alpha) = S_G(l, 1/(4*alpha))
     513             : ! **************************************************************************************************
     514     1352614 :    PURE SUBROUTINE pgf_sum_2c_rspace_1d(S_R, R, alpha, lgth, R_c)
     515             :       REAL(KIND=dp), DIMENSION(0:), INTENT(OUT)          :: S_R
     516             :       REAL(KIND=dp), INTENT(IN)                          :: R, alpha, lgth, R_c
     517             : 
     518             :       INTEGER                                            :: l, l_max, rr
     519             :       REAL(KIND=dp)                                      :: dR, exp_tot, R_pow_l, Rp
     520     1352614 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: h_to_c
     521             : 
     522     1352614 :       dR = lgth
     523     2705228 :       l_max = UBOUND(S_R, 1)
     524             : 
     525             :       ! 1) compute sum over C(l, alpha, R - R') instead of H(l, alpha, R - R')
     526     5842547 :       S_R(:) = 0.0_dp
     527     1352614 :       Rp = R - R_c*dR
     528     5440994 :       DO rr = CEILING(-R_c - R/dR), FLOOR(R_c - R/dR)
     529     4088380 :          Rp = R + rr*dR
     530     4088380 :          exp_tot = EXP(-alpha*Rp**2) ! cost: exp_w flops
     531     4088380 :          R_pow_l = 1.0_dp
     532    19424628 :          DO l = 0, l_max
     533    13983634 :             S_R(l) = S_R(l) + R_pow_l*exp_tot ! cost: 2 flops
     534    18072014 :             R_pow_l = R_pow_l*Rp ! cost: 1 flop
     535             :          ENDDO
     536             :       ENDDO
     537             : 
     538             :       ! 2) C --> H
     539             :       CALL create_hermite_to_cartesian(alpha, l_max, h_to_c)
     540     5842547 :       S_R = MATMUL(TRANSPOSE(h_to_c(0:l_max, 0:l_max)), S_R)*SQRT(alpha/pi)
     541     2705228 :    END SUBROUTINE pgf_sum_2c_rspace_1d
     542             : 
     543             : ! **************************************************************************************************
     544             : !> \brief Roughly estimated number of floating point operations
     545             : !> \param ns_R ...
     546             : !> \param l ...
     547             : !> \param m ...
     548             : !> \return ...
     549             : ! **************************************************************************************************
     550       31248 :    PURE FUNCTION nsum_2c_rspace_3d(ns_R, l, m)
     551             :       REAL(KIND=dp), INTENT(IN)                          :: ns_R
     552             :       INTEGER, INTENT(IN)                                :: l, m
     553             :       INTEGER(KIND=int_8)                                :: nsum_2c_rspace_3d
     554             : 
     555       31248 :       nsum_2c_rspace_3d = NINT(ns_R*(exp_w + ncoset(l + m)*(4 + ncoset(l + m)*4)), KIND=int_8)
     556       31248 :    END FUNCTION
     557             : 
     558             : ! **************************************************************************************************
     559             : !> \brief As pgf_sum_2c_rspace_1d but 3d sum required for non-orthorhombic cells
     560             : !> \param S_R ...
     561             : !> \param l_max ...
     562             : !> \param R ...
     563             : !> \param alpha ...
     564             : !> \param hmat ...
     565             : !> \param h_inv ...
     566             : !> \param R_c ...
     567             : !> \param R_rad ...
     568             : !> \note  MMME Method is not very efficient for non-orthorhombic cells
     569             : ! **************************************************************************************************
     570       35163 :    PURE SUBROUTINE pgf_sum_2c_rspace_3d(S_R, l_max, R, alpha, hmat, h_inv, R_c, R_rad)
     571             :       REAL(KIND=dp), DIMENSION(:), INTENT(OUT)           :: S_R
     572             :       INTEGER, INTENT(IN)                                :: l_max
     573             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: R
     574             :       REAL(KIND=dp), INTENT(IN)                          :: alpha
     575             :       REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN)         :: hmat, h_inv
     576             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: R_c
     577             :       REAL(KIND=dp), INTENT(IN)                          :: R_rad
     578             : 
     579             :       INTEGER                                            :: k, l, lco, llx, lly, llz, lx, ly, lz, &
     580             :                                                             sx, sy, sz
     581             :       REAL(KIND=dp)                                      :: exp_tot, R_rad_sq, R_sq
     582       35163 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: h_to_c
     583             :       REAL(KIND=dp), DIMENSION(3)                        :: R_l, R_r, Rp, Rx, Ry, Rz, s_shift
     584       70326 :       REAL(KIND=dp), DIMENSION(3, 0:l_max)               :: R_pow_l
     585       35163 :       REAL(KIND=dp), DIMENSION(ncoset(l_max))            :: S_R_C
     586             : 
     587      445064 :       S_R(:) = 0.0_dp
     588      445064 :       S_R_C(:) = 0.0_dp
     589             : 
     590      562608 :       s_shift = MATMUL(h_inv, -R)
     591      140652 :       R_l = -R_c + s_shift
     592      140652 :       R_r = R_c + s_shift
     593             : 
     594       35163 :       R_rad_sq = R_rad**2
     595             : 
     596       95444 :       DO sx = CEILING(R_l(1)), FLOOR(R_r(1))
     597      241124 :          Rx = hmat(:, 1)*sx
     598      251589 :          DO sy = CEILING(R_l(2)), FLOOR(R_r(2))
     599      624580 :             Ry = hmat(:, 2)*sy
     600      782873 :             DO sz = CEILING(R_l(3)), FLOOR(R_r(3))
     601     2265788 :                Rz = hmat(:, 3)*sz
     602     2265788 :                Rp = Rx + Ry + Rz
     603      566447 :                R_sq = (Rp(1) + R(1))**2 + (Rp(2) + R(2))**2 + (Rp(3) + R(3))**2
     604      566447 :                IF (R_sq .GT. R_rad_sq) CYCLE
     605      296776 :                exp_tot = EXP(-alpha*R_sq) ! cost: exp_w flops
     606     1187104 :                DO k = 1, 3
     607      890328 :                   R_pow_l(k, 0) = 1.0_dp
     608     2374189 :                   DO l = 1, l_max
     609     2077413 :                      R_pow_l(k, l) = R_pow_l(k, l - 1)*(Rp(k) + R(k))
     610             :                   ENDDO
     611             :                ENDDO
     612     2942472 :                DO lco = 1, ncoset(l_max)
     613     2489551 :                   CALL get_l(lco, l, lx, ly, lz)
     614     3055998 :                   S_R_C(coset(lx, ly, lz)) = S_R_C(coset(lx, ly, lz)) + R_pow_l(1, lx)*R_pow_l(2, ly)*R_pow_l(3, lz)*exp_tot ! cost: 4 flops
     615             :                ENDDO
     616             :             ENDDO
     617             :          ENDDO
     618             :       ENDDO
     619             : 
     620             :       CALL create_hermite_to_cartesian(alpha, l_max, h_to_c)
     621             : 
     622      445064 :       DO lco = 1, ncoset(l_max)
     623      409901 :          CALL get_l(lco, l, lx, ly, lz)
     624             : 
     625     1195340 :          DO llx = 0, lx
     626     2453751 :          DO lly = 0, ly
     627     4170544 :          DO llz = 0, lz
     628             :             S_R(coset(lx, ly, lz)) = S_R(coset(lx, ly, lz)) + & ! cost: 4 flops
     629             :                                      h_to_c(llx, lx)*h_to_c(lly, ly)*h_to_c(llz, lz)* &
     630     3420268 :                                      S_R_C(coset(llx, lly, llz))
     631             :          ENDDO
     632             :          ENDDO
     633             :          ENDDO
     634             :       ENDDO
     635      445064 :       S_R(:) = S_R(:)*(alpha/pi)**1.5_dp
     636             : 
     637       70326 :    END SUBROUTINE pgf_sum_2c_rspace_3d
     638             : 
     639             : ! **************************************************************************************************
     640             : !> \brief Compute 1d sum
     641             : !>        S_G(l, alpha) = inv_lgth*sum_G( C(l, alpha, delta_l, G) ) with
     642             : !>          C(l, alpha, delta_l, G) = prefactor*|G|^(l-delta_l) exp(-alpha*G^2)
     643             : !>          if G not equal 0
     644             : !>          C(l = 0, alpha, delta_l, 0) = 1, C(l>0, alpha, delta_l, 0) = 0
     645             : !>        dG = inv_lgth*twopi and G = -G_bound*dG, (-G_bound + 1)*dG, ..., G_bound*dG
     646             : !>        for all l < = l_max.
     647             : !> \param S_G ...
     648             : !> \param alpha ...
     649             : !> \param inv_lgth ...
     650             : !> \param G_min ...
     651             : !> \param G_c ...
     652             : !> \param delta_l ...
     653             : !> \param prefactor ...
     654             : !> \note  needed for cutoff error estimate
     655             : ! **************************************************************************************************
     656      262644 :    PURE SUBROUTINE pgf_sum_2c_gspace_1d_deltal(S_G, alpha, inv_lgth, G_min, G_c, delta_l, prefactor)
     657             :       REAL(KIND=dp), DIMENSION(0:), INTENT(OUT)          :: S_G
     658             :       REAL(KIND=dp), INTENT(IN)                          :: alpha, inv_lgth
     659             :       INTEGER, INTENT(IN)                                :: G_min, G_c
     660             :       REAL(KIND=dp), INTENT(IN)                          :: delta_l, prefactor
     661             : 
     662             :       INTEGER                                            :: k, k0, k1, l, l_max
     663             :       REAL(KIND=dp)                                      :: dG, exp_tot, G, prefac
     664             : 
     665      262644 :       prefac = prefactor*inv_lgth
     666      262644 :       dG = inv_lgth*twopi ! positive
     667      525288 :       l_max = UBOUND(S_G, 1)
     668             : 
     669     1490184 :       S_G(:) = 0.0_dp
     670      262644 :       IF (0 .LT. G_min) THEN
     671             :          k0 = G_min; k1 = 0
     672       65661 :       ELSE IF (G_c .LT. 0) THEN
     673             :          k0 = 0; k1 = G_c
     674             :       ELSE ! Gmin <= 0 /\ 0 <= Gc
     675       65661 :          S_G(0) = prefac
     676             :          k0 = 1; k1 = -1
     677             :       ENDIF
     678             :       ! positive range
     679             :       IF (0 .LT. k0) THEN
     680    13363636 :          DO k = k0, G_c
     681    13100992 :             G = k*dG; exp_tot = EXP(-alpha*G**2)*prefac
     682    78655736 :             DO l = 0, l_max
     683    78393092 :                S_G(l) = S_G(l) + G**(l - delta_l)*exp_tot
     684             :             ENDDO
     685             :          ENDDO
     686             :       ENDIF
     687             :       ! negative range
     688      262644 :       IF (k1 .LT. 0) THEN
     689     4437282 :          DO k = G_min, k1
     690     4371621 :             G = k*dG; exp_tot = EXP(-alpha*G**2)*prefac
     691    25755453 :             DO l = 0, l_max
     692    25689792 :                S_G(l) = S_G(l) + ABS(G)**(l - delta_l)*exp_tot
     693             :             ENDDO
     694             :          ENDDO
     695             :       ENDIF
     696      262644 :    END SUBROUTINE pgf_sum_2c_gspace_1d_deltal
     697             : 
     698             : ! **************************************************************************************************
     699             : !> \brief Compute Ewald-like sum for 3-center integrals in 1 dimension
     700             : !>        S_G(l, m, n, alpha, beta, gamma) = i^(l+m+n)*(-1)^(l+m)*inv_lgth^2*
     701             : !>                                           sum_G sum_G'( exp(i G R1)
     702             : !>                                           C(l,alpha,G) C(m,beta,G'-G) C(n,gamma,G') exp(i G' R2) )
     703             : !>        for all l < = l_max, m <= m_max, n <= n_max.
     704             : !>        a_mm is the minimax exponent.
     705             : !>        alpha =  1/(4 zeta), beta = 1/(4 zetb), gamma = 1/(4 zetc) + a_mm
     706             : !>        R1 = RB-RA; R2 = RC-RB
     707             : !>        Note on method / order arguments:
     708             : !>        Three equivalent methods (Poisson summation) to compute this sum over
     709             : !>        Cartesian Gaussians C or Hermite Gaussians H and
     710             : !>        reciprocal lattice vectors G or direct lattice vectors R:
     711             : !>        - method 1: sum_G sum_G' C(G) C(G,G') C(G')
     712             : !>        - method 2: sum_G sum_R C(G) C(R)
     713             : !>        - method 3: sum_R sum_R' H(R, R')
     714             : !>        The order parameter selects the Gaussian functions over which the sum is performed
     715             : !>        method 1: order = 1, 2, 3
     716             : !>        method 2: order = 1, 2, 3
     717             : !>        method 3: order = 1
     718             : !>        If method and order are not present, the method / order that converges fastest is
     719             : !>        automatically chosen.
     720             : !> \param S_G ...
     721             : !> \param RA ...
     722             : !> \param RB ...
     723             : !> \param RC ...
     724             : !> \param zeta ...
     725             : !> \param zetb ...
     726             : !> \param zetc ...
     727             : !> \param a_mm ...
     728             : !> \param lgth ...
     729             : !> \param G_bounds_1 ...
     730             : !> \param R_bounds_2 ...
     731             : !> \param R_bounds_3 ...
     732             : !> \param method ...
     733             : !> \param method_out ...
     734             : !> \param order ...
     735             : ! **************************************************************************************************
     736     3796920 :    SUBROUTINE pgf_sum_3c_1d(S_G, RA, RB, RC, zeta, zetb, zetc, a_mm, lgth, R_bounds_3)
     737             :       REAL(KIND=dp), DIMENSION(0:, 0:, 0:), INTENT(OUT)  :: S_G
     738             :       REAL(KIND=dp), INTENT(IN)                          :: RA, RB, RC, zeta, zetb, zetc, a_mm, lgth
     739             :       REAL(KIND=dp), DIMENSION(2), INTENT(IN)            :: R_bounds_3
     740             : 
     741             :       INTEGER                                            :: l_max, m_max, n_max
     742             :       INTEGER                                            :: nR1, nR2
     743             :       INTEGER                                            :: prop_exp
     744             : 
     745     3796920 :       l_max = UBOUND(S_G, 1)
     746     3796920 :       m_max = UBOUND(S_G, 2)
     747     3796920 :       n_max = UBOUND(S_G, 3)
     748             : 
     749     3796920 :       nR1 = 2*FLOOR(R_bounds_3(1)) + 1
     750     3796920 :       nR2 = 2*FLOOR(R_bounds_3(2)) + 1
     751             : 
     752     3796920 :       IF (nR1*nR2 > 1 + nR1*2) THEN
     753             :          prop_exp = 1
     754             :       ELSE
     755     2261286 :          prop_exp = 0
     756             :       ENDIF
     757             : 
     758    15187680 :       IF (MAXVAL([l_max, m_max, n_max]) .GT. ${lmax_unroll}$) THEN
     759       92160 :          CALL pgf_sum_3c_rspace_1d_generic(S_G, RA, RB, RC, zeta, zetb, zetc, a_mm, lgth, R_bounds_3)
     760             :       ELSE
     761             : #:for l_max in range(0,lmax_unroll+1)
     762     3704760 :          IF (l_max == ${l_max}$) THEN
     763             : #:for m_max in range(0,lmax_unroll+1)
     764     3704760 :             IF (m_max == ${m_max}$) THEN
     765             : #:for n_max in range(0, lmax_unroll+1)
     766     3704760 :                IF (n_max == ${n_max}$) THEN
     767             : #:for prop_exp in range(0,2)
     768     3704760 :                   IF (prop_exp == ${prop_exp}$) THEN
     769             :                      CALL pgf_sum_3c_rspace_1d_${l_max}$_${m_max}$_${n_max}$_exp_${prop_exp}$ (S_G, RA, RB, RC, &
     770     3704760 :                                                                                            zeta, zetb, zetc, a_mm, lgth, R_bounds_3)
     771     3704760 :                      RETURN
     772             :                   ENDIF
     773             : #:endfor
     774             :                ENDIF
     775             : #:endfor
     776             :             ENDIF
     777             : #:endfor
     778             :          ENDIF
     779             : #:endfor
     780             :       ENDIF
     781             : 
     782             :    END SUBROUTINE pgf_sum_3c_1d
     783             : 
     784             : ! **************************************************************************************************
     785             : !> \brief Roughly estimated number of floating point operations
     786             : !> \param ns_G1 ...
     787             : !> \param ns_G2 ...
     788             : !> \param l ...
     789             : !> \param m ...
     790             : !> \param n ...
     791             : !> \return ...
     792             : ! **************************************************************************************************
     793           0 :    PURE FUNCTION nsum_3c_gspace_1d()
     794             :       INTEGER(KIND=int_8)                                :: nsum_3c_gspace_1d
     795             : 
     796           0 :       nsum_3c_gspace_1d = 15
     797           0 :    END FUNCTION
     798             : 
     799             : ! **************************************************************************************************
     800             : !> \brief Roughly estimated number of floating point operations
     801             : !> \param ns_G ...
     802             : !> \param ns_R ...
     803             : !> \param l ...
     804             : !> \param m ...
     805             : !> \param n ...
     806             : !> \return ...
     807             : ! **************************************************************************************************
     808           0 :    PURE FUNCTION nsum_product_3c_gspace_1d(ns_G, ns_R)
     809             :       REAL(KIND=dp), INTENT(IN) :: ns_G, ns_R
     810             :       INTEGER(KIND=int_8)                                :: nsum_product_3c_gspace_1d
     811             : 
     812           0 :       nsum_product_3c_gspace_1d = MIN(19, NINT(ns_G*(3 + ns_R*2)))
     813           0 :    END FUNCTION
     814             : 
     815             : ! **************************************************************************************************
     816             : !> \brief Roughly estimated number of floating point operations
     817             : !> \param ns_R1 ...
     818             : !> \param ns_R2 ...
     819             : !> \param l ...
     820             : !> \param m ...
     821             : !> \param n ...
     822             : !> \return ...
     823             : ! **************************************************************************************************
     824     3796920 :    PURE FUNCTION nsum_3c_rspace_1d(ns_R1, ns_R2)
     825             :       REAL(KIND=dp), INTENT(IN)                          :: ns_R1, ns_R2
     826             :       INTEGER(KIND=int_8)                                :: nsum_3c_rspace_1d
     827             : 
     828     3796920 :       nsum_3c_rspace_1d = NINT(MIN((4 + ns_R1*2), ns_R1*(ns_R2 + 1)), KIND=int_8)
     829     3796920 :    END FUNCTION
     830             : 
     831             : ! **************************************************************************************************
     832             : !> \brief Helper routine: compute SQRT(alpha/pi) (-1)^n sum_(R, R') sum_{t=0}^{l+m} E(t,l,m) H(RC - P(R) - R', t + n, alpha)
     833             : !> with alpha = 1.0_dp/((a + b + c)/((a + b)*c) + 4.0_dp*a_mm),
     834             : !> P(R) = (a*(RA + R) + b*RB)/(a + b)
     835             : !> \param S_R ...
     836             : !> \param RA ...
     837             : !> \param RB ...
     838             : !> \param RC ...
     839             : !> \param zeta ...
     840             : !> \param zetb ...
     841             : !> \param zetc ...
     842             : !> \param a_mm ...
     843             : !> \param lgth ...
     844             : !> \param R_c ...
     845             : ! **************************************************************************************************
     846       92160 :    PURE SUBROUTINE pgf_sum_3c_rspace_1d_generic(S_R, RA, RB, RC, zeta, zetb, zetc, a_mm, lgth, R_c)
     847             :       REAL(KIND=dp), DIMENSION(0:, 0:, 0:), INTENT(OUT)  :: S_R
     848             :       REAL(KIND=dp), INTENT(IN)                          :: RA, RB, RC, zeta, zetb, zetc, a_mm, lgth
     849             :       REAL(KIND=dp), DIMENSION(2), INTENT(IN)            :: R_c
     850             : 
     851             :       INTEGER                                            :: ll, mm, l, k, l_max, m, m_max, n, n_max, rr1, rr2, t, l_tot_max
     852             :       REAL(KIND=dp)                                      :: alpha, dR, exp_tot, R, R1, R2, R_offset, &
     853             :                                                             R_pow_t, R_tmp, rr1_delta, rr2_delta, c1, c2, c3
     854             :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: S_R_t
     855       92160 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: h_to_c
     856       92160 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: E
     857             : 
     858       92160 :       dR = lgth
     859       92160 :       alpha = 1.0_dp/((zeta + zetb + zetc)/((zeta + zetb)*zetc) + 4.0_dp*a_mm)
     860       92160 :       l_max = UBOUND(S_R, 1)
     861       92160 :       m_max = UBOUND(S_R, 2)
     862       92160 :       n_max = UBOUND(S_R, 3)
     863       92160 :       l_tot_max = l_max + m_max + n_max
     864             : 
     865      276480 :       ALLOCATE (S_R_t(0:l_max + m_max + n_max))
     866      460800 :       ALLOCATE (E(-1:l_max + m_max + 1, -1:l_max, -1:m_max))
     867             : 
     868     2903040 :       S_R(:, :, :) = 0.0_dp
     869             : 
     870       92160 :       R_offset = RC - (zeta*RA + zetb*RB)/(zeta + zetb)
     871             : 
     872             :       ! inline CALL create_hermite_to_cartesian(alpha, l_tot_max, h_to_c)
     873      368640 :       ALLOCATE (h_to_c(-1:l_tot_max + 1, 0:l_tot_max))
     874     8202240 :       h_to_c(:, :) = 0.0_dp
     875       92160 :       h_to_c(0, 0) = 1.0_dp
     876      737280 :       DO l = 0, l_tot_max - 1
     877     3962880 :          DO k = 0, l + 1
     878     3870720 :             h_to_c(k, l + 1) = -(k + 1)*h_to_c(k + 1, l) + 2.0_dp*alpha*h_to_c(k - 1, l)
     879             :          ENDDO
     880             :       ENDDO
     881             : 
     882       92160 :       rr1_delta = (RA - RB)/dR
     883     1428480 :       DO rr1 = CEILING(-R_c(1) + rr1_delta), FLOOR(R_c(1) + rr1_delta)
     884    12026880 :          S_R_t(:) = 0.0_dp
     885     1336320 :          R1 = rr1*dR
     886     1336320 :          R_tmp = R_offset + R1*zeta/(zeta + zetb)
     887     1336320 :          rr2_delta = -R_tmp/dR
     888    15398400 :          DO rr2 = CEILING(-R_c(2) + rr2_delta), FLOOR(R_c(2) + rr2_delta)
     889    14062080 :             R2 = rr2*dR
     890    14062080 :             R = R_tmp + R2
     891    14062080 :             exp_tot = EXP(-alpha*R**2) ! cost: exp_w flops
     892    14062080 :             R_pow_t = 1.0_dp
     893   127895040 :             DO t = 0, l_max + m_max + n_max
     894   112496640 :                S_R_t(t) = S_R_t(t) + R_pow_t*exp_tot ! cost: 2 flops
     895   126558720 :                R_pow_t = R_pow_t*R ! cost: 1 flop
     896             :             ENDDO
     897             :          ENDDO
     898             : 
     899             :          ! C --> H
     900    12026880 :          S_R_t(:) = MATMUL(TRANSPOSE(h_to_c(0:l_max + m_max + n_max, 0:l_max + m_max + n_max)), S_R_t)*SQRT(alpha/pi)
     901             : 
     902             :          ! H --> HH
     903             :          !inline CALL create_gaussian_overlap_dist_to_hermite(l_max, m_max, zeta, zetb, RA-R1, RB, 2, E)
     904             : 
     905   183191040 :          E(:, :, :) = 0.0_dp
     906     1336320 :          E(0, 0, 0) = EXP(-zeta*zetb/(zeta + zetb)*(RA - R1 - RB)**2)
     907             : 
     908     1336320 :          c1 = 1.0_dp/(zeta + zetb)
     909     1336320 :          c2 = 2.0_dp*(zetb/(zeta + zetb))*(RB - (RA - R1))
     910     1336320 :          c3 = 2.0_dp*(zeta/(zeta + zetb))*(RA - R1 - RB)
     911             : 
     912     7948800 :          DO mm = 0, m_max - 1
     913     8179200 :             DO ll = 0, l_max - 1
     914     7879680 :                DO t = 0, ll + mm + 1
     915             :                   E(t, ll + 1, mm) = zeta*(c1*E(t - 1, ll, mm) + &
     916             :                                            c2*E(t, ll, mm) + &
     917             :                                            2*(t + 1)*E(t + 1, ll, mm) - &
     918     1036800 :                                            2*ll*E(t, ll - 1, mm))
     919             :                   E(t, ll, mm + 1) = zetb*(c1*E(t - 1, ll, mm) + &
     920             :                                            c3*E(t, ll, mm) + &
     921             :                                            2*(t + 1)*E(t + 1, ll, mm) - &
     922     1267200 :                                            2*mm*E(t, ll, mm - 1))
     923             :                ENDDO
     924             :             ENDDO
     925             :          ENDDO
     926             : 
     927     1451520 :          DO ll = 0, l_max - 1
     928     2142720 :             DO t = 0, ll + m_max + 1
     929             :                E(t, ll + 1, m_max) = zeta*(c1*E(t - 1, ll, m_max) + &
     930             :                                            c2*E(t, ll, m_max) + &
     931             :                                            2*(t + 1)*E(t + 1, ll, m_max) - &
     932      806400 :                                            2*ll*E(t, ll - 1, m_max))
     933             :             ENDDO
     934             :          ENDDO
     935             : 
     936     7948800 :          DO mm = 0, m_max - 1
     937    34560000 :             DO t = 0, l_max + mm + 1
     938             :                E(t, l_max, mm + 1) = zetb*(c1*E(t - 1, l_max, mm) + &
     939             :                                            c3*E(t, l_max, mm) + &
     940             :                                            2*(t + 1)*E(t + 1, l_max, mm) - &
     941    33223680 :                                            2*mm*E(t, l_max, mm - 1))
     942             :             ENDDO
     943             :          ENDDO
     944             : 
     945     5391360 :          DO n = 0, n_max
     946    29007360 :             DO m = 0, m_max
     947    51724800 :                DO l = 0, l_max
     948   132364800 :                   DO t = 0, l + m
     949   108656640 :                      S_R(l, m, n) = S_R(l, m, n) + E(t, l, m)*(-1)**n*S_R_t(t + n) ! cost: 5 flops
     950             :                   ENDDO
     951             :                ENDDO
     952             :             ENDDO
     953             :          ENDDO
     954             :       ENDDO
     955             : 
     956     2903040 :       S_R = S_R*pi**(-0.5_dp)*((zeta + zetb)/(zeta*zetb))**(-0.5_dp)
     957       92160 :    END SUBROUTINE
     958             : 
     959             : ! **************************************************************************************************
     960             : !> \brief Helper routine: compute SQRT(alpha/pi) (-1)^n sum_(R, R') sum_{t=0}^{l+m} E(t,l,m) H(RC - P(R) - R', t + n, alpha)
     961             : !> with alpha = 1.0_dp/((a + b + c)/((a + b)*c) + 4.0_dp*a_mm),
     962             : !> P(R) = (a*(RA + R) + b*RB)/(a + b)
     963             : !> \param S_R ...
     964             : !> \param RA ...
     965             : !> \param RB ...
     966             : !> \param RC ...
     967             : !> \param zeta ...
     968             : !> \param zetb ...
     969             : !> \param zetc ...
     970             : !> \param a_mm ...
     971             : !> \param lgth ...
     972             : !> \param R_c ...
     973             : ! **************************************************************************************************
     974             : #:for prop_exp in range(0,2)
     975             : #:for l_max in range(0, lmax_unroll+1)
     976             : #:for m_max in range(0, lmax_unroll+1)
     977             : #:for n_max in range(0, lmax_unroll+1)
     978             : #:set l_tot_max = l_max + m_max + n_max
     979     3704760 :    PURE SUBROUTINE pgf_sum_3c_rspace_1d_${l_max}$_${m_max}$_${n_max}$_exp_${prop_exp}$ ( &
     980     3704760 :       S_R, RA, RB, RC, zeta, zetb, zetc, a_mm, lgth, R_c)
     981             :       REAL(KIND=dp), DIMENSION(0:, 0:, 0:), INTENT(OUT)  :: S_R
     982             :       REAL(KIND=dp), INTENT(IN)                          :: RA, RB, RC, zeta, zetb, zetc, a_mm, lgth
     983             :       REAL(KIND=dp), DIMENSION(2), INTENT(IN)            :: R_c
     984             :       INTEGER                                            :: rr1, rr2, rr2_l, rr2_r, rr1_l, rr1_r
     985             :       REAL(KIND=dp)                                      :: alpha, alpha_E, dR, exp2_Rsq, R, R1, R_offset, &
     986             :                                                             R_pow_t, R_tmp, rr1_delta, rr2_delta
     987             : 
     988             : #:if l_tot_max > 0
     989             :       REAL(KIND=dp)                                      :: c1, c2, c3
     990             : #:endif
     991             :       REAL(KIND=dp) :: ${", ".join(["S_R_t_{}".format(t) for t in range(0,l_tot_max+1)])}$
     992             :       REAL(KIND=dp) :: ${", ".join(["S_R_t2_{}".format(t) for t in range(0,l_tot_max+1)])}$
     993             :   REAL(KIND=dp) :: ${", ".join([", ".join(["h_to_c_{}_{}".format(l1,l2) for l1 in range(0,l2+1)]) for l2 in range(0,l_tot_max+1)])}$
     994             :       REAL(KIND=dp) :: ${", ".join([", ".join([", ".join(["E_{}_{}_{}".format(t,l,m) for t in range(0,l+m+1)]) for l in range(0,l_max+1)]) for m in range(0,m_max+1)])}$
     995             : 
     996             : #:if prop_exp
     997             :       REAL(KIND=dp) :: exp2_2RdR, exp_dRsq, exp_2dRsq !exp_E_dRsq, exp_E_2dRsq, exp_E_2RdR, exp_E_Rsq
     998             : #:endif
     999             : 
    1000     3704760 :       dR = lgth
    1001     3704760 :       alpha = 1.0_dp/((zeta + zetb + zetc)/((zeta + zetb)*zetc) + 4.0_dp*a_mm)
    1002             : 
    1003    71708472 :       S_R(:, :, :) = 0.0_dp
    1004             : 
    1005     3704760 :       R_offset = RC - (zeta*RA + zetb*RB)/(zeta + zetb)
    1006             : 
    1007     3704760 :       h_to_c_0_0 = SQRT(alpha/pi)
    1008             : 
    1009             : #:for l in range(0, l_tot_max)
    1010             : #:for k in range(0, l+2)
    1011             : #:if k<l or k>0
    1012     3345336 :       h_to_c_${k}$_${l+1}$ = #{if k<l}#-${k+1}$*h_to_c_${k+1}$_${l}$#{endif}# #{if k > 0}#+2*alpha*h_to_c_${k-1}$_${l}$#{endif}#
    1013             : #:else
    1014     3345336 :       h_to_c_${k}$_${l+1}$ = 0.0_dp
    1015             : #:endif
    1016             : #:endfor
    1017             : #:endfor
    1018             : 
    1019             : #:if prop_exp
    1020     1477266 :       exp_dRsq = exp(-alpha*dR*dR)
    1021     1477266 :       exp_2dRsq = exp_dRsq*exp_dRsq
    1022             : #:endif
    1023             : 
    1024     3704760 :       rr1_delta = (RA - RB)/dR
    1025             : 
    1026     3704760 :       rr1_l = CEILING(-R_c(1) + rr1_delta)
    1027     3704760 :       rr1_r = FLOOR(R_c(1) + rr1_delta)
    1028             : 
    1029     3704760 :       R1 = rr1_l*dR
    1030             : 
    1031     3704760 :       alpha_E = zeta*zetb/(zeta + zetb)
    1032             : 
    1033    11061672 :       DO rr1 = rr1_l, rr1_r
    1034             : #:for t in range(0, l_tot_max + 1)
    1035     7356912 :          S_R_t_${t}$ = 0.0_dp
    1036     7356912 :          S_R_t2_${t}$ = 0.0_dp
    1037             : #:endfor
    1038     7356912 :          R_tmp = R_offset + R1*zeta/(zeta + zetb)
    1039     7356912 :          rr2_delta = -R_tmp/dR
    1040             : 
    1041     7356912 :          rr2_l = CEILING(-R_c(2) + rr2_delta)
    1042     7356912 :          rr2_r = FLOOR(R_c(2) + rr2_delta)
    1043             : 
    1044     7356912 :          R = R_tmp + (rr2_l)*dR
    1045             : 
    1046             : #:if prop_exp
    1047     3805288 :          exp2_2RdR = exp(-2*alpha*R*dR)
    1048     3805288 :          exp2_Rsq = exp(-alpha*R*R)
    1049             : #:endif
    1050             : 
    1051    35134181 :          DO rr2 = rr2_l, rr2_r
    1052    27777269 :             R_pow_t = 1.0_dp
    1053             : #:if not prop_exp
    1054     4982766 :             exp2_Rsq = exp(-alpha*R*R)
    1055             : #:endif
    1056             : #:for t in range(0, l_tot_max + 1)
    1057    27777269 :             S_R_t_${t}$ = S_R_t_${t}$+R_pow_t*exp2_Rsq
    1058             : #:if t<l_tot_max
    1059    25662401 :             R_pow_t = R_pow_t*R
    1060             : #:endif
    1061             : #:endfor
    1062             : 
    1063             : #:if prop_exp
    1064    22794503 :             exp2_Rsq = exp2_Rsq*exp_dRsq*exp2_2RdR
    1065    22794503 :             exp2_2RdR = exp2_2RdR*exp_2dRsq
    1066             : #:endif
    1067    35134181 :             R = R + dR
    1068             :          ENDDO
    1069             : 
    1070             :          ! C --> H
    1071             : #:for l in range(0, l_tot_max+1)
    1072             : #:for k in range(0, l+1)
    1073     7356912 :          S_R_t2_${l}$ = S_R_t2_${l}$+h_to_c_${k}$_${l}$*S_R_t_${k}$
    1074             : #:endfor
    1075             : #:endfor
    1076             : 
    1077             :          ! H --> HH
    1078     7356912 :          E_0_0_0 = exp(-alpha_E*(RA - RB - R1)*(RA - RB - R1))
    1079             : 
    1080             : #:if l_tot_max > 0
    1081     6688752 :          c1 = 1.0_dp/(zeta + zetb)
    1082     6688752 :          c2 = 2.0_dp*(zetb/(zeta + zetb))*(RB - (RA - R1))
    1083     6688752 :          c3 = 2.0_dp*(zeta/(zeta + zetb))*(RA - R1 - RB)
    1084             : #:endif
    1085             : 
    1086             : #:for m in range(0,m_max+1)
    1087             : #:for l in range(0,l_max+1)
    1088             : #:for t in range(0,l+m+2)
    1089             : #:if l < l_max
    1090             :          E_${t}$_${l+1}$_${m}$ = zeta*(#{if t>0}# c1*E_${t-1}$_${l}$_${m}$#{endif}# &
    1091             :          #{if t<=l+m}# +c2*E_${t}$_${l}$_${m}$&#{endif}#
    1092             :          #{if t<l+m}# +${2*(t+1)}$*E_${t+1}$_${l}$_${m}$ &#{endif}#
    1093     4914576 :          #{if l>0 and t<=l-1+m}#-${2*l}$*E_${t}$_${l-1}$_${m}$#{endif}#)
    1094             : #:endif
    1095             : #:if m < m_max
    1096             :          E_${t}$_${l}$_${m+1}$ = zetb*(#{if t>0}# c1*E_${t-1}$_${l}$_${m}$#{endif}# &
    1097             :          #{if t<=l+m}#+c3*E_${t}$_${l}$_${m}$&#{endif}#
    1098             :          #{if t<l+m}# +${2*(t+1)}$*E_${t+1}$_${l}$_${m}$ &#{endif}#
    1099     4914576 :          #{if m>0 and t<=m-1+l}#-${2*m}$*E_${t}$_${l}$_${m-1}$#{endif}#)
    1100             : #:endif
    1101             : #:endfor
    1102             : #:endfor
    1103             : #:endfor
    1104             : 
    1105             : #:for n in range(0, n_max+1)
    1106             : #:for m in range(0, m_max+1)
    1107             : #:for l in range(0, l_max+1)
    1108             : #:for t in range(0, l+m+1)
    1109     7356912 :          S_R(${l}$, ${m}$, ${n}$) = S_R(${l}$, ${m}$, ${n}$) + E_${t}$_${l}$_${m}$*(${(-1)**n}$)*S_R_t2_${t+n}$ ! cost: 5 flops
    1110             : #:endfor
    1111             : #:endfor
    1112             : #:endfor
    1113             : #:endfor
    1114    11061672 :          R1 = R1 + dR
    1115             :       ENDDO
    1116             : 
    1117    71708472 :       S_R = S_R*pi**(-0.5_dp)*((zeta + zetb)/(zeta*zetb))**(-0.5_dp)
    1118     3704760 :    END SUBROUTINE
    1119             : #:endfor
    1120             : #:endfor
    1121             : #:endfor
    1122             : #:endfor
    1123             : 
    1124             : ! **************************************************************************************************
    1125             : !> \brief As pgf_sum_3c_1d but 3d sum required for non-orthorhombic cells
    1126             : !> \param S_G ...
    1127             : !> \param la_max ...
    1128             : !> \param lb_max ...
    1129             : !> \param lc_max ...
    1130             : !> \param RA ...
    1131             : !> \param RB ...
    1132             : !> \param RC ...
    1133             : !> \param zeta ...
    1134             : !> \param zetb ...
    1135             : !> \param zetc ...
    1136             : !> \param a_mm ...
    1137             : !> \param hmat ...
    1138             : !> \param h_inv ...
    1139             : !> \param vol ...
    1140             : !> \param G_bounds_1 ...
    1141             : !> \param R_bounds_2 ...
    1142             : !> \param R_bounds_3 ...
    1143             : !> \param G_rads_1 ...
    1144             : !> \param R_rads_2 ...
    1145             : !> \param R_rads_3 ...
    1146             : !> \param method ...
    1147             : !> \param method_out ...
    1148             : !> \param order ...
    1149             : ! **************************************************************************************************
    1150      148500 :    SUBROUTINE pgf_sum_3c_3d(S_G, la_max, lb_max, lc_max, RA, RB, RC, &
    1151             :                             zeta, zetb, zetc, a_mm, hmat, h_inv, vol, &
    1152             :                             G_bounds_1, R_bounds_2, R_bounds_3, &
    1153             :                             G_rads_1, R_rads_2, R_rads_3, &
    1154             :                             method, method_out, order)
    1155             :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT)     :: S_G
    1156             :       INTEGER, INTENT(IN)                                :: la_max, lb_max, lc_max
    1157             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: RA, RB, RC
    1158             :       REAL(KIND=dp), INTENT(IN)                          :: zeta, zetb, zetc, a_mm
    1159             :       REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN)         :: hmat, h_inv
    1160             :       REAL(KIND=dp), INTENT(IN)                          :: vol
    1161             :       REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN)         :: G_bounds_1, R_bounds_2
    1162             :       REAL(KIND=dp), DIMENSION(2, 3), INTENT(IN)         :: R_bounds_3
    1163             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: G_rads_1, R_rads_2
    1164             :       REAL(KIND=dp), DIMENSION(2), INTENT(IN)            :: R_rads_3
    1165             :       INTEGER, INTENT(IN)                                :: method
    1166             :       INTEGER, INTENT(OUT), OPTIONAL                     :: method_out
    1167             :       INTEGER, INTENT(IN), OPTIONAL                      :: order
    1168             : 
    1169             :       INTEGER                                            :: l_max, m_max, n_max, sum_method, &
    1170             :                                                             sum_order
    1171             :       LOGICAL                                            :: assert_stm
    1172             :       REAL(KIND=dp)                                      :: alpha, beta, G_rad, gamma, R_rad
    1173      148500 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: S_G_tmp
    1174             :       REAL(KIND=dp), DIMENSION(3)                        :: G_bound, R1, R2, R_bound
    1175             :       REAL(KIND=dp), DIMENSION(3, 3)                     :: ht
    1176             : 
    1177      148500 :       IF (PRESENT(order)) THEN
    1178           0 :          sum_order = order
    1179             :       ELSE
    1180      148500 :          sum_order = 0
    1181             :       ENDIF
    1182             : 
    1183      148500 :       sum_method = method
    1184             : 
    1185      148500 :       alpha = 0.25_dp/zeta
    1186      148500 :       beta = 0.25_dp/zetb
    1187      148500 :       gamma = 0.25_dp/zetc + a_mm
    1188             : 
    1189      148500 :       l_max = la_max
    1190      148500 :       m_max = lb_max
    1191      148500 :       n_max = lc_max
    1192             : 
    1193      594000 :       R1 = RB - RA
    1194      594000 :       R2 = RC - RB
    1195             : 
    1196     1930500 :       ht = twopi*TRANSPOSE(h_inv)
    1197             : 
    1198           0 :       SELECT CASE (sum_method)
    1199             :       CASE (1) ! sum_G sum_G' C(G) C(G,G') C(G')
    1200             : 
    1201           0 :          IF (sum_order .EQ. 0) THEN
    1202           0 :             sum_order = MAXLOC(G_bounds_1(:, 1), DIM=1)
    1203             :             assert_stm = MINLOC(G_bounds_1(:, 1), DIM=1) == &
    1204             :                MINLOC(G_bounds_1(:, 2), DIM=1) .AND. &
    1205             :                MINLOC(G_bounds_1(:, 1), DIM=1) == &
    1206           0 :                MINLOC(G_bounds_1(:, 3), DIM=1)
    1207           0 :             CPASSERT(assert_stm)
    1208             :          ENDIF
    1209             : 
    1210           0 :          CALL pgf_sum_3c_gspace_3d(S_G, l_max, m_max, n_max, R1, R2, alpha, beta, gamma, ht, vol, G_bounds_1, G_rads_1, sum_order)
    1211             : 
    1212             :       CASE (2) ! sum_G sum_R C(G) C(R)
    1213        5846 :          IF (sum_order .EQ. 0) THEN
    1214       23384 :             sum_order = MINLOC(G_bounds_1(:, 1), DIM=1)
    1215             :             assert_stm = MINLOC(G_bounds_1(:, 1), DIM=1) == &
    1216             :                MINLOC(G_bounds_1(:, 2), DIM=1) .AND. &
    1217             :                MINLOC(G_bounds_1(:, 1), DIM=1) == &
    1218       93536 :                MINLOC(G_bounds_1(:, 3), DIM=1)
    1219           0 :             CPASSERT(assert_stm)
    1220             :          ENDIF
    1221        5846 :          R_rad = R_rads_2(sum_order)
    1222        5846 :          G_rad = G_rads_1(sum_order)
    1223       23384 :          R_bound = R_bounds_2(sum_order, :)
    1224       23384 :          G_bound = G_bounds_1(sum_order, :)
    1225      142654 :          SELECT CASE (sum_order)
    1226             :          CASE (1)
    1227           0 :             ALLOCATE (S_G_tmp(ncoset(l_max), ncoset(m_max), ncoset(n_max)))
    1228             :             CALL pgf_sum_product_3c_gspace_3d(S_G_tmp, l_max, m_max, n_max, R1, R2, alpha, beta, gamma, hmat, h_inv, vol, &
    1229           0 :                                               R_bound, G_bound, R_rad, G_rad)
    1230           0 :             S_G = RESHAPE(S_G_tmp, SHAPE(S_G), order=[1, 2, 3])
    1231             :          CASE (2)
    1232           0 :             ALLOCATE (S_G_tmp(ncoset(m_max), ncoset(l_max), ncoset(n_max)))
    1233             :             CALL pgf_sum_product_3c_gspace_3d(S_G_tmp, m_max, l_max, n_max, -R1, R1 + R2, beta, alpha, gamma, hmat, h_inv, vol, &
    1234           0 :                                               R_bound, G_bound, R_rad, G_rad)
    1235           0 :             S_G = RESHAPE(S_G_tmp, SHAPE(S_G), order=[2, 1, 3])
    1236             :          CASE (3)
    1237       29230 :             ALLOCATE (S_G_tmp(ncoset(n_max), ncoset(m_max), ncoset(l_max)))
    1238             :             CALL pgf_sum_product_3c_gspace_3d(S_G_tmp, n_max, m_max, l_max, -R2, -R1, gamma, beta, alpha, hmat, h_inv, vol, &
    1239       40922 :                                               R_bound, G_bound, R_rad, G_rad)
    1240       23384 :             S_G = RESHAPE(S_G_tmp, SHAPE(S_G), order=[3, 2, 1])
    1241             :          END SELECT
    1242             :       CASE (3) ! sum_R sum_R' H(R, R')
    1243             :          CALL pgf_sum_3c_rspace_3d(S_G, l_max, m_max, n_max, RA, RB, RC, zeta, zetb, zetc, a_mm, hmat, h_inv, &
    1244      142654 :                                    R_bounds_3, R_rads_3)
    1245    24412738 :          S_G = S_G*pi**(-1.5_dp)*((zeta + zetb)/(zeta*zetb))**(-1.5_dp)
    1246             :       END SELECT
    1247             : 
    1248      148500 :       IF (PRESENT(method_out)) method_out = sum_method
    1249             : 
    1250      148500 :    END SUBROUTINE pgf_sum_3c_3d
    1251             : 
    1252             : ! **************************************************************************************************
    1253             : !> \brief Roughly estimated number of floating point operations
    1254             : !> \param ns_G1 ...
    1255             : !> \param ns_G2 ...
    1256             : !> \param l ...
    1257             : !> \param m ...
    1258             : !> \param n ...
    1259             : !> \return ...
    1260             : ! **************************************************************************************************
    1261      148500 :    PURE FUNCTION nsum_3c_gspace_3d(ns_G1, ns_G2, l, m, n)
    1262             :       REAL(KIND=dp), INTENT(IN)                          :: ns_G1, ns_G2
    1263             :       INTEGER, INTENT(IN)                                :: l, m, n
    1264             :       INTEGER(KIND=int_8)                                :: nsum_3c_gspace_3d
    1265             : 
    1266      148500 :       nsum_3c_gspace_3d = NINT(ns_G1*ns_G2*(5*exp_w + ncoset(l)*ncoset(m)*ncoset(n)*4), KIND=int_8)
    1267             : 
    1268      148500 :    END FUNCTION
    1269             : 
    1270             : ! **************************************************************************************************
    1271             : !> \brief ...
    1272             : !> \param S_G ...
    1273             : !> \param l_max ...
    1274             : !> \param m_max ...
    1275             : !> \param n_max ...
    1276             : !> \param R1 ...
    1277             : !> \param R2 ...
    1278             : !> \param alpha ...
    1279             : !> \param beta ...
    1280             : !> \param gamma ...
    1281             : !> \param ht ...
    1282             : !> \param vol ...
    1283             : !> \param G_c ...
    1284             : !> \param G_rad ...
    1285             : !> \param sum_order ...
    1286             : !> \param coulomb ...
    1287             : ! **************************************************************************************************
    1288           0 :    PURE SUBROUTINE pgf_sum_3c_gspace_3d(S_G, l_max, m_max, n_max, R1, R2, alpha, beta, gamma, &
    1289             :                                         ht, vol, G_c, G_rad, sum_order, coulomb)
    1290             :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT)     :: S_G
    1291             :       INTEGER, INTENT(IN)                                :: l_max, m_max, n_max
    1292             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: R1, R2
    1293             :       REAL(KIND=dp), INTENT(IN)                          :: alpha, beta, gamma
    1294             :       REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN)         :: ht
    1295             :       REAL(KIND=dp), INTENT(IN)                          :: vol
    1296             :       REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN)         :: G_c
    1297             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: G_rad
    1298             :       INTEGER, INTENT(IN)                                :: sum_order
    1299             :       LOGICAL, INTENT(IN), OPTIONAL                      :: coulomb
    1300             : 
    1301             :       INTEGER, DIMENSION(3)                              :: G1c, G2c, G3c
    1302             :       INTEGER                                            :: g1x, g1y, g1z, g2x, g2y, g2z, g3x, g3y, &
    1303             :                                                             g3z
    1304             :       COMPLEX(KIND=dp), DIMENSION(ncoset(l_max), ncoset( &
    1305           0 :                                   m_max), ncoset(n_max))                          :: S_G_c
    1306             :       LOGICAL                                            :: use_coulomb
    1307             :       REAL(KIND=dp)                                      :: G1_sq, G2_sq, G3_sq
    1308             :       REAL(KIND=dp), DIMENSION(3)                        :: G1, G1_x, G1_y, G1_z, G2, G2_x, G2_y, &
    1309             :                                                             G2_z, G3, G3_x, G3_y, G3_z, G_rads_sq
    1310             : 
    1311           0 :       S_G_c(:, :, :) = 0.0_dp
    1312             : 
    1313           0 :       G1c = FLOOR(G_c(1, :))
    1314           0 :       G2c = FLOOR(G_c(2, :))
    1315           0 :       G3c = FLOOR(G_c(3, :))
    1316             : 
    1317             :       ! we can not precompute exponentials for general cell
    1318             :       ! Perform double G sum
    1319           0 :       G_rads_sq = G_rad**2
    1320             : 
    1321           0 :       IF (PRESENT(coulomb)) THEN
    1322           0 :          use_coulomb = coulomb
    1323             :       ELSE
    1324           0 :          use_coulomb = .FALSE.
    1325             :       ENDIF
    1326             : 
    1327           0 :       SELECT CASE (sum_order)
    1328             :       CASE (1)
    1329           0 :          DO g2x = -G2c(1), G2c(1)
    1330           0 :             G2_x = ht(:, 1)*g2x
    1331           0 :             DO g2y = -G2c(2), G2c(2)
    1332           0 :                G2_y = ht(:, 2)*g2y
    1333           0 :                DO g2z = -G2c(3), G2c(3)
    1334           0 :                   G2_z = ht(:, 3)*g2z
    1335           0 :                   G2 = G2_x + G2_y + G2_z
    1336           0 :                   G2_sq = G2(1)**2 + G2(2)**2 + G2(3)**2
    1337           0 :                   IF (G2_sq .GT. G_rads_sq(2)) CYCLE
    1338           0 :                   DO g3x = -G3c(1), G3c(1)
    1339           0 :                      G3_x = ht(:, 1)*g3x
    1340           0 :                      DO g3y = -G3c(2), G3c(2)
    1341           0 :                         G3_y = ht(:, 2)*g3y
    1342           0 :                         DO g3z = -G3c(3), G3c(3)
    1343           0 :                            G3_z = ht(:, 3)*g3z
    1344           0 :                            G3 = G3_x + G3_y + G3_z
    1345           0 :                            G3_sq = G3(1)**2 + G3(2)**2 + G3(3)**2
    1346           0 :                            IF (G3_sq .GT. G_rads_sq(3)) CYCLE
    1347           0 :                            G1 = G3 - G2
    1348           0 :                            G1_sq = G1(1)**2 + G1(2)**2 + G1(3)**2
    1349           0 :                            IF (G1_sq .GT. G_rads_sq(1)) CYCLE
    1350           0 :                            IF (use_coulomb) THEN
    1351           0 :                               IF (g3x == 0 .AND. g3y == 0 .AND. g3z == 0) CYCLE
    1352             :                            ENDIF
    1353             :                            CALL pgf_product_3c_gspace_3d(S_G_c, G1, G1_sq, G2, G2_sq, G3, G3_sq, l_max, m_max, n_max, &
    1354           0 :                                                          alpha, beta, gamma, R1, R2, use_coulomb)
    1355             :                         ENDDO
    1356             :                      ENDDO
    1357             :                   ENDDO
    1358             :                ENDDO
    1359             :             ENDDO
    1360             :          ENDDO
    1361             :       CASE (2)
    1362           0 :          DO g1x = -G1c(1), G1c(1)
    1363           0 :             G1_x = ht(:, 1)*g1x
    1364           0 :             DO g1y = -G1c(2), G1c(2)
    1365           0 :                G1_y = ht(:, 2)*g1y
    1366           0 :                DO g1z = -G1c(3), G1c(3)
    1367           0 :                   G1_z = ht(:, 3)*g1z
    1368           0 :                   G1 = G1_x + G1_y + G1_z
    1369           0 :                   G1_sq = G1(1)**2 + G1(2)**2 + G1(3)**2
    1370           0 :                   IF (G1_sq .GT. G_rads_sq(1)) CYCLE
    1371           0 :                   DO g3x = -G3c(1), G3c(1)
    1372           0 :                      G3_x = ht(:, 1)*g3x
    1373           0 :                      DO g3y = -G3c(2), G3c(2)
    1374           0 :                         G3_y = ht(:, 2)*g3y
    1375           0 :                         DO g3z = -G3c(3), G3c(3)
    1376           0 :                            G3_z = ht(:, 3)*g3z
    1377           0 :                            G3 = G3_x + G3_y + G3_z
    1378           0 :                            G3_sq = G3(1)**2 + G3(2)**2 + G3(3)**2
    1379           0 :                            IF (G3_sq .GT. G_rads_sq(3)) CYCLE
    1380           0 :                            G2 = G3 - G1
    1381           0 :                            G2_sq = G2(1)**2 + G2(2)**2 + G2(3)**2
    1382           0 :                            IF (G2_sq .GT. G_rads_sq(2)) CYCLE
    1383           0 :                            IF (use_coulomb) THEN
    1384           0 :                               IF (g3x == 0 .AND. g3y == 0 .AND. g3z == 0) CYCLE
    1385             :                            ENDIF
    1386             :                            CALL pgf_product_3c_gspace_3d(S_G_c, G1, G1_sq, G2, G2_sq, G3, G3_sq, l_max, m_max, n_max, &
    1387           0 :                                                          alpha, beta, gamma, R1, R2, use_coulomb)
    1388             :                         ENDDO
    1389             :                      ENDDO
    1390             :                   ENDDO
    1391             :                ENDDO
    1392             :             ENDDO
    1393             :          ENDDO
    1394             :       CASE (3)
    1395           0 :          DO g1x = -G1c(1), G1c(1)
    1396           0 :             G1_x = ht(:, 1)*g1x
    1397           0 :             DO g1y = -G1c(2), G1c(2)
    1398           0 :                G1_y = ht(:, 2)*g1y
    1399           0 :                DO g1z = -G1c(3), G1c(3)
    1400           0 :                   G1_z = ht(:, 3)*g1z
    1401           0 :                   G1 = G1_x + G1_y + G1_z
    1402           0 :                   G1_sq = G1(1)**2 + G1(2)**2 + G1(3)**2
    1403           0 :                   IF (G1_sq .GT. G_rads_sq(1)) CYCLE
    1404           0 :                   DO g2x = -G2c(1), G2c(1)
    1405           0 :                      G2_x = ht(:, 1)*g2x
    1406           0 :                      DO g2y = -G2c(2), G2c(2)
    1407           0 :                         G2_y = ht(:, 2)*g2y
    1408           0 :                         DO g2z = -G2c(3), G2c(3)
    1409           0 :                            G2_z = ht(:, 3)*g2z
    1410           0 :                            G2 = G2_x + G2_y + G2_z
    1411           0 :                            G2_sq = G2(1)**2 + G2(2)**2 + G2(3)**2
    1412           0 :                            IF (G2_sq .GT. G_rads_sq(2)) CYCLE
    1413           0 :                            G3 = G1 + G2
    1414           0 :                            G3_sq = G3(1)**2 + G3(2)**2 + G3(3)**2
    1415           0 :                            IF (G3_sq .GT. G_rads_sq(3)) CYCLE
    1416           0 :                            IF (use_coulomb) THEN
    1417           0 :                               IF (g1x + g2x == 0 .AND. g1y + g2y == 0 .AND. g1z + g2z == 0) CYCLE
    1418             :                            ENDIF
    1419             :                            CALL pgf_product_3c_gspace_3d(S_G_c, G1, G1_sq, G2, G2_sq, G3, G3_sq, l_max, m_max, n_max, &
    1420           0 :                                                          alpha, beta, gamma, R1, R2, use_coulomb)
    1421             :                         ENDDO
    1422             :                      ENDDO
    1423             :                   ENDDO
    1424             :                ENDDO
    1425             :             ENDDO
    1426             :          ENDDO
    1427             :       END SELECT
    1428           0 :       S_G = REAL(S_G_c, KIND=dp)/vol**2
    1429             : 
    1430           0 :    END SUBROUTINE
    1431             : 
    1432             : ! **************************************************************************************************
    1433             : !> \brief ...
    1434             : !> \param S_G ...
    1435             : !> \param G1 ...
    1436             : !> \param G1_sq ...
    1437             : !> \param G2 ...
    1438             : !> \param G2_sq ...
    1439             : !> \param G3 ...
    1440             : !> \param G3_sq ...
    1441             : !> \param l_max ...
    1442             : !> \param m_max ...
    1443             : !> \param n_max ...
    1444             : !> \param alpha ...
    1445             : !> \param beta ...
    1446             : !> \param gamma ...
    1447             : !> \param R1 ...
    1448             : !> \param R2 ...
    1449             : !> \param coulomb ...
    1450             : ! **************************************************************************************************
    1451           0 :    PURE SUBROUTINE pgf_product_3c_gspace_3d(S_G, G1, G1_sq, G2, G2_sq, G3, G3_sq, l_max, m_max, n_max, &
    1452             :                                             alpha, beta, gamma, R1, R2, coulomb)
    1453             : 
    1454             :       COMPLEX(KIND=dp), DIMENSION(:, :, :), INTENT(OUT)  :: S_G
    1455             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: G1
    1456             :       REAL(KIND=dp), INTENT(IN)                          :: G1_sq
    1457             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: G2
    1458             :       REAL(KIND=dp), INTENT(IN)                          :: G2_sq
    1459             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: G3
    1460             :       REAL(KIND=dp), INTENT(IN)                          :: G3_sq
    1461             :       INTEGER, INTENT(IN)                                :: l_max, m_max, n_max
    1462             :       REAL(KIND=dp), INTENT(IN)                          :: alpha, beta, gamma
    1463             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: R1, R2
    1464             :       LOGICAL, INTENT(IN)                                :: coulomb
    1465             : 
    1466             :       COMPLEX(KIND=dp)                                   :: exp_tot
    1467             :       INTEGER                                            :: k, l, lco, lx, ly, lz, m, mco, mx, my, &
    1468             :                                                             mz, n, nco, nx, ny, nz
    1469           0 :       COMPLEX(KIND=dp), DIMENSION(ncoset(n_max))         :: S_G_n
    1470           0 :       COMPLEX(KIND=dp), DIMENSION(ncoset(m_max))         :: S_G_m
    1471           0 :       COMPLEX(KIND=dp), DIMENSION(ncoset(l_max))         :: S_G_l
    1472           0 :       REAL(KIND=dp), DIMENSION(3, 0:l_max)               :: G1_pow_l
    1473           0 :       REAL(KIND=dp), DIMENSION(3, 0:m_max)               :: G2_pow_m
    1474           0 :       REAL(KIND=dp), DIMENSION(3, 0:n_max)               :: G3_pow_n
    1475             : 
    1476             :       exp_tot = EXP(gaussi*DOT_PRODUCT(G1, R1))*EXP(-alpha*G1_sq)* & ! cost: 5 exp_w flops
    1477             :                 EXP(-beta*G2_sq)* &
    1478           0 :                 EXP(-gamma*G3_sq)*EXP(gaussi*DOT_PRODUCT(G3, R2))
    1479             : 
    1480           0 :       IF (coulomb) exp_tot = exp_tot/G3_sq
    1481             : 
    1482           0 :       DO k = 1, 3
    1483           0 :          G1_pow_l(k, 0) = 1.0_dp
    1484           0 :          DO l = 1, l_max
    1485           0 :             G1_pow_l(k, l) = G1_pow_l(k, l - 1)*G1(k)
    1486             :          ENDDO
    1487           0 :          G2_pow_m(k, 0) = 1.0_dp
    1488           0 :          DO m = 1, m_max
    1489           0 :             G2_pow_m(k, m) = G2_pow_m(k, m - 1)*G2(k)
    1490             :          ENDDO
    1491           0 :          G3_pow_n(k, 0) = 1.0_dp
    1492           0 :          DO n = 1, n_max
    1493           0 :             G3_pow_n(k, n) = G3_pow_n(k, n - 1)*G3(k)
    1494             :          ENDDO
    1495             :       ENDDO
    1496             : 
    1497           0 :       DO lco = 1, ncoset(l_max)
    1498           0 :          CALL get_l(lco, l, lx, ly, lz)
    1499           0 :          S_G_l(lco) = G1_pow_l(1, lx)*G1_pow_l(2, ly)*G1_pow_l(3, lz)*(-1.0_dp)**l*i_pow(l)
    1500             :       ENDDO
    1501             : 
    1502           0 :       DO mco = 1, ncoset(m_max)
    1503           0 :          CALL get_l(mco, m, mx, my, mz)
    1504           0 :          S_G_m(mco) = G2_pow_m(1, mx)*G2_pow_m(2, my)*G2_pow_m(3, mz)*(-1.0_dp)**m*i_pow(m)
    1505             :       ENDDO
    1506             : 
    1507           0 :       DO nco = 1, ncoset(n_max)
    1508           0 :          CALL get_l(nco, n, nx, ny, nz)
    1509           0 :          S_G_n(nco) = G3_pow_n(1, nx)*G3_pow_n(2, ny)*G3_pow_n(3, nz)*i_pow(n)
    1510             :       ENDDO
    1511             : 
    1512           0 :       DO nco = 1, ncoset(n_max)
    1513           0 :          DO mco = 1, ncoset(m_max)
    1514           0 :             DO lco = 1, ncoset(l_max)
    1515             :                S_G(lco, mco, nco) = S_G(lco, mco, nco) + & ! cost: 4 flops
    1516             :                                     S_G_l(lco)*S_G_m(mco)*S_G_n(nco)* &
    1517           0 :                                     exp_tot
    1518             :             ENDDO
    1519             :          ENDDO
    1520             :       ENDDO
    1521             : 
    1522           0 :    END SUBROUTINE pgf_product_3c_gspace_3d
    1523             : 
    1524             : ! **************************************************************************************************
    1525             : !> \brief Roughly estimated number of floating point operations
    1526             : !> \param ns_G ...
    1527             : !> \param ns_R ...
    1528             : !> \param l ...
    1529             : !> \param m ...
    1530             : !> \param n ...
    1531             : !> \return ...
    1532             : ! **************************************************************************************************
    1533      148500 :    PURE FUNCTION nsum_product_3c_gspace_3d(ns_G, ns_R, l, m, n)
    1534             :       REAL(KIND=dp), INTENT(IN)                          :: ns_G, ns_R
    1535             :       INTEGER, INTENT(IN)                                :: l, m, n
    1536             :       INTEGER(KIND=int_8)                                :: nsum_product_3c_gspace_3d
    1537             : 
    1538             :       nsum_product_3c_gspace_3d = &
    1539             :          NINT( &
    1540             :          ns_G*( &
    1541             :          (exp_w*2) + &
    1542             :          ns_R*(exp_w*2 + ncoset(l + m)*7) + &
    1543             :          3*nsum_gaussian_overlap(l, m, 1) + &
    1544             :          ncoset(l)*ncoset(m)*(ncoset(l + m)*4 + ncoset(n)*8)), &
    1545      148500 :          KIND=int_8)
    1546      148500 :    END FUNCTION
    1547             : 
    1548             : ! **************************************************************************************************
    1549             : !> \brief ...
    1550             : !> \param S_G ...
    1551             : !> \param l_max ...
    1552             : !> \param m_max ...
    1553             : !> \param n_max ...
    1554             : !> \param R1 ...
    1555             : !> \param R2 ...
    1556             : !> \param alpha ...
    1557             : !> \param beta ...
    1558             : !> \param gamma ...
    1559             : !> \param hmat ...
    1560             : !> \param h_inv ...
    1561             : !> \param vol ...
    1562             : !> \param R_c ...
    1563             : !> \param G_c ...
    1564             : !> \param R_rad ...
    1565             : !> \param G_rad ...
    1566             : ! **************************************************************************************************
    1567        5846 :    PURE SUBROUTINE pgf_sum_product_3c_gspace_3d(S_G, l_max, m_max, n_max, R1, R2, alpha, beta, gamma, &
    1568             :                                                 hmat, h_inv, vol, R_c, G_c, R_rad, G_rad)
    1569             :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT)     :: S_G
    1570             :       INTEGER, INTENT(IN)                                :: l_max, m_max, n_max
    1571             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: R1, R2
    1572             :       REAL(KIND=dp), INTENT(IN)                          :: alpha, beta, gamma
    1573             :       REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN)         :: hmat, h_inv
    1574             :       REAL(KIND=dp), INTENT(IN)                          :: vol
    1575             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: R_c, G_c
    1576             :       REAL(KIND=dp), INTENT(IN)                          :: R_rad, G_rad
    1577             : 
    1578             :       COMPLEX(KIND=dp)                                   :: exp_tot
    1579             :       INTEGER                                            :: gx, gy, gz, k, l, lco, lnco, lx, ly, lz, &
    1580             :                                                             m, mco, n, nco
    1581             :       COMPLEX(KIND=dp), &
    1582       11692 :          DIMENSION(ncoset(m_max), ncoset(n_max))         :: S_R
    1583             :       COMPLEX(KIND=dp), DIMENSION(ncoset(l_max), ncoset( &
    1584       11692 :                                   m_max), ncoset(n_max))                          :: S_G_c
    1585             :       REAL(KIND=dp)                                      :: G_rad_sq, G_sq, R_rad_sq
    1586             :       REAL(KIND=dp), DIMENSION(3)                        :: G, G_x, G_y, G_z
    1587       11692 :       REAL(KIND=dp), DIMENSION(3, 0:l_max)               :: G_pow_l
    1588             :       REAL(KIND=dp), DIMENSION(3, 3)                     :: ht
    1589        5846 :       REAL(KIND=dp), DIMENSION(ncoset(l_max))            :: S_G_c_l
    1590             : 
    1591      276845 :       S_G_c(:, :, :) = 0.0_dp
    1592             : 
    1593        5846 :       G_rad_sq = G_rad**2
    1594        5846 :       R_rad_sq = R_rad**2
    1595             : 
    1596        5846 :       lnco = ncoset(l_max)
    1597             : 
    1598       75998 :       ht = twopi*TRANSPOSE(h_inv)
    1599       36078 :       DO gx = -FLOOR(G_c(1)), FLOOR(G_c(1))
    1600      120928 :          G_x = ht(:, 1)*gx
    1601      194252 :          DO gy = -FLOOR(G_c(2)), FLOOR(G_c(2))
    1602      632696 :             G_y = ht(:, 2)*gy
    1603     1028374 :             DO gz = -FLOOR(G_c(3)), FLOOR(G_c(3))
    1604     3359872 :                G_z = ht(:, 3)*gz
    1605     3359872 :                G = G_x + G_y + G_z
    1606      839968 :                G_sq = G(1)**2 + G(2)**2 + G(3)**2
    1607      839968 :                IF (G_sq .GT. G_rad_sq) CYCLE
    1608             : 
    1609     1350880 :                exp_tot = EXP(-alpha*G_sq)*EXP(gaussi*DOT_PRODUCT(G, R1)) ! cost: exp_w*2 flops
    1610             : 
    1611     1350880 :                DO k = 1, 3
    1612     1013160 :                   G_pow_l(k, 0) = 1.0_dp
    1613     2236258 :                   DO l = 1, l_max
    1614     1898538 :                      G_pow_l(k, l) = G_pow_l(k, l - 1)*G(k)
    1615             :                   ENDDO
    1616             :                ENDDO
    1617             : 
    1618      337720 :                CALL pgf_sum_product_3c_rspace_3d(S_R, m_max, n_max, G, R2, beta, gamma, hmat, h_inv, vol, R_c, R_rad_sq)
    1619             : 
    1620     1947011 :                DO lco = 1, ncoset(l_max)
    1621     1609291 :                   CALL get_l(lco, l, lx, ly, lz)
    1622     1947011 :                   S_G_c_l(lco) = G_pow_l(1, lx)*G_pow_l(2, ly)*G_pow_l(3, lz)*(-1.0_dp)**l
    1623             :                ENDDO
    1624             : 
    1625     1519228 :                DO nco = 1, ncoset(n_max)
    1626     1023334 :                   CALL get_l(nco, n)
    1627     4612179 :                   DO mco = 1, ncoset(m_max)
    1628     2748877 :                      CALL get_l(mco, m)
    1629    15176745 :                      DO lco = 1, ncoset(l_max)
    1630    11404534 :                         CALL get_l(lco, l)
    1631             :                         S_G_c(lco, mco, nco) = & ! cost: 8 flops
    1632             :                            S_G_c(lco, mco, nco) + &
    1633             :                            S_G_c_l(lco)* &
    1634    14153411 :                            exp_tot*i_pow(l + m + n)*(-1.0_dp)**m*S_R(mco, nco)
    1635             :                      ENDDO
    1636             :                   ENDDO
    1637             :                ENDDO
    1638             :             ENDDO
    1639             :          ENDDO
    1640             :       ENDDO
    1641      276845 :       S_G = REAL(S_G_c, KIND=dp)/vol**2
    1642             : 
    1643        5846 :    END SUBROUTINE pgf_sum_product_3c_gspace_3d
    1644             : 
    1645             : ! **************************************************************************************************
    1646             : !> \brief ...
    1647             : !> \param S_R ...
    1648             : !> \param l_max ...
    1649             : !> \param m_max ...
    1650             : !> \param G ...
    1651             : !> \param R ...
    1652             : !> \param alpha ...
    1653             : !> \param beta ...
    1654             : !> \param hmat ...
    1655             : !> \param h_inv ...
    1656             : !> \param vol ...
    1657             : !> \param R_c ...
    1658             : !> \param R_rad_sq ...
    1659             : ! **************************************************************************************************
    1660      337720 :    PURE SUBROUTINE pgf_sum_product_3c_rspace_3d(S_R, l_max, m_max, G, R, alpha, beta, hmat, h_inv, vol, R_c, R_rad_sq)
    1661             :       COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(OUT)     :: S_R
    1662             :       INTEGER, INTENT(IN)                                :: l_max, m_max
    1663             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: G, R
    1664             :       REAL(KIND=dp), INTENT(IN)                          :: alpha, beta
    1665             :       REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN)         :: hmat, h_inv
    1666             :       REAL(KIND=dp), INTENT(IN)                          :: vol
    1667             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: R_c
    1668             :       REAL(KIND=dp), INTENT(IN)                          :: R_rad_sq
    1669             : 
    1670             :       COMPLEX(KIND=dp)                                   :: exp_tot
    1671             :       INTEGER                                            :: k, l, lco, lx, ly, lz, m, mco, mx, my, &
    1672             :                                                             mz, sx, sy, sz, t, tco, tx, ty, tz
    1673      675440 :       COMPLEX(KIND=dp), DIMENSION(ncoset(l_max + m_max))   :: S_R_t
    1674             :       REAL(KIND=dp)                                      :: c1, c2, Rp_sq
    1675             :       REAL(KIND=dp), &
    1676      675440 :          DIMENSION(-1:l_max + m_max + 1, -1:l_max, -1:m_max) :: E1, E2, E3
    1677             :       REAL(KIND=dp), DIMENSION(3)                        :: R_l, R_r, Rp, Rx, Ry, Rz, s_shift
    1678      675440 :       REAL(KIND=dp), DIMENSION(3, 0:l_max + m_max)         :: R_pow_t
    1679             : 
    1680      337720 :       c1 = 0.25_dp/(alpha + beta)
    1681      337720 :       c2 = alpha/(alpha + beta)
    1682             : 
    1683     2198096 :       S_R_t(:) = 0.0_dp
    1684     4109931 :       S_R(:, :) = 0.0_dp
    1685             : 
    1686     4390360 :       s_shift = MATMUL(h_inv, R)
    1687     1350880 :       R_l = -R_c + s_shift
    1688     1350880 :       R_r = R_c + s_shift
    1689             : 
    1690     1084696 :       DO sx = CEILING(R_l(1)), FLOOR(R_r(1))
    1691     2987904 :          Rx = hmat(:, 1)*sx
    1692     3046393 :          DO sy = CEILING(R_l(2)), FLOOR(R_r(2))
    1693     7846788 :             Ry = hmat(:, 2)*sy
    1694     8205714 :             DO sz = CEILING(R_l(3)), FLOOR(R_r(3))
    1695    21988164 :                Rz = hmat(:, 3)*sz
    1696    21988164 :                Rp = Rx + Ry + Rz - R
    1697     5497041 :                Rp_sq = Rp(1)**2 + Rp(2)**2 + Rp(3)**2
    1698     5497041 :                IF (Rp_sq .GT. R_rad_sq) CYCLE
    1699             : 
    1700     8361484 :                exp_tot = EXP(-c1*Rp_sq)*EXP(-gaussi*c2*DOT_PRODUCT(Rp, G)) ! cost: exp_w*2 flops
    1701     8361484 :                DO k = 1, 3
    1702     6271113 :                   R_pow_t(k, 0) = 1.0_dp
    1703    13292668 :                   DO t = 1, l_max + m_max
    1704    11202297 :                      R_pow_t(k, t) = R_pow_t(k, t - 1)*Rp(k)
    1705             :                   ENDDO
    1706             :                ENDDO
    1707             : 
    1708    13088657 :                DO tco = 1, ncoset(l_max + m_max)
    1709     9036589 :                   CALL get_l(tco, t, tx, ty, tz)
    1710             :                   S_R_t(tco) = S_R_t(tco) + & ! cost: 7 flops
    1711             :                                R_pow_t(1, tx)*R_pow_t(2, ty)*R_pow_t(3, tz)* &
    1712    14533630 :                                (-1.0_dp)**t*i_pow(t)*exp_tot
    1713             :                ENDDO
    1714             : 
    1715             :             ENDDO
    1716             :          ENDDO
    1717             :       ENDDO
    1718             : 
    1719      337720 :       CALL create_gaussian_overlap_dist_to_hermite(l_max, m_max, alpha, beta, G(1), 0.0_dp, 1, E1)
    1720      337720 :       CALL create_gaussian_overlap_dist_to_hermite(l_max, m_max, alpha, beta, G(2), 0.0_dp, 1, E2)
    1721      337720 :       CALL create_gaussian_overlap_dist_to_hermite(l_max, m_max, alpha, beta, G(3), 0.0_dp, 1, E3)
    1722             : 
    1723     1361054 :       DO mco = 1, ncoset(m_max)
    1724     1023334 :          CALL get_l(mco, m, mx, my, mz)
    1725     4109931 :          DO lco = 1, ncoset(l_max)
    1726     2748877 :             CALL get_l(lco, l, lx, ly, lz)
    1727     8010112 :             DO tx = 0, lx + mx
    1728    13267665 :             DO ty = 0, ly + my
    1729    19515063 :             DO tz = 0, lz + mz
    1730     8996275 :                tco = coset(tx, ty, tz)
    1731             :                S_R(lco, mco) = S_R(lco, mco) + & ! cost: 4 flops
    1732    15277162 :                                E1(tx, lx, mx)*E2(ty, ly, my)*E3(tz, lz, mz)*S_R_t(tco)
    1733             : 
    1734             :             ENDDO
    1735             :             ENDDO
    1736             :             ENDDO
    1737             :          ENDDO
    1738             :       ENDDO
    1739             : 
    1740     4109931 :       S_R(:, :) = S_R(:, :)*vol/(twopi)**3*(pi/(alpha + beta))**1.5_dp
    1741             : 
    1742      337720 :    END SUBROUTINE pgf_sum_product_3c_rspace_3d
    1743             : 
    1744             : ! **************************************************************************************************
    1745             : !> \brief Roughly estimated number of floating point operations
    1746             : !> \param ns_R1 ...
    1747             : !> \param ns_R2 ...
    1748             : !> \param l ...
    1749             : !> \param m ...
    1750             : !> \param n ...
    1751             : !> \return ...
    1752             : ! **************************************************************************************************
    1753      148500 :    PURE FUNCTION nsum_3c_rspace_3d(ns_R1, ns_R2, l, m, n)
    1754             :       REAL(KIND=dp), INTENT(IN)                          :: ns_R1, ns_R2
    1755             :       INTEGER, INTENT(IN)                                :: l, m, n
    1756             :       INTEGER(KIND=int_8)                                :: nsum_3c_rspace_3d
    1757             : 
    1758             :       nsum_3c_rspace_3d = &
    1759             :          NINT( &
    1760             :          ns_R1*( &
    1761             :          ns_R2*(exp_w + ncoset(l + m + n)*4) + &
    1762             :          3*nsum_gaussian_overlap(l, m, 2) + &
    1763             :          ncoset(m)*ncoset(l)*( &
    1764             :          ncoset(l + m)*2 + ncoset(n)*ncoset(l + m)*4)), &
    1765      148500 :          KIND=int_8)
    1766             : 
    1767      148500 :    END FUNCTION
    1768             : 
    1769             : ! **************************************************************************************************
    1770             : !> \brief ...
    1771             : !> \param S_R ...
    1772             : !> \param l_max ...
    1773             : !> \param m_max ...
    1774             : !> \param n_max ...
    1775             : !> \param RA ...
    1776             : !> \param RB ...
    1777             : !> \param RC ...
    1778             : !> \param zeta ...
    1779             : !> \param zetb ...
    1780             : !> \param zetc ...
    1781             : !> \param a_mm ...
    1782             : !> \param hmat ...
    1783             : !> \param h_inv ...
    1784             : !> \param R_c ...
    1785             : !> \param R_rad ...
    1786             : ! **************************************************************************************************
    1787      142654 :    SUBROUTINE pgf_sum_3c_rspace_3d(S_R, l_max, m_max, n_max, RA, RB, RC, zeta, zetb, zetc, a_mm, &
    1788             :                                    hmat, h_inv, R_c, R_rad)
    1789             :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT)     :: S_R
    1790             :       INTEGER, INTENT(IN)                                :: l_max, m_max, n_max
    1791             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: RA, RB, RC
    1792             :       REAL(KIND=dp), INTENT(IN)                          :: zeta, zetb, zetc, a_mm
    1793             :       REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN)         :: hmat, h_inv
    1794             :       REAL(KIND=dp), DIMENSION(2, 3), INTENT(IN)         :: R_c
    1795             :       REAL(KIND=dp), DIMENSION(2), INTENT(IN)            :: R_rad
    1796             : 
    1797             :       INTEGER                                            :: k, l, lco, lx, ly, lz, m, mco, mx, my, &
    1798             :                                                             mz, n, nco, nx, ny, nz, s1x, s1y, s1z, &
    1799             :                                                             s2x, s2y, s2z, t, tco, tnco, ttco, &
    1800             :                                                             ttx, tty, ttz, tx, ty, tz
    1801             :       REAL(KIND=dp)                                      :: alpha, exp_tot, R1_sq, R_sq
    1802      142654 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: h_to_c
    1803             :       REAL(KIND=dp), &
    1804      285308 :          DIMENSION(-1:l_max + m_max + 1, -1:l_max, -1:m_max) :: E1, E2, E3
    1805             :       REAL(KIND=dp), DIMENSION(2)                        :: R_rad_sq
    1806             :       REAL(KIND=dp), DIMENSION(3)                        :: R, R1, R1_l, R1_r, R1_tmp, R1x, R1y, &
    1807             :                                                             R1z, R2_l, R2_r, R2x, R2y, R2z, &
    1808             :                                                             R_offset, R_tmp, s1_shift, s2_shift
    1809      285308 :       REAL(KIND=dp), DIMENSION(3, 0:l_max + m_max + n_max)   :: R_pow_t
    1810             :       REAL(KIND=dp), DIMENSION(ncoset(l_max + m_max), &
    1811      285308 :                                ncoset(l_max), ncoset(m_max))                   :: Eco
    1812             :       REAL(KIND=dp), &
    1813      285308 :          DIMENSION(ncoset(l_max + m_max + n_max))            :: S_R_t
    1814             :       REAL(KIND=dp), DIMENSION(ncoset(l_max + m_max + n_max) &
    1815      285308 :                                , ncoset(l_max + m_max + n_max))                    :: h_to_c_tco
    1816             : 
    1817      142654 :       alpha = 1.0_dp/((zeta + zetb + zetc)/((zeta + zetb)*zetc) + 4.0_dp*a_mm)
    1818             : 
    1819    71800162 :       Eco(:, :, :) = 0.0_dp
    1820    24264238 :       S_R(:, :, :) = 0.0_dp
    1821   157194456 :       h_to_c_tco(:, :) = 0.0_dp
    1822             : 
    1823      570616 :       R_offset = RC - (zeta*RA + zetb*RB)/(zeta + zetb)
    1824             : 
    1825      142654 :       CALL create_hermite_to_cartesian(alpha, l_max + m_max + n_max, h_to_c)
    1826             : 
    1827     2682070 :       DO tco = 1, ncoset(l_max + m_max + n_max)
    1828     2539416 :          CALL get_l(tco, t, tx, ty, tz)
    1829     8078511 :          DO ttx = 0, tx
    1830    18802080 :          DO tty = 0, ty
    1831    37176724 :          DO ttz = 0, tz
    1832    20914060 :             ttco = coset(ttx, tty, ttz)
    1833    31780283 :             h_to_c_tco(ttco, tco) = h_to_c(ttx, tx)*h_to_c(tty, ty)*h_to_c(ttz, tz)
    1834             :          ENDDO
    1835             :          ENDDO
    1836             :          ENDDO
    1837             :       ENDDO
    1838             : 
    1839     2282464 :       s1_shift = MATMUL(h_inv, RA - RB)
    1840      570616 :       R1_l = -R_c(1, :) + s1_shift
    1841      570616 :       R1_r = R_c(1, :) + s1_shift
    1842             : 
    1843      427962 :       R_rad_sq = R_rad**2
    1844             : 
    1845      348524 :       DO s1x = CEILING(R1_l(1)), FLOOR(R1_r(1))
    1846      823480 :          R1x = hmat(:, 1)*s1x
    1847      741859 :          DO s1y = CEILING(R1_l(2)), FLOOR(R1_r(2))
    1848     1573340 :             R1y = hmat(:, 2)*s1y
    1849     1532744 :             DO s1z = CEILING(R1_l(3)), FLOOR(R1_r(3))
    1850     3734156 :                R1z = hmat(:, 3)*s1z
    1851     3734156 :                R1 = R1x + R1y + R1z
    1852    13798961 :                S_R_t(:) = 0.0_dp
    1853     3734156 :                R1_tmp = R1 - (RA - RB)
    1854      933539 :                R1_sq = R1_tmp(1)**2 + R1_tmp(2)**2 + R1_tmp(3)**2
    1855             : 
    1856      933539 :                IF (R1_sq .GT. R_rad_sq(1)) CYCLE
    1857             : 
    1858     1625580 :                R_tmp = R_offset + R1*zeta/(zeta + zetb)
    1859     6502320 :                s2_shift = -MATMUL(h_inv, R_tmp)
    1860     1625580 :                R2_l = -R_c(2, :) + s2_shift
    1861     1625580 :                R2_r = R_c(2, :) + s2_shift
    1862     1468530 :                DO s2x = CEILING(R2_l(1)), FLOOR(R2_r(1))
    1863     4248540 :                   R2x = hmat(:, 1)*s2x
    1864     6145785 :                   DO s2y = CEILING(R2_l(2)), FLOOR(R2_r(2))
    1865    18709020 :                      R2y = hmat(:, 2)*s2y
    1866    34911377 :                      DO s2z = CEILING(R2_l(3)), FLOOR(R2_r(3))
    1867   116687948 :                         R2z = hmat(:, 3)*s2z
    1868   116687948 :                         R = R_tmp + R2x + R2y + R2z
    1869    29171987 :                         R_sq = R(1)**2 + R(2)**2 + R(3)**2
    1870             : 
    1871    29171987 :                         IF (R_sq .GT. R_rad_sq(2)) CYCLE
    1872             : 
    1873    14475027 :                         exp_tot = EXP(-alpha*R_sq) ! cost: exp_w flops
    1874    57900108 :                         DO k = 1, 3
    1875    43425081 :                            R_pow_t(k, 0) = 1.0_dp
    1876   156067857 :                            DO t = 1, l_max + m_max + n_max
    1877   141592830 :                               R_pow_t(k, t) = R_pow_t(k, t - 1)*R(k)
    1878             :                            ENDDO
    1879             :                         ENDDO
    1880   272529341 :                         DO tco = 1, ncoset(l_max + m_max + n_max)
    1881   253377059 :                            CALL get_l(tco, t, tx, ty, tz)
    1882   282549046 :                            S_R_t(tco) = S_R_t(tco) + R_pow_t(1, tx)*R_pow_t(2, ty)*R_pow_t(3, tz)*exp_tot ! cost: 4 flops
    1883             :                         ENDDO
    1884             :                      ENDDO
    1885             :                   ENDDO
    1886             :                ENDDO
    1887             : 
    1888     5981328 :                S_R_t(:) = MATMUL(TRANSPOSE(h_to_c_tco), S_R_t)*(alpha/pi)**1.5_dp
    1889             : 
    1890      406395 :                CALL create_gaussian_overlap_dist_to_hermite(l_max, m_max, zeta, zetb, RA(1) - R1(1), RB(1), 2, E1)
    1891      406395 :                CALL create_gaussian_overlap_dist_to_hermite(l_max, m_max, zeta, zetb, RA(2) - R1(2), RB(2), 2, E2)
    1892      406395 :                CALL create_gaussian_overlap_dist_to_hermite(l_max, m_max, zeta, zetb, RA(3) - R1(3), RB(3), 2, E3)
    1893             : 
    1894     1416942 :                DO mco = 1, ncoset(m_max)
    1895     1010547 :                   CALL get_l(mco, m, mx, my, mz)
    1896     5054628 :                   DO lco = 1, ncoset(l_max)
    1897     3637686 :                      CALL get_l(lco, l, lx, ly, lz)
    1898    10890282 :                      DO tx = 0, lx + mx
    1899    20345458 :                      DO ty = 0, ly + my
    1900    33814276 :                      DO tz = 0, lz + mz
    1901    17106504 :                         tco = coset(tx, ty, tz)
    1902    27572227 :                         Eco(tco, lco, mco) = E1(tx, lx, mx)*E2(ty, ly, my)*E3(tz, lz, mz) ! cost: 2 flops
    1903             :                      ENDDO
    1904             :                      ENDDO
    1905             :                      ENDDO
    1906             :                   ENDDO
    1907             :                ENDDO
    1908             : 
    1909     3040511 :                DO nco = 1, ncoset(n_max)
    1910     2240781 :                   CALL get_l(nco, n, nx, ny, nz)
    1911    17172306 :                   DO tco = 1, ncoset(l_max + m_max)
    1912    14525130 :                      CALL get_l(tco, t, tx, ty, tz)
    1913    14525130 :                      tnco = coset(tx + nx, ty + ny, tz + nz)
    1914             :                      S_R(:, :, nco) = S_R(:, :, nco) + & ! cost: 4 flops
    1915             :                                       Eco(tco, :, :)* &
    1916  1070960721 :                                       (-1)**n*S_R_t(tnco)
    1917             : 
    1918             :                   ENDDO
    1919             :                ENDDO
    1920             :             ENDDO
    1921             :          ENDDO
    1922             :       ENDDO
    1923      142654 :    END SUBROUTINE pgf_sum_3c_rspace_3d
    1924             : 
    1925             : ! **************************************************************************************************
    1926             : !> \brief Compute bounding box for ellipsoid. This is needed in order to find summation bounds for
    1927             : !>        sphere for sums over non-orthogonal lattice vectors.
    1928             : !> \param s_rad sphere radius
    1929             : !> \param s_to_e sphere to ellipsoid trafo
    1930             : !> \return ...
    1931             : ! **************************************************************************************************
    1932    12295200 :    PURE FUNCTION ellipsoid_bounds(s_rad, s_to_e)
    1933             :       REAL(KIND=dp), INTENT(IN)                          :: s_rad
    1934             :       REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN)         :: s_to_e
    1935             :       REAL(KIND=dp), DIMENSION(3)                        :: ellipsoid_bounds
    1936             : 
    1937             :       INTEGER                                            :: i_xyz
    1938             : 
    1939    49180800 :       DO i_xyz = 1, 3
    1940    49180800 :          ellipsoid_bounds(i_xyz) = SQRT(s_to_e(i_xyz, 1)**2 + s_to_e(i_xyz, 2)**2 + s_to_e(i_xyz, 3)**2)*s_rad
    1941             :       ENDDO
    1942             : 
    1943             :    END FUNCTION ellipsoid_bounds
    1944             : 
    1945             : ! **************************************************************************************************
    1946             : !> \brief Roughly estimated number of floating point operations
    1947             : !> \param l ...
    1948             : !> \param m ...
    1949             : !> \param H_or_C_product ...
    1950             : !> \return ...
    1951             : ! **************************************************************************************************
    1952      297000 :    PURE FUNCTION nsum_gaussian_overlap(l, m, H_or_C_product)
    1953             :       INTEGER, INTENT(IN)                                :: l, m, H_or_C_product
    1954             :       INTEGER                                            :: nsum_gaussian_overlap
    1955             : 
    1956             :       INTEGER                                            :: loop
    1957             : 
    1958      297000 :       nsum_gaussian_overlap = exp_w
    1959      297000 :       loop = (m + 1)*(l + 1)*(m + l + 2)
    1960      297000 :       IF (H_or_C_product .EQ. 1) THEN
    1961      148500 :          nsum_gaussian_overlap = nsum_gaussian_overlap + loop*16
    1962             :       ELSE
    1963      148500 :          nsum_gaussian_overlap = nsum_gaussian_overlap + loop*32
    1964             :       ENDIF
    1965      297000 :    END FUNCTION
    1966             : 
    1967             : ! **************************************************************************************************
    1968             : !> \brief ...
    1969             : !> \param lco ...
    1970             : !> \param l ...
    1971             : !> \param lx ...
    1972             : !> \param ly ...
    1973             : !> \param lz ...
    1974             : ! **************************************************************************************************
    1975   669680053 :    PURE ELEMENTAL SUBROUTINE get_l(lco, l, lx, ly, lz)
    1976             :       INTEGER, INTENT(IN)                                :: lco
    1977             :       INTEGER, INTENT(OUT)                               :: l
    1978             :       INTEGER, INTENT(OUT), OPTIONAL                     :: lx, ly, lz
    1979             : 
    1980  2678720212 :       l = SUM(indco(:, lco))
    1981   669680053 :       IF (PRESENT(lx)) lx = indco(1, lco)
    1982   669680053 :       IF (PRESENT(ly)) ly = indco(2, lco)
    1983   669680053 :       IF (PRESENT(lz)) lz = indco(3, lco)
    1984   669680053 :    END SUBROUTINE
    1985             : 
    1986             : ! **************************************************************************************************
    1987             : !> \brief ...
    1988             : !> \param i ...
    1989             : !> \return ...
    1990             : ! **************************************************************************************************
    1991   333871597 :    PURE ELEMENTAL FUNCTION i_pow(i)
    1992             :       INTEGER, INTENT(IN)                                :: i
    1993             :       COMPLEX(KIND=dp)                                   :: i_pow
    1994             : 
    1995             :       COMPLEX(KIND=dp), DIMENSION(0:3), PARAMETER :: &
    1996             :          ip = (/(1.0_dp, 0.0_dp), (0.0_dp, 1.0_dp), (-1.0_dp, 0.0_dp), (0.0_dp, -1.0_dp)/)
    1997             : 
    1998   333871597 :       i_pow = ip(MOD(i, 4))
    1999             : 
    2000   333871597 :    END FUNCTION
    2001             : 
    2002     3095329 : END MODULE eri_mme_lattice_summation

Generated by: LCOV version 1.15