LCOV - code coverage report
Current view: top level - src/eri_mme - eri_mme_lattice_summation.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 79.7 % 681 543
Test Date: 2025-07-25 12:55:17 Functions: 54.1 % 281 152

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \brief 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       948038 :    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       948038 :       l_max = la_max + lb_max
      83       948038 :       alpha_G = a_mm + 0.25_dp/zeta + 0.25_dp/zetb
      84       948038 :       alpha_R = 0.25_dp/alpha_G
      85              : 
      86       948038 :       G_res = 0.5_dp*G_min
      87       948038 :       R_res = 0.5_dp*R_min
      88              : 
      89       948038 :       IF (PRESENT(G_rad)) G_rad = exp_radius(l_max, alpha_G, sum_precision, 1.0_dp, epsabs=G_res)
      90       948038 :       IF (PRESENT(R_rad)) R_rad = exp_radius(l_max, alpha_R, sum_precision, 1.0_dp, epsabs=R_res)
      91              : 
      92       948038 :    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      2411094 :    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      2411094 :       alpha = 0.25_dp/zeta
     123      2411094 :       beta = 0.25_dp/zetb
     124      2411094 :       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      2411094 :       IF (PRESENT(G_rads_1)) THEN
     129      2411094 :          G_res = 0.5_dp*G_min
     130      2411094 :          G_rads_1(1) = exp_radius(la_max, alpha, sum_precision, 1.0_dp, G_res)
     131      2411094 :          G_rads_1(2) = exp_radius(lb_max, beta, sum_precision, 1.0_dp, G_res)
     132      2411094 :          G_rads_1(3) = exp_radius(lc_max, gamma, sum_precision, 1.0_dp, G_res)
     133              :       END IF
     134              : 
     135              :       ! sum method 2
     136      2411094 :       IF (PRESENT(R_rads_2)) THEN
     137      2256540 :          R_res = 0.5_dp*R_min
     138      2256540 :          R_rads_2(1) = exp_radius(lb_max + lc_max, 0.25_dp/(beta + gamma), sum_precision, 1.0_dp, epsabs=R_res)
     139      2256540 :          R_rads_2(2) = exp_radius(lc_max + la_max, 0.25_dp/(alpha + gamma), sum_precision, 1.0_dp, epsabs=R_res)
     140      2256540 :          R_rads_2(3) = exp_radius(lb_max + la_max, 0.25_dp/(alpha + beta), sum_precision, 1.0_dp, epsabs=R_res)
     141              :       END IF
     142              : 
     143              :       ! sum method 3
     144              : 
     145      2411094 :       IF (PRESENT(R_rads_3)) THEN
     146      2256540 :          R_res = 0.5_dp*R_min
     147      2256540 :          alpha_R = 1.0_dp/((zeta + zetb + zetc)/((zeta + zetb)*zetc) + 4.0_dp*a_mm)
     148      2256540 :          R_rads_3(1) = exp_radius(la_max + lb_max, zeta*zetb/(zeta + zetb), sum_precision, 1.0_dp, R_res)
     149      2256540 :          R_rads_3(2) = exp_radius(la_max + lb_max + lc_max, alpha_R, sum_precision, 1.0_dp, R_res)
     150              :       END IF
     151              : 
     152      2411094 :    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       892104 :    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       892104 :       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     11597352 :       G_bounds = ellipsoid_bounds(G_rad, TRANSPOSE(hmat)/(2.0_dp*pi))
     197       892104 :       R_bounds = ellipsoid_bounds(R_rad, h_inv)
     198              : 
     199       892104 :       vol_G = twopi**3/vol
     200              : 
     201       892104 :       IF (is_ortho) THEN
     202      2981984 :          DO i_xyz = 1, 3
     203      2236488 :             ns_G = 2.0_dp*G_bounds(i_xyz)
     204      2236488 :             ns_R = 2.0_dp*R_bounds(i_xyz)
     205      2236488 :             n_sum_1d(1, i_xyz) = nsum_2c_gspace_1d(ns_G, la_max, lb_max)
     206      2981984 :             n_sum_1d(2, i_xyz) = nsum_2c_rspace_1d(ns_R, la_max, lb_max)
     207              :          END DO
     208              :       ELSE
     209       146608 :          ns_G = 4._dp/3._dp*pi*G_rad**3/vol_G
     210       146608 :          ns_R = 4._dp/3._dp*pi*R_rad**3/vol
     211       146608 :          n_sum_3d(1) = nsum_2c_gspace_3d(ns_G, la_max, lb_max)
     212       146608 :          n_sum_3d(2) = nsum_2c_rspace_3d(ns_R, la_max, lb_max)
     213              :       END IF
     214              : 
     215       892104 :    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      2256540 :    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      2256540 :                                G_rads_1=G_rads_1, R_rads_2=R_rads_2, R_rads_3=R_rads_3)
     265              : 
     266      2256540 :       vol_G = twopi**3/vol
     267              : 
     268      9026160 :       order_1 = MAXLOC(G_rads_1, DIM=1)
     269      9026160 :       order_2 = MINLOC(G_rads_1, DIM=1)
     270              : 
     271      9026160 :       DO i = 1, 3
     272    110570460 :          G_bounds_1(i, :) = ellipsoid_bounds(G_rads_1(i), TRANSPOSE(hmat)/(2.0_dp*pi))
     273              :       END DO
     274              : 
     275      9026160 :       DO i = 1, 3
     276     29335020 :          R_bounds_2(i, :) = ellipsoid_bounds(R_rads_2(i), h_inv)
     277              :       END DO
     278              : 
     279      6769620 :       DO i = 1, 2
     280     20308860 :          R_bounds_3(i, :) = ellipsoid_bounds(R_rads_3(i), h_inv)
     281              :       END DO
     282              : 
     283      2256540 :       IF (is_ortho) THEN
     284      8432160 :          DO i_xyz = 1, 3
     285              : 
     286      6324120 :             ns3_R1 = 2.0_dp*R_bounds_3(1, i_xyz)
     287      6324120 :             ns3_R2 = 2.0_dp*R_bounds_3(2, i_xyz)
     288              : 
     289      8432160 :             n_sum_1d(3, i_xyz) = nsum_3c_rspace_1d(ns3_R1, ns3_R2)
     290              :          END DO
     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              :       END IF
     328              : 
     329      2256540 :    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      2236488 :    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      2236488 :       nsum_2c_gspace_1d = NINT(ns_G*(2*exp_w + (l + m + 1)*5), KIND=int_8)
     344      2236488 :    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       202238 :    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       202238 :       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       202238 :       dG = inv_lgth*twopi
     369       202238 :       l_max = UBOUND(S_G, 1)
     370              : 
     371       606714 :       ALLOCATE (S_G_c(0:l_max))
     372      1025491 :       S_G_c(:) = 0.0_dp
     373      1079898 :       DO gg = -FLOOR(G_c), FLOOR(G_c)
     374       877660 :          G = gg*dG
     375       877660 :          exp_tot = EXP(-alpha*G**2)*EXP(gaussi*G*R) ! cost: 2 exp_w flops
     376       877660 :          G_pow_l = 1.0_dp
     377      4490891 :          DO l = 0, l_max
     378      3410993 :             S_G_c(l) = S_G_c(l) + G_pow_l*(-1.0_dp)**l*exp_tot ! cost: 4 flops
     379      4288653 :             G_pow_l = G_pow_l*G ! cost: 1 flop
     380              :          END DO
     381              :       END DO
     382              : 
     383      1848744 :       S_G(:) = REAL(S_G_c(0:l_max)*i_pow((/(l, l=0, l_max)/)))*inv_lgth
     384       202238 :    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       146608 :    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       146608 :       nsum_2c_gspace_3d = NINT(ns_G*(2*exp_w + ncoset(l + m)*7), KIND=int_8)
     399       146608 :    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        27581 :    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        55162 :       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        27581 :       REAL(KIND=dp), DIMENSION(3, 0:l_max)               :: G_pow_l
     432              :       REAL(KIND=dp), DIMENSION(3, 3)                     :: ht
     433              : 
     434       358553 :       ht = twopi*TRANSPOSE(h_inv)
     435       660268 :       Ig(:) = 0.0_dp
     436              : 
     437        27581 :       G_rads_sq = G_rad**2
     438              : 
     439       231084 :       DO gx = -FLOOR(G_c(1)), FLOOR(G_c(1))
     440       814012 :          G_x = ht(:, 1)*gx
     441      1913121 :          DO gy = -FLOOR(G_c(2)), FLOOR(G_c(2))
     442      6728148 :             G_y = ht(:, 2)*gy
     443     18506843 :             DO gz = -FLOOR(G_c(3)), FLOOR(G_c(3))
     444     66485212 :                G_z = ht(:, 3)*gz
     445              : 
     446     66485212 :                G = G_x + G_y + G_z
     447     16621303 :                G_sq = G(1)**2 + G(2)**2 + G(3)**2
     448     16621303 :                IF (G_sq .GT. G_rads_sq) CYCLE
     449              : 
     450      6396767 :                IF (PRESENT(potential)) THEN
     451      1033920 :                   IF (gx .EQ. 0 .AND. gy .EQ. 0 .AND. gz .EQ. 0) CYCLE
     452              :                END IF
     453              : 
     454     25586300 :                exp_tot = EXP(-alpha*G_sq)*EXP(gaussi*DOT_PRODUCT(G, R)) ! cost: 2 exp_w flops
     455      6396575 :                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              :                END IF
     466     25586300 :                DO k = 1, 3
     467     19189725 :                   G_pow_l(k, 0) = 1.0_dp
     468    110227076 :                   DO l = 1, l_max
     469    103830501 :                      G_pow_l(k, l) = G_pow_l(k, l - 1)*G(k)
     470              :                   END DO
     471              :                END DO
     472    466764837 :                DO lco = 1, ncoset(l_max)
     473    458686225 :                   CALL get_l(lco, l, lx, ly, lz)
     474   1834744900 :                   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    475307528 :                                           exp_tot*(-1.0_dp)**l*i_pow(l)
     478              :                END DO
     479              :             END DO
     480              :          END DO
     481              :       END DO
     482       660268 :       S_G(:) = REAL(Ig(:), KIND=dp)/vol
     483        27581 :    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      2236488 :    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      2236488 :       nsum_2c_rspace_1d = NINT(ns_R*(exp_w + (l + m + 1)*3), KIND=int_8)
     498      2236488 :    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      2212684 :    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      2212684 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: h_to_c
     521              : 
     522      2212684 :       dR = lgth
     523      2212684 :       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      9453605 :       S_R(:) = 0.0_dp
     527      2212684 :       Rp = R - R_c*dR
     528      9268511 :       DO rr = CEILING(-R_c - R/dR), FLOOR(R_c - R/dR)
     529      7055827 :          Rp = R + rr*dR
     530      7055827 :          exp_tot = EXP(-alpha*Rp**2) ! cost: exp_w flops
     531      7055827 :          R_pow_l = 1.0_dp
     532     32695273 :          DO l = 0, l_max
     533     23426762 :             S_R(l) = S_R(l) + R_pow_l*exp_tot ! cost: 2 flops
     534     30482589 :             R_pow_l = R_pow_l*Rp ! cost: 1 flop
     535              :          END DO
     536              :       END DO
     537              : 
     538              :       ! 2) C --> H
     539              :       CALL create_hermite_to_cartesian(alpha, l_max, h_to_c)
     540      9453605 :       S_R = MATMUL(TRANSPOSE(h_to_c(0:l_max, 0:l_max)), S_R)*SQRT(alpha/pi)
     541      4425368 :    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       146608 :    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       146608 :       nsum_2c_rspace_3d = NINT(ns_R*(exp_w + ncoset(l + m)*(4 + ncoset(l + m)*4)), KIND=int_8)
     556       146608 :    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       128219 :    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       128219 :       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       256438 :       REAL(KIND=dp), DIMENSION(3, 0:l_max)               :: R_pow_l
     585       128219 :       REAL(KIND=dp), DIMENSION(ncoset(l_max))            :: S_R_C
     586              : 
     587      2057776 :       S_R(:) = 0.0_dp
     588      2057776 :       S_R_C(:) = 0.0_dp
     589              : 
     590      2051504 :       s_shift = MATMUL(h_inv, -R)
     591       512876 :       R_l = -R_c + s_shift
     592       512876 :       R_r = R_c + s_shift
     593              : 
     594       128219 :       R_rad_sq = R_rad**2
     595              : 
     596       432234 :       DO sx = CEILING(R_l(1)), FLOOR(R_r(1))
     597      1216060 :          Rx = hmat(:, 1)*sx
     598      1470931 :          DO sy = CEILING(R_l(2)), FLOOR(R_r(2))
     599      4154788 :             Ry = hmat(:, 2)*sy
     600      5870059 :             DO sz = CEILING(R_l(3)), FLOOR(R_r(3))
     601     18109388 :                Rz = hmat(:, 3)*sz
     602     18109388 :                Rp = Rx + Ry + Rz
     603      4527347 :                R_sq = (Rp(1) + R(1))**2 + (Rp(2) + R(2))**2 + (Rp(3) + R(3))**2
     604      4527347 :                IF (R_sq .GT. R_rad_sq) CYCLE
     605      1814156 :                exp_tot = EXP(-alpha*R_sq) ! cost: exp_w flops
     606      7256624 :                DO k = 1, 3
     607      5442468 :                   R_pow_l(k, 0) = 1.0_dp
     608     16096739 :                   DO l = 1, l_max
     609     14282583 :                      R_pow_l(k, l) = R_pow_l(k, l - 1)*(Rp(k) + R(k))
     610              :                   END DO
     611              :                END DO
     612     22199618 :                DO lco = 1, ncoset(l_max)
     613     19346765 :                   CALL get_l(lco, l, lx, ly, lz)
     614     23874112 :                   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              :                END DO
     616              :             END DO
     617              :          END DO
     618              :       END DO
     619              : 
     620              :       CALL create_hermite_to_cartesian(alpha, l_max, h_to_c)
     621              : 
     622      2057776 :       DO lco = 1, ncoset(l_max)
     623      1929557 :          CALL get_l(lco, l, lx, ly, lz)
     624              : 
     625      5692372 :          DO llx = 0, lx
     626     11963703 :          DO lly = 0, ly
     627     20729884 :          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     17095288 :                                      S_R_C(coset(llx, lly, llz))
     631              :          END DO
     632              :          END DO
     633              :          END DO
     634              :       END DO
     635      2057776 :       S_R(:) = S_R(:)*(alpha/pi)**1.5_dp
     636              : 
     637       256438 :    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       295044 :    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       295044 :       prefac = prefactor*inv_lgth
     666       295044 :       dG = inv_lgth*twopi ! positive
     667       295044 :       l_max = UBOUND(S_G, 1)
     668              : 
     669      1660104 :       S_G(:) = 0.0_dp
     670       295044 :       IF (0 .LT. G_min) THEN
     671              :          k0 = G_min; k1 = 0
     672        73761 :       ELSE IF (G_c .LT. 0) THEN
     673              :          k0 = 0; k1 = G_c
     674              :       ELSE ! Gmin <= 0 /\ 0 <= Gc
     675        73761 :          S_G(0) = prefac
     676              :          k0 = 1; k1 = -1
     677              :       END IF
     678              :       ! positive range
     679              :       IF (0 .LT. k0) THEN
     680     14185303 :          DO k = k0, G_c
     681     13890259 :             G = k*dG; exp_tot = EXP(-alpha*G**2)*prefac
     682     83043248 :             DO l = 0, l_max
     683     82748204 :                S_G(l) = S_G(l) + G**(l - delta_l)*exp_tot
     684              :             END DO
     685              :          END DO
     686              :       END IF
     687              :       ! negative range
     688       295044 :       IF (k1 .LT. 0) THEN
     689      4816941 :          DO k = G_min, k1
     690      4743180 :             G = k*dG; exp_tot = EXP(-alpha*G**2)*prefac
     691     27823869 :             DO l = 0, l_max
     692     27750108 :                S_G(l) = S_G(l) + ABS(G)**(l - delta_l)*exp_tot
     693              :             END DO
     694              :          END DO
     695              :       END IF
     696       295044 :    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      6324120 :    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      6324120 :       l_max = UBOUND(S_G, 1)
     746      6324120 :       m_max = UBOUND(S_G, 2)
     747      6324120 :       n_max = UBOUND(S_G, 3)
     748              : 
     749      6324120 :       nR1 = 2*FLOOR(R_bounds_3(1)) + 1
     750      6324120 :       nR2 = 2*FLOOR(R_bounds_3(2)) + 1
     751              : 
     752      6324120 :       IF (nR1*nR2 > 1 + nR1*2) THEN
     753              :          prop_exp = 1
     754              :       ELSE
     755      2997690 :          prop_exp = 0
     756              :       END IF
     757              : 
     758     25296480 :       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      6231960 :             IF (l_max == ${l_max}$) THEN
     763              :                #:for m_max in range(0,lmax_unroll+1)
     764      6231960 :                   IF (m_max == ${m_max}$) THEN
     765              :                      #:for n_max in range(0, lmax_unroll+1)
     766      6231960 :                         IF (n_max == ${n_max}$) THEN
     767              :                            #:for prop_exp in range(0,2)
     768      6231960 :                               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      6231960 :                                                                                            zeta, zetb, zetc, a_mm, lgth, R_bounds_3)
     771      6231960 :                                  RETURN
     772              :                               END IF
     773              :                            #:endfor
     774              :                         END IF
     775              :                      #:endfor
     776              :                   END IF
     777              :                #:endfor
     778              :             END IF
     779              :          #:endfor
     780              :       END IF
     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      6324120 :    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      6324120 :       nsum_3c_rspace_1d = NINT(MIN((4 + ns_R1*2), ns_R1*(ns_R2 + 1)), KIND=int_8)
     829      6324120 :    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              :          END DO
     880              :       END DO
     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              :             END DO
     897              :          END DO
     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              :                END DO
     924              :             END DO
     925              :          END DO
     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              :             END DO
     934              :          END DO
     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              :             END DO
     943              :          END DO
     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              :                   END DO
     951              :                END DO
     952              :             END DO
     953              :          END DO
     954              :       END DO
     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      6231960 :                PURE SUBROUTINE pgf_sum_3c_rspace_1d_${l_max}$_${m_max}$_${n_max}$_exp_${prop_exp}$ ( &
     980      6231960 :                   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      6231960 :                   dR = lgth
    1001      6231960 :                   alpha = 1.0_dp/((zeta + zetb + zetc)/((zeta + zetb)*zetc) + 4.0_dp*a_mm)
    1002              : 
    1003    113325552 :                   S_R(:, :, :) = 0.0_dp
    1004              : 
    1005      6231960 :                   R_offset = RC - (zeta*RA + zetb*RB)/(zeta + zetb)
    1006              : 
    1007      6231960 :                   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      5509656 :                            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      5509656 :                            h_to_c_${k}$_${l+1}$ = 0.0_dp
    1015              :                         #:endif
    1016              :                      #:endfor
    1017              :                   #:endfor
    1018              : 
    1019              :                   #:if prop_exp
    1020      3268062 :                      exp_dRsq = exp(-alpha*dR*dR)
    1021      3268062 :                      exp_2dRsq = exp_dRsq*exp_dRsq
    1022              :                   #:endif
    1023              : 
    1024      6231960 :                   rr1_delta = (RA - RB)/dR
    1025              : 
    1026      6231960 :                   rr1_l = CEILING(-R_c(1) + rr1_delta)
    1027      6231960 :                   rr1_r = FLOOR(R_c(1) + rr1_delta)
    1028              : 
    1029      6231960 :                   R1 = rr1_l*dR
    1030              : 
    1031      6231960 :                   alpha_E = zeta*zetb/(zeta + zetb)
    1032              : 
    1033     20231352 :                   DO rr1 = rr1_l, rr1_r
    1034              :                      #:for t in range(0, l_tot_max + 1)
    1035     13999392 :                         S_R_t_${t}$ = 0.0_dp
    1036     13999392 :                         S_R_t2_${t}$ = 0.0_dp
    1037              :                      #:endfor
    1038     13999392 :                      R_tmp = R_offset + R1*zeta/(zeta + zetb)
    1039     13999392 :                      rr2_delta = -R_tmp/dR
    1040              : 
    1041     13999392 :                      rr2_l = CEILING(-R_c(2) + rr2_delta)
    1042     13999392 :                      rr2_r = FLOOR(R_c(2) + rr2_delta)
    1043              : 
    1044     13999392 :                      R = R_tmp + (rr2_l)*dR
    1045              : 
    1046              :                      #:if prop_exp
    1047      9119092 :                         exp2_2RdR = exp(-2*alpha*R*dR)
    1048      9119092 :                         exp2_Rsq = exp(-alpha*R*R)
    1049              :                      #:endif
    1050              : 
    1051     71598485 :                      DO rr2 = rr2_l, rr2_r
    1052     57599093 :                         R_pow_t = 1.0_dp
    1053              :                         #:if not prop_exp
    1054      7179779 :                            exp2_Rsq = exp(-alpha*R*R)
    1055              :                         #:endif
    1056              :                         #:for t in range(0, l_tot_max + 1)
    1057     57599093 :                            S_R_t_${t}$ = S_R_t_${t}$+R_pow_t*exp2_Rsq
    1058              :                            #:if t<l_tot_max
    1059     51468255 :                               R_pow_t = R_pow_t*R
    1060              :                            #:endif
    1061              :                         #:endfor
    1062              : 
    1063              :                         #:if prop_exp
    1064     50419314 :                            exp2_Rsq = exp2_Rsq*exp_dRsq*exp2_2RdR
    1065     50419314 :                            exp2_2RdR = exp2_2RdR*exp_2dRsq
    1066              :                         #:endif
    1067     71598485 :                         R = R + dR
    1068              :                      END DO
    1069              : 
    1070              :                      ! C --> H
    1071              :                      #:for l in range(0, l_tot_max+1)
    1072              :                         #:for k in range(0, l+1)
    1073     13999392 :                            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     13999392 :                      E_0_0_0 = exp(-alpha_E*(RA - RB - R1)*(RA - RB - R1))
    1079              : 
    1080              :                      #:if l_tot_max > 0
    1081     12382032 :                         c1 = 1.0_dp/(zeta + zetb)
    1082     12382032 :                         c2 = 2.0_dp*(zetb/(zeta + zetb))*(RB - (RA - R1))
    1083     12382032 :                         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      8760456 :                                  #{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      8761656 :                                  #{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     13999392 :               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     20231352 :                      R1 = R1 + dR
    1115              :                   END DO
    1116              : 
    1117    113325552 :                   S_R = S_R*pi**(-0.5_dp)*((zeta + zetb)/(zeta*zetb))**(-0.5_dp)
    1118      6231960 :                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              :       END IF
    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              :          END IF
    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              :          END IF
    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        29230 :             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              :       END IF
    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              :                            END IF
    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              :                         END DO
    1356              :                      END DO
    1357              :                   END DO
    1358              :                END DO
    1359              :             END DO
    1360              :          END DO
    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              :                            END IF
    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              :                         END DO
    1389              :                      END DO
    1390              :                   END DO
    1391              :                END DO
    1392              :             END DO
    1393              :          END DO
    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              :                            END IF
    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              :                         END DO
    1422              :                      END DO
    1423              :                   END DO
    1424              :                END DO
    1425              :             END DO
    1426              :          END DO
    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              :          END DO
    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              :          END DO
    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              :          END DO
    1495              :       END DO
    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              :       END DO
    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              :       END DO
    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              :       END DO
    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              :             END DO
    1519              :          END DO
    1520              :       END DO
    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              :                   END DO
    1616              :                END DO
    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              :                END DO
    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              :                      END DO
    1636              :                   END DO
    1637              :                END DO
    1638              :             END DO
    1639              :          END DO
    1640              :       END DO
    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              :                   END DO
    1706              :                END DO
    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              :                END DO
    1714              : 
    1715              :             END DO
    1716              :          END DO
    1717              :       END DO
    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              :             END DO
    1735              :             END DO
    1736              :             END DO
    1737              :          END DO
    1738              :       END DO
    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              :          END DO
    1835              :          END DO
    1836              :          END DO
    1837              :       END DO
    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              :                            END DO
    1879              :                         END DO
    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              :                         END DO
    1884              :                      END DO
    1885              :                   END DO
    1886              :                END DO
    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              :                      END DO
    1904              :                      END DO
    1905              :                      END DO
    1906              :                   END DO
    1907              :                END DO
    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              :                   END DO
    1919              :                END DO
    1920              :             END DO
    1921              :          END DO
    1922              :       END DO
    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     19854720 :    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     79418880 :       DO i_xyz = 1, 3
    1940     79418880 :          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              :       END DO
    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              :       END IF
    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    843147983 :    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   3372591932 :       l = SUM(indco(:, lco))
    1981    843147983 :       IF (PRESENT(lx)) lx = indco(1, lco)
    1982    843147983 :       IF (PRESENT(ly)) ly = indco(2, lco)
    1983    843147983 :       IF (PRESENT(lz)) lz = indco(3, lco)
    1984    843147983 :    END SUBROUTINE
    1985              : 
    1986              : ! **************************************************************************************************
    1987              : !> \brief ...
    1988              : !> \param i ...
    1989              : !> \return ...
    1990              : ! **************************************************************************************************
    1991    479950601 :    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    479950601 :       i_pow = ip(MOD(i, 4))
    1999              : 
    2000    479950601 :    END FUNCTION
    2001              : 
    2002      3955399 : END MODULE eri_mme_lattice_summation
        

Generated by: LCOV version 2.0-1