LCOV - code coverage report
Current view: top level - src/shg_int - generic_shg_integrals.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 100.0 % 511 511
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 13 13

            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 Calculation of contracted, spherical Gaussian integrals using the solid harmonic
      10              : !>        Gaussian (SHG) integral scheme. Routines for the following two-center integrals:
      11              : !>        i)  (a|O(r12)|b) where O(r12) is the overlap, coulomb operator etc.
      12              : !>        ii) (aba) and (abb) s-overlaps
      13              : !> \par Literature
      14              : !>      T.J. Giese and D. M. York, J. Chem. Phys, 128, 064104 (2008)
      15              : !>      T. Helgaker, P Joergensen, J. Olsen, Molecular Electronic-Structure
      16              : !>                                           Theory, Wiley
      17              : !> \par History
      18              : !>      created [05.2016]
      19              : !> \author Dorothea Golze
      20              : ! **************************************************************************************************
      21              : MODULE generic_shg_integrals
      22              :    USE basis_set_types,                 ONLY: gto_basis_set_type
      23              :    USE constants_operator,              ONLY: operator_coulomb,&
      24              :                                               operator_gauss,&
      25              :                                               operator_verf,&
      26              :                                               operator_verfc,&
      27              :                                               operator_vgauss
      28              :    USE construct_shg,                   ONLY: &
      29              :         construct_dev_shg_ab, construct_int_shg_ab, construct_overlap_shg_aba, &
      30              :         construct_overlap_shg_abb, dev_overlap_shg_aba, dev_overlap_shg_abb, get_W_matrix, &
      31              :         get_dW_matrix, get_real_scaled_solid_harmonic
      32              :    USE kinds,                           ONLY: dp
      33              :    USE orbital_pointers,                ONLY: nsoset
      34              :    USE s_contract_shg,                  ONLY: &
      35              :         contract_s_overlap_aba, contract_s_overlap_abb, contract_s_ra2m_ab, &
      36              :         contract_sint_ab_chigh, contract_sint_ab_clow, s_coulomb_ab, s_gauss_ab, s_overlap_ab, &
      37              :         s_overlap_abb, s_ra2m_ab, s_verf_ab, s_verfc_ab, s_vgauss_ab
      38              : #include "../base/base_uses.f90"
      39              : 
      40              :    IMPLICIT NONE
      41              : 
      42              :    PRIVATE
      43              : 
      44              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'generic_shg_integrals'
      45              : 
      46              :    PUBLIC :: int_operators_r12_ab_shg, int_overlap_ab_shg, &
      47              :              int_ra2m_ab_shg, int_overlap_aba_shg, int_overlap_abb_shg, &
      48              :              get_abb_same_kind, lri_precalc_angular_shg_part, &
      49              :              int_overlap_ab_shg_low, int_ra2m_ab_shg_low, int_overlap_aba_shg_low, &
      50              :              int_overlap_abb_shg_low
      51              : 
      52              :    ABSTRACT INTERFACE
      53              : ! **************************************************************************************************
      54              : !> \brief Interface for the calculation of integrals over s-functions and their scalar derivatives
      55              : !>        with respect to rab2
      56              : !> \param la_max ...
      57              : !> \param npgfa ...
      58              : !> \param zeta ...
      59              : !> \param lb_max ...
      60              : !> \param npgfb ...
      61              : !> \param zetb ...
      62              : !> \param omega ...
      63              : !> \param rab ...
      64              : !> \param v matrix storing the integrals and scalar derivatives
      65              : !> \param calculate_forces ...
      66              : ! **************************************************************************************************
      67              :       SUBROUTINE ab_sint_shg(la_max, npgfa, zeta, lb_max, npgfb, zetb, omega, rab, v, calculate_forces)
      68              :          USE kinds, ONLY: dp
      69              :       INTEGER, INTENT(IN)                                :: la_max, npgfa
      70              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: zeta
      71              :       INTEGER, INTENT(IN)                                :: lb_max, npgfb
      72              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: zetb
      73              :       REAL(KIND=dp), INTENT(IN)                          :: omega
      74              :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rab
      75              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT)   :: v
      76              :       LOGICAL, INTENT(IN)                                :: calculate_forces
      77              : 
      78              :       END SUBROUTINE ab_sint_shg
      79              :    END INTERFACE
      80              : 
      81              : ! **************************************************************************************************
      82              : 
      83              : CONTAINS
      84              : 
      85              : ! **************************************************************************************************
      86              : !> \brief Calcululates the two-center integrals of the type (a|O(r12)|b) using the SHG scheme
      87              : !> \param r12_operator the integral operator, which depends on r12=|r1-r2|
      88              : !> \param vab integral matrix of spherical contracted Gaussian functions
      89              : !> \param dvab derivative of the integrals
      90              : !> \param rab distance vector between center A and B
      91              : !> \param fba basis at center A
      92              : !> \param fbb basis at center B
      93              : !> \param scona_shg SHG contraction matrix for A
      94              : !> \param sconb_shg SHG contraction matrix for B
      95              : !> \param omega parameter in the operator
      96              : !> \param calculate_forces ...
      97              : ! **************************************************************************************************
      98      2528004 :    SUBROUTINE int_operators_r12_ab_shg(r12_operator, vab, dvab, rab, fba, fbb, scona_shg, sconb_shg, &
      99              :                                        omega, calculate_forces)
     100              : 
     101              :       INTEGER, INTENT(IN)                                :: r12_operator
     102              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT)      :: vab
     103              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT), &
     104              :          OPTIONAL                                        :: dvab
     105              :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rab
     106              :       TYPE(gto_basis_set_type), POINTER                  :: fba, fbb
     107              :       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: scona_shg, sconb_shg
     108              :       REAL(KIND=dp), INTENT(IN), OPTIONAL                :: omega
     109              :       LOGICAL, INTENT(IN)                                :: calculate_forces
     110              : 
     111              :       INTEGER                                            :: la_max, lb_max
     112              :       REAL(KIND=dp)                                      :: my_omega
     113      2528004 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: Waux_mat
     114      2528004 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :)  :: dWaux_mat
     115              : 
     116              :       PROCEDURE(ab_sint_shg), POINTER                    :: s_operator_ab
     117              : 
     118      2528004 :       NULLIFY (s_operator_ab)
     119              : 
     120      7385868 :       la_max = MAXVAL(fba%lmax)
     121      7386668 :       lb_max = MAXVAL(fbb%lmax)
     122              : 
     123              :       CALL precalc_angular_shg_part(la_max, lb_max, rab, Waux_mat, dWaux_mat, calculate_forces)
     124      2528004 :       my_omega = 1.0_dp
     125              : 
     126      5055880 :       SELECT CASE (r12_operator)
     127              :       CASE (operator_coulomb)
     128      2527876 :          s_operator_ab => s_coulomb_ab
     129              :       CASE (operator_verf)
     130           32 :          s_operator_ab => s_verf_ab
     131           32 :          IF (PRESENT(omega)) my_omega = omega
     132              :       CASE (operator_verfc)
     133           32 :          s_operator_ab => s_verfc_ab
     134           32 :          IF (PRESENT(omega)) my_omega = omega
     135              :       CASE (operator_vgauss)
     136           32 :          s_operator_ab => s_vgauss_ab
     137           32 :          IF (PRESENT(omega)) my_omega = omega
     138              :       CASE (operator_gauss)
     139           32 :          s_operator_ab => s_gauss_ab
     140           32 :          IF (PRESENT(omega)) my_omega = omega
     141              :       CASE DEFAULT
     142      2528004 :          CPABORT("Operator not available")
     143              :       END SELECT
     144              : 
     145              :       CALL int_operator_ab_shg_low(s_operator_ab, vab, dvab, rab, fba, fbb, scona_shg, sconb_shg, &
     146      5055848 :                                    my_omega, Waux_mat, dWaux_mat, calculate_forces)
     147              : 
     148      2528004 :       DEALLOCATE (Waux_mat, dWaux_mat)
     149              : 
     150      2528004 :    END SUBROUTINE int_operators_r12_ab_shg
     151              : 
     152              : ! **************************************************************************************************
     153              : !> \brief calculate overlap integrals (a,b)
     154              : !> \param vab integral (a,b)
     155              : !> \param dvab derivative of sab
     156              : !> \param rab distance vector
     157              : !> \param fba basis at center A
     158              : !> \param fbb basis at center B
     159              : !> \param scona_shg contraction matrix A
     160              : !> \param sconb_shg contraxtion matrix B
     161              : !> \param calculate_forces ...
     162              : ! **************************************************************************************************
     163           32 :    SUBROUTINE int_overlap_ab_shg(vab, dvab, rab, fba, fbb, scona_shg, sconb_shg, &
     164              :                                  calculate_forces)
     165              : 
     166              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
     167              :          INTENT(INOUT)                                   :: vab
     168              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
     169              :          INTENT(INOUT)                                   :: dvab
     170              :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rab
     171              :       TYPE(gto_basis_set_type), POINTER                  :: fba, fbb
     172              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN)      :: scona_shg, sconb_shg
     173              :       LOGICAL, INTENT(IN)                                :: calculate_forces
     174              : 
     175              :       INTEGER                                            :: la_max, lb_max
     176           32 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: Waux_mat
     177           32 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :)  :: dWaux_mat
     178              : 
     179          352 :       la_max = MAXVAL(fba%lmax)
     180          512 :       lb_max = MAXVAL(fbb%lmax)
     181              : 
     182              :       CALL precalc_angular_shg_part(la_max, lb_max, rab, Waux_mat, dWaux_mat, calculate_forces)
     183              : 
     184              :       CALL int_overlap_ab_shg_low(vab, dvab, rab, fba, fbb, scona_shg, sconb_shg, &
     185           32 :                                   Waux_mat, dWaux_mat, .TRUE., calculate_forces, contraction_high=.TRUE.)
     186              : 
     187           32 :       DEALLOCATE (Waux_mat, dWaux_mat)
     188              : 
     189           32 :    END SUBROUTINE int_overlap_ab_shg
     190              : 
     191              : ! **************************************************************************************************
     192              : !> \brief Calcululates the two-center integrals of the type (a|(r-Ra)^(2m)|b) using the SHG scheme
     193              : !> \param vab integral matrix of spherical contracted Gaussian functions
     194              : !> \param dvab derivative of the integrals
     195              : !> \param rab distance vector between center A and B
     196              : !> \param fba basis at center A
     197              : !> \param fbb basis at center B
     198              : !> \param scon_ra2m contraction matrix for A including the combinatorial factors
     199              : !> \param sconb_shg SHG contraction matrix for B
     200              : !> \param m exponent in (r-Ra)^(2m) operator
     201              : !> \param calculate_forces ...
     202              : ! **************************************************************************************************
     203           32 :    SUBROUTINE int_ra2m_ab_shg(vab, dvab, rab, fba, fbb, scon_ra2m, sconb_shg, &
     204              :                               m, calculate_forces)
     205              : 
     206              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT)      :: vab
     207              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT)   :: dvab
     208              :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rab
     209              :       TYPE(gto_basis_set_type), POINTER                  :: fba, fbb
     210              :       REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN)   :: scon_ra2m
     211              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN)      :: sconb_shg
     212              :       INTEGER, INTENT(IN)                                :: m
     213              :       LOGICAL, INTENT(IN)                                :: calculate_forces
     214              : 
     215              :       INTEGER                                            :: la_max, lb_max
     216           32 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: Waux_mat
     217           32 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :)  :: dWaux_mat
     218              : 
     219          352 :       la_max = MAXVAL(fba%lmax)
     220          512 :       lb_max = MAXVAL(fbb%lmax)
     221              : 
     222              :       CALL precalc_angular_shg_part(la_max, lb_max, rab, Waux_mat, dWaux_mat, calculate_forces)
     223              :       CALL int_ra2m_ab_shg_low(vab, dvab, rab, fba, fbb, sconb_shg, scon_ra2m, &
     224           32 :                                m, Waux_mat, dWaux_mat, calculate_forces)
     225              : 
     226           32 :       DEALLOCATE (Waux_mat, dWaux_mat)
     227              : 
     228           32 :    END SUBROUTINE int_ra2m_ab_shg
     229              : 
     230              : ! **************************************************************************************************
     231              : !> \brief calculate integrals (a,b,fa)
     232              : !> \param saba integral [aba]
     233              : !> \param dsaba derivative of [aba]
     234              : !> \param rab distance vector between A and B
     235              : !> \param oba orbital basis at center A
     236              : !> \param obb orbital basis at center B
     237              : !> \param fba auxiliary basis set at center A
     238              : !> \param scon_obb contraction matrix for orb bas on B
     239              : !> \param scona_mix mixed contraction matrix orb + ri basis on A
     240              : !> \param oba_index orbital basis index for scona_mix
     241              : !> \param fba_index ri basis index for scona_mix
     242              : !> \param cg_coeff Clebsch-Gordon coefficients
     243              : !> \param cg_none0_list list of none-zero Clebsch-Gordon coefficients
     244              : !> \param ncg_none0 number of non-zero Clebsch-Gordon coefficients
     245              : !> \param calculate_forces ...
     246              : ! **************************************************************************************************
     247          144 :    SUBROUTINE int_overlap_aba_shg(saba, dsaba, rab, oba, obb, fba, scon_obb, &
     248          144 :                                   scona_mix, oba_index, fba_index, &
     249          144 :                                   cg_coeff, cg_none0_list, ncg_none0, &
     250              :                                   calculate_forces)
     251              : 
     252              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
     253              :          INTENT(INOUT)                                   :: saba
     254              :       REAL(KIND=dp), ALLOCATABLE, &
     255              :          DIMENSION(:, :, :, :), INTENT(INOUT)            :: dsaba
     256              :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rab
     257              :       TYPE(gto_basis_set_type), POINTER                  :: oba, obb, fba
     258              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN)      :: scon_obb
     259              :       REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN)   :: scona_mix
     260              :       INTEGER, DIMENSION(:, :, :), INTENT(IN)            :: oba_index, fba_index
     261              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN)      :: cg_coeff
     262              :       INTEGER, DIMENSION(:, :, :), INTENT(IN)            :: cg_none0_list
     263              :       INTEGER, DIMENSION(:, :), INTENT(IN)               :: ncg_none0
     264              :       LOGICAL, INTENT(IN)                                :: calculate_forces
     265              : 
     266              :       INTEGER                                            :: laa_max, lb_max
     267          144 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: Waux_mat
     268          144 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :)  :: dWaux_mat
     269              : 
     270         1948 :       laa_max = MAXVAL(oba%lmax) + MAXVAL(fba%lmax)
     271          332 :       lb_max = MAXVAL(obb%lmax)
     272              : 
     273      1884366 :       saba = 0.0_dp
     274       951984 :       IF (calculate_forces) dsaba = 0.0_dp
     275              :       CALL precalc_angular_shg_part(laa_max, lb_max, rab, Waux_mat, dWaux_mat, &
     276              :                                     calculate_forces)
     277              :       CALL int_overlap_aba_shg_low(saba, dsaba, rab, oba, obb, fba, &
     278              :                                    scon_obb, scona_mix, oba_index, fba_index, &
     279              :                                    cg_coeff, cg_none0_list, ncg_none0, &
     280          144 :                                    Waux_mat, dWaux_mat, .TRUE., calculate_forces)
     281              : 
     282          144 :       DEALLOCATE (Waux_mat, dWaux_mat)
     283              : 
     284          144 :    END SUBROUTINE int_overlap_aba_shg
     285              : 
     286              : ! **************************************************************************************************
     287              : !> \brief calculate integrals (a,b,fb)
     288              : !> \param sabb integral [abb]
     289              : !> \param dsabb derivative of [abb]
     290              : !> \param rab distance vector between A and B
     291              : !> \param oba orbital basis at center A
     292              : !> \param obb orbital basis at center B
     293              : !> \param fbb auxiliary basis set at center B
     294              : !> \param scon_oba contraction matrix for orb bas on A
     295              : !> \param sconb_mix mixed contraction matrix orb + ri basis on B
     296              : !> \param obb_index orbital basis index for sconb_mix
     297              : !> \param fbb_index ri basis index for sconb_mix
     298              : !> \param cg_coeff Clebsch-Gordon coefficients
     299              : !> \param cg_none0_list list of none-zero Clebsch-Gordon coefficients
     300              : !> \param ncg_none0 number of non-zero Clebsch-Gordon coefficients
     301              : !> \param calculate_forces ...
     302              : ! **************************************************************************************************
     303           16 :    SUBROUTINE int_overlap_abb_shg(sabb, dsabb, rab, oba, obb, fbb, scon_oba, &
     304           16 :                                   sconb_mix, obb_index, fbb_index, &
     305           16 :                                   cg_coeff, cg_none0_list, ncg_none0, &
     306              :                                   calculate_forces)
     307              : 
     308              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
     309              :          INTENT(INOUT)                                   :: sabb
     310              :       REAL(KIND=dp), ALLOCATABLE, &
     311              :          DIMENSION(:, :, :, :), INTENT(INOUT)            :: dsabb
     312              :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rab
     313              :       TYPE(gto_basis_set_type), POINTER                  :: oba, obb, fbb
     314              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN)      :: scon_oba
     315              :       REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN)   :: sconb_mix
     316              :       INTEGER, DIMENSION(:, :, :), INTENT(IN)            :: obb_index, fbb_index
     317              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN)      :: cg_coeff
     318              :       INTEGER, DIMENSION(:, :, :), INTENT(IN)            :: cg_none0_list
     319              :       INTEGER, DIMENSION(:, :), INTENT(IN)               :: ncg_none0
     320              :       LOGICAL, INTENT(IN)                                :: calculate_forces
     321              : 
     322              :       INTEGER                                            :: la_max, lbb_max
     323           16 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: Waux_mat
     324           16 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :)  :: dWaux_mat
     325              : 
     326           32 :       la_max = MAXVAL(oba%lmax)
     327          272 :       lbb_max = MAXVAL(obb%lmax) + MAXVAL(fbb%lmax)
     328              : 
     329       317280 :       sabb = 0.0_dp
     330       951856 :       IF (calculate_forces) dsabb = 0.0_dp
     331              :       CALL precalc_angular_shg_part(lbb_max, la_max, rab, Waux_mat, dWaux_mat, &
     332              :                                     calculate_forces)
     333              :       CALL int_overlap_abb_shg_low(sabb, dsabb, rab, oba, obb, fbb, &
     334              :                                    scon_oba, sconb_mix, obb_index, fbb_index, &
     335              :                                    cg_coeff, cg_none0_list, ncg_none0, &
     336           16 :                                    Waux_mat, dWaux_mat, .TRUE., calculate_forces)
     337              : 
     338           16 :       DEALLOCATE (Waux_mat, dWaux_mat)
     339              : 
     340           16 :    END SUBROUTINE int_overlap_abb_shg
     341              : 
     342              : ! **************************************************************************************************
     343              : !> \brief precalculates the angular part of the SHG integrals
     344              : !> \param la_max ...
     345              : !> \param lb_max ...
     346              : !> \param rab distance vector between a and b
     347              : !> \param Waux_mat W matrix that contains the angular-dependent part
     348              : !> \param dWaux_mat derivative of the W matrix
     349              : !> \param calculate_forces ...
     350              : ! **************************************************************************************************
     351      2528228 :    SUBROUTINE precalc_angular_shg_part(la_max, lb_max, rab, Waux_mat, dWaux_mat, calculate_forces)
     352              : 
     353              :       INTEGER, INTENT(IN)                                :: la_max, lb_max
     354              :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rab
     355              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
     356              :          INTENT(OUT)                                     :: Waux_mat
     357              :       REAL(KIND=dp), ALLOCATABLE, &
     358              :          DIMENSION(:, :, :, :), INTENT(OUT)              :: dWaux_mat
     359              :       LOGICAL, INTENT(IN)                                :: calculate_forces
     360              : 
     361              :       INTEGER                                            :: lmax, mdim(3)
     362              :       INTEGER, DIMENSION(:), POINTER                     :: la_max_all
     363              :       REAL(KIND=dp)                                      :: rab2
     364      2528228 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: Rc, Rs
     365              : 
     366      2528228 :       NULLIFY (la_max_all)
     367      2528228 :       rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
     368      2528228 :       lmax = MAX(la_max, lb_max)
     369              : 
     370      7584684 :       ALLOCATE (la_max_all(0:lb_max))
     371     15169368 :       ALLOCATE (Rc(0:lmax, -2*lmax:2*lmax), Rs(0:lmax, -2*lmax:2*lmax))
     372     37555622 :       Rc = 0._dp
     373     37555622 :       Rs = 0._dp
     374      2528228 :       mdim(1) = MIN(la_max, lb_max) + 1
     375      2528228 :       mdim(2) = nsoset(la_max) + 1
     376      2528228 :       mdim(3) = nsoset(lb_max) + 1
     377     12641140 :       ALLOCATE (Waux_mat(mdim(1), mdim(2), mdim(3)))
     378     12641140 :       ALLOCATE (dWaux_mat(3, mdim(1), mdim(2), mdim(3)))
     379              : 
     380      6391782 :       la_max_all(0:lb_max) = la_max
     381              :       !*** -rab, since Eq. in Ref. use Ra-Rb, not Rb-Ra
     382     10112912 :       CALL get_real_scaled_solid_harmonic(Rc, Rs, lmax, -rab, rab2)
     383      2528228 :       CALL get_W_matrix(la_max_all, lb_max, lmax, Rc, Rs, Waux_mat)
     384      2528228 :       IF (calculate_forces) THEN
     385          256 :          CALL get_dW_matrix(la_max_all, lb_max, Waux_mat, dWaux_mat)
     386              :       END IF
     387              : 
     388      2528228 :       DEALLOCATE (Rc, Rs, la_max_all)
     389              : 
     390      2528228 :    END SUBROUTINE precalc_angular_shg_part
     391              : 
     392              : ! **************************************************************************************************
     393              : !> \brief calculate integrals (a|O(r12)|b)
     394              : !> \param s_operator_ab procedure pointer for the respective operator. The integral evaluation
     395              : !>        differs only in the calculation of the [s|O(r12)|s] integrals and their scalar
     396              : !>        derivatives.
     397              : !> \param vab integral matrix of spherical contracted Gaussian functions
     398              : !> \param dvab derivative of the integrals
     399              : !> \param rab distance vector between center A and B
     400              : !> \param fba basis at center A
     401              : !> \param fbb basis at center B
     402              : !> \param scona_shg SHG contraction matrix for A
     403              : !> \param sconb_shg SHG contraction matrix for B
     404              : !> \param omega parameter in the operator
     405              : !> \param Waux_mat W matrix that contains the angular-dependent part
     406              : !> \param dWaux_mat derivative of the W matrix
     407              : !> \param calculate_forces ...
     408              : ! **************************************************************************************************
     409      2528004 :    SUBROUTINE int_operator_ab_shg_low(s_operator_ab, vab, dvab, rab, fba, fbb, scona_shg, sconb_shg, &
     410      2528004 :                                       omega, Waux_mat, dWaux_mat, calculate_forces)
     411              : 
     412              :       PROCEDURE(ab_sint_shg), POINTER                    :: s_operator_ab
     413              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT)      :: vab
     414              :       REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, INTENT(INOUT)   :: dvab
     415              :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rab
     416              :       TYPE(gto_basis_set_type), POINTER                  :: fba, fbb
     417              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT)   :: scona_shg, sconb_shg
     418              :       REAL(KIND=dp), INTENT(IN)                          :: omega
     419              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN)      :: Waux_mat
     420              :       REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN)   :: dWaux_mat
     421              :       LOGICAL, INTENT(IN)                                :: calculate_forces
     422              : 
     423              :       INTEGER :: iset, jset, la_max_set, lb_max_set, ndev, nds, nds_max, npgfa_set, &
     424              :                  npgfb_set, nseta, nsetb, nsgfa_set, nsgfb_set, nshella_set, nshellb_set
     425      2528004 :       INTEGER, DIMENSION(:), POINTER                     :: la_max, lb_max, npgfa, npgfb, nsgfa, &
     426      2528004 :                                                             nsgfb, nshella, nshellb
     427      2528004 :       INTEGER, DIMENSION(:, :), POINTER                  :: first_sgfa, first_sgfb, la, lb
     428              :       REAL(KIND=dp)                                      :: dab
     429      2528004 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: set_radius_a, set_radius_b
     430      2528004 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: zeta, zetb
     431              :       REAL(KIND=dp), DIMENSION(:, :, :), ALLOCATABLE     :: swork, swork_cont
     432              : 
     433      2528004 :       NULLIFY (la_max, lb_max, npgfa, npgfb, first_sgfa, first_sgfb, set_radius_a, &
     434      2528004 :                set_radius_b, zeta, zetb)
     435              : 
     436              :       ! basis ikind
     437      2528004 :       first_sgfa => fba%first_sgf
     438      2528004 :       la_max => fba%lmax
     439      2528004 :       la => fba%l
     440      2528004 :       npgfa => fba%npgf
     441      2528004 :       nsgfa => fba%nsgf_set
     442      2528004 :       nseta = fba%nset
     443      2528004 :       set_radius_a => fba%set_radius
     444      2528004 :       zeta => fba%zet
     445      2528004 :       nshella => fba%nshell
     446              :       ! basis jkind
     447      2528004 :       first_sgfb => fbb%first_sgf
     448      2528004 :       lb_max => fbb%lmax
     449      2528004 :       lb => fbb%l
     450      2528004 :       npgfb => fbb%npgf
     451      2528004 :       nsgfb => fbb%nsgf_set
     452      2528004 :       nsetb = fbb%nset
     453      2528004 :       set_radius_b => fbb%set_radius
     454      2528004 :       zetb => fbb%zet
     455      2528004 :       nshellb => fbb%nshell
     456              : 
     457     10112016 :       dab = SQRT(SUM(rab**2))
     458              : 
     459      7385868 :       la_max_set = MAXVAL(la_max)
     460      7386668 :       lb_max_set = MAXVAL(lb_max)
     461              : 
     462              :       ! allocate some work matrices
     463      7385868 :       npgfa_set = MAXVAL(npgfa)
     464      7386668 :       npgfb_set = MAXVAL(npgfb)
     465      7385868 :       nshella_set = MAXVAL(nshella)
     466      7386668 :       nshellb_set = MAXVAL(nshellb)
     467      7385868 :       nsgfa_set = MAXVAL(nsgfa)
     468      7386668 :       nsgfb_set = MAXVAL(nsgfb)
     469      2528004 :       ndev = 0
     470      2528004 :       IF (calculate_forces) ndev = 1
     471      2528004 :       nds_max = la_max_set + lb_max_set + ndev + 1
     472     12640020 :       ALLOCATE (swork(npgfa_set, npgfb_set, nds_max))
     473     12640020 :       ALLOCATE (swork_cont(nds_max, nshella_set, nshellb_set))
     474              : 
     475    158059000 :       vab = 0.0_dp
     476     16986244 :       IF (calculate_forces) dvab = 0.0_dp
     477              : 
     478      7385868 :       DO iset = 1, nseta
     479              : 
     480     27697272 :          DO jset = 1, nsetb
     481              : 
     482     20311404 :             nds = la_max(iset) + lb_max(jset) + ndev + 1
     483    177975240 :             swork(1:npgfa(iset), 1:npgfb(jset), 1:nds) = 0.0_dp
     484              :             CALL s_operator_ab(la_max(iset), npgfa(iset), zeta(:, iset), &
     485              :                                lb_max(jset), npgfb(jset), zetb(:, jset), &
     486     20311404 :                                omega, rab, swork, calculate_forces)
     487              :             CALL contract_sint_ab_chigh(npgfa(iset), nshella(iset), &
     488              :                                         scona_shg(1:npgfa(iset), 1:nshella(iset), iset), &
     489              :                                         npgfb(jset), nshellb(jset), &
     490              :                                         sconb_shg(1:npgfb(jset), 1:nshellb(jset), jset), &
     491              :                                         nds, swork(1:npgfa(iset), 1:npgfb(jset), 1:nds), &
     492     20311404 :                                         swork_cont(1:nds, 1:nshella(iset), 1:nshellb(jset)))
     493              :             CALL construct_int_shg_ab(la(:, iset), first_sgfa(:, iset), nshella(iset), &
     494              :                                       lb(:, jset), first_sgfb(:, jset), nshellb(jset), &
     495     20311404 :                                       swork_cont, Waux_mat, vab)
     496     25169268 :             IF (calculate_forces) THEN
     497              :                !*** -rab, since Eq. in Ref. use Ra-Rb, not Rb-Ra
     498              :                CALL construct_dev_shg_ab(la(:, iset), first_sgfa(:, iset), nshella(iset), &
     499              :                                          lb(:, jset), first_sgfb(:, jset), nshellb(jset), &
     500        96000 :                                          -rab, swork_cont, Waux_mat, dWaux_mat, dvab)
     501              :             END IF
     502              :          END DO
     503              :       END DO
     504              : 
     505      2528004 :       DEALLOCATE (swork, swork_cont)
     506              : 
     507      2528004 :    END SUBROUTINE int_operator_ab_shg_low
     508              : 
     509              : ! **************************************************************************************************
     510              : !> \brief calculate overlap integrals (a,b); requires angular-dependent part as input
     511              : !> \param sab integral (a,b)
     512              : !> \param dsab derivative of sab
     513              : !> \param rab distance vector
     514              : !> \param fba basis at center A
     515              : !> \param fbb basis at center B
     516              : !> \param scona_shg contraction matrix A
     517              : !> \param sconb_shg contraxtion matrix B
     518              : !> \param Waux_mat W matrix that contains the angular-dependent part
     519              : !> \param dWaux_mat derivative of the W matrix
     520              : !> \param calculate_ints ...
     521              : !> \param calculate_forces ...
     522              : !> \param contraction_high ...
     523              : ! **************************************************************************************************
     524        16834 :    SUBROUTINE int_overlap_ab_shg_low(sab, dsab, rab, fba, fbb, scona_shg, sconb_shg, Waux_mat, dWaux_mat, &
     525              :                                      calculate_ints, calculate_forces, contraction_high)
     526              : 
     527              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
     528              :          INTENT(INOUT)                                   :: sab
     529              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
     530              :          INTENT(INOUT)                                   :: dsab
     531              :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rab
     532              :       TYPE(gto_basis_set_type), POINTER                  :: fba, fbb
     533              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN)      :: scona_shg, sconb_shg, Waux_mat
     534              :       REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN)   :: dWaux_mat
     535              :       LOGICAL, INTENT(IN)                                :: calculate_ints, calculate_forces
     536              :       LOGICAL, INTENT(IN), OPTIONAL                      :: contraction_high
     537              : 
     538              :       INTEGER                                            :: iset, jset, la_max_set, lb_max_set, &
     539              :                                                             ndev, nds, nds_max, npgfa_set, &
     540              :                                                             npgfb_set, nseta, nsetb, nshella_set, &
     541              :                                                             nshellb_set
     542        16834 :       INTEGER, DIMENSION(:), POINTER                     :: la_max, lb_max, npgfa, npgfb, nshella, &
     543        16834 :                                                             nshellb
     544        16834 :       INTEGER, DIMENSION(:, :), POINTER                  :: first_sgfa, first_sgfb, la, lb
     545              :       LOGICAL                                            :: my_contraction_high
     546              :       REAL(KIND=dp)                                      :: dab
     547              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: swork, swork_cont
     548        16834 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: set_radius_a, set_radius_b
     549        16834 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: zeta, zetb
     550              : 
     551        16834 :       NULLIFY (la_max, lb_max, npgfa, npgfb, first_sgfa, first_sgfb, set_radius_a, &
     552        16834 :                set_radius_b, zeta, zetb)
     553              : 
     554              :       ! basis ikind
     555        16834 :       first_sgfa => fba%first_sgf
     556        16834 :       la_max => fba%lmax
     557        16834 :       la => fba%l
     558        16834 :       npgfa => fba%npgf
     559        16834 :       nseta = fba%nset
     560        16834 :       set_radius_a => fba%set_radius
     561        16834 :       zeta => fba%zet
     562        16834 :       nshella => fba%nshell
     563              :       ! basis jkind
     564        16834 :       first_sgfb => fbb%first_sgf
     565        16834 :       lb_max => fbb%lmax
     566        16834 :       lb => fbb%l
     567        16834 :       npgfb => fbb%npgf
     568        16834 :       nsetb = fbb%nset
     569        16834 :       set_radius_b => fbb%set_radius
     570        16834 :       zetb => fbb%zet
     571        16834 :       nshellb => fbb%nshell
     572              : 
     573        67336 :       dab = SQRT(SUM(rab**2))
     574              : 
     575        46529 :       la_max_set = MAXVAL(la_max)
     576        46689 :       lb_max_set = MAXVAL(lb_max)
     577              : 
     578              :       ! allocate some work matrices
     579        46529 :       npgfa_set = MAXVAL(npgfa)
     580        46689 :       npgfb_set = MAXVAL(npgfb)
     581        46529 :       nshella_set = MAXVAL(nshella)
     582        46689 :       nshellb_set = MAXVAL(nshellb)
     583        16834 :       ndev = 0
     584        16834 :       IF (calculate_forces) ndev = 1
     585        16834 :       nds_max = la_max_set + lb_max_set + ndev + 1
     586        84170 :       ALLOCATE (swork(npgfa_set, npgfb_set, nds_max))
     587        84170 :       ALLOCATE (swork_cont(nds_max, nshella_set, nshellb_set))
     588              : 
     589     10392820 :       IF (calculate_ints) sab = 0.0_dp
     590     16449973 :       IF (calculate_forces) dsab = 0.0_dp
     591              : 
     592        16834 :       my_contraction_high = .TRUE.
     593        16834 :       IF (PRESENT(contraction_high)) my_contraction_high = contraction_high
     594              : 
     595        46529 :       DO iset = 1, nseta
     596              : 
     597       193520 :          DO jset = 1, nsetb
     598              : 
     599       146991 :             IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
     600              : 
     601       120606 :             nds = la_max(iset) + lb_max(jset) + ndev + 1
     602      4407240 :             swork(1:npgfa(iset), 1:npgfb(jset), 1:nds) = 0.0_dp
     603              :             CALL s_overlap_ab(la_max(iset), npgfa(iset), zeta(:, iset), &
     604              :                               lb_max(jset), npgfb(jset), zetb(:, jset), &
     605       120606 :                               rab, swork, calculate_forces)
     606       120606 :             IF (my_contraction_high) THEN
     607              :                CALL contract_sint_ab_chigh(npgfa(iset), nshella(iset), &
     608              :                                            scona_shg(1:npgfa(iset), 1:nshella(iset), iset), &
     609              :                                            npgfb(jset), nshellb(jset), &
     610              :                                            sconb_shg(1:npgfb(jset), 1:nshellb(jset), jset), &
     611              :                                            nds, swork(1:npgfa(iset), 1:npgfb(jset), 1:nds), &
     612        19819 :                                            swork_cont(1:nds, 1:nshella(iset), 1:nshellb(jset)))
     613              :             ELSE
     614              :                CALL contract_sint_ab_clow(la(:, iset), npgfa(iset), nshella(iset), &
     615              :                                           scona_shg(:, :, iset), &
     616              :                                           lb(:, jset), npgfb(jset), nshellb(jset), &
     617              :                                           sconb_shg(:, :, jset), &
     618       100787 :                                           swork, swork_cont, calculate_forces)
     619              :             END IF
     620       120606 :             IF (calculate_ints) THEN
     621              :                CALL construct_int_shg_ab(la(:, iset), first_sgfa(:, iset), nshella(iset), &
     622              :                                          lb(:, jset), first_sgfb(:, jset), nshellb(jset), &
     623        76964 :                                          swork_cont, Waux_mat, sab)
     624              :             END IF
     625       150301 :             IF (calculate_forces) THEN
     626              :                !*** -rab, since Eq. in Ref. use Ra-Rb, not Rb-Ra
     627              :                CALL construct_dev_shg_ab(la(:, iset), first_sgfa(:, iset), nshella(iset), &
     628              :                                          lb(:, jset), first_sgfb(:, jset), nshellb(jset), &
     629       193768 :                                          -rab, swork_cont, Waux_mat, dWaux_mat, dsab)
     630              :             END IF
     631              :          END DO
     632              :       END DO
     633              : 
     634        16834 :       DEALLOCATE (swork, swork_cont)
     635              : 
     636        16834 :    END SUBROUTINE int_overlap_ab_shg_low
     637              : 
     638              : ! **************************************************************************************************
     639              : !> \brief calculate integrals (a|ra^2m)|b); requires angular-dependent part as input
     640              : !> \param vab integral matrix of spherical contracted Gaussian functions
     641              : !> \param dvab derivative of the integrals
     642              : !> \param rab distance vector between center A and B
     643              : !> \param fba basis at center A
     644              : !> \param fbb basis at center B
     645              : !> \param sconb_shg SHG contraction matrix for B
     646              : !> \param scon_ra2m contraction matrix for A including the combinatorial factors
     647              : !> \param m exponent in (r-Ra)^(2m) operator
     648              : !> \param Waux_mat W matrix that contains the angular-dependent part
     649              : !> \param dWaux_mat derivative of the W matrix
     650              : !> \param calculate_forces ...
     651              : ! **************************************************************************************************
     652           32 :    SUBROUTINE int_ra2m_ab_shg_low(vab, dvab, rab, fba, fbb, sconb_shg, scon_ra2m, m, Waux_mat, dWaux_mat, &
     653              :                                   calculate_forces)
     654              : 
     655              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT)      :: vab
     656              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT)   :: dvab
     657              :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rab
     658              :       TYPE(gto_basis_set_type), POINTER                  :: fba, fbb
     659              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN)      :: sconb_shg
     660              :       REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN)   :: scon_ra2m
     661              :       INTEGER, INTENT(IN)                                :: m
     662              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN)      :: Waux_mat
     663              :       REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN)   :: dWaux_mat
     664              :       LOGICAL, INTENT(IN)                                :: calculate_forces
     665              : 
     666              :       INTEGER                                            :: iset, jset, la_max_set, lb_max_set, &
     667              :                                                             ndev, nds, nds_max, npgfa_set, &
     668              :                                                             npgfb_set, nseta, nsetb, nshella_set, &
     669              :                                                             nshellb_set
     670           32 :       INTEGER, DIMENSION(:), POINTER                     :: la_max, lb_max, npgfa, npgfb, nsgfa, &
     671           32 :                                                             nsgfb, nshella, nshellb
     672           32 :       INTEGER, DIMENSION(:, :), POINTER                  :: first_sgfa, first_sgfb, la, lb
     673              :       REAL(KIND=dp)                                      :: dab
     674              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: swork_cont
     675              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :)  :: swork
     676           32 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: zeta, zetb
     677              : 
     678           32 :       NULLIFY (la_max, lb_max, npgfa, npgfb, first_sgfa, first_sgfb, zeta, zetb)
     679              : 
     680              :       ! basis ikind
     681           32 :       first_sgfa => fba%first_sgf
     682           32 :       la_max => fba%lmax
     683           32 :       la => fba%l
     684           32 :       npgfa => fba%npgf
     685           32 :       nsgfa => fba%nsgf_set
     686           32 :       nseta = fba%nset
     687           32 :       zeta => fba%zet
     688           32 :       nshella => fba%nshell
     689              :       ! basis jkind
     690           32 :       first_sgfb => fbb%first_sgf
     691           32 :       lb_max => fbb%lmax
     692           32 :       lb => fbb%l
     693           32 :       npgfb => fbb%npgf
     694           32 :       nsgfb => fbb%nsgf_set
     695           32 :       nsetb = fbb%nset
     696           32 :       zetb => fbb%zet
     697           32 :       nshellb => fbb%nshell
     698              : 
     699          128 :       dab = SQRT(SUM(rab**2))
     700              : 
     701          352 :       la_max_set = MAXVAL(la_max)
     702          512 :       lb_max_set = MAXVAL(lb_max)
     703              : 
     704              :       ! allocate some work matrices
     705          352 :       npgfa_set = MAXVAL(npgfa)
     706          512 :       npgfb_set = MAXVAL(npgfb)
     707          352 :       nshella_set = MAXVAL(nshella)
     708          512 :       nshellb_set = MAXVAL(nshellb)
     709           32 :       ndev = 0
     710           32 :       IF (calculate_forces) ndev = 1
     711           32 :       nds_max = la_max_set + lb_max_set + ndev + 1
     712          192 :       ALLOCATE (swork(npgfa_set, npgfb_set, 1:m + 1, nds_max))
     713          160 :       ALLOCATE (swork_cont(nds_max, nshella_set, nshellb_set))
     714              : 
     715       963872 :       vab = 0.0_dp
     716      2891680 :       IF (calculate_forces) dvab = 0.0_dp
     717              : 
     718          352 :       DO iset = 1, nseta
     719              : 
     720         5152 :          DO jset = 1, nsetb
     721              : 
     722         4800 :             nds = la_max(iset) + lb_max(jset) + ndev + 1
     723       447840 :             swork(1:npgfa(iset), 1:npgfb(jset), 1:m + 1, 1:nds) = 0.0_dp
     724              :             CALL s_ra2m_ab(la_max(iset), npgfa(iset), zeta(:, iset), &
     725              :                            lb_max(jset), npgfb(jset), zetb(:, jset), &
     726         4800 :                            m, rab, swork, calculate_forces)
     727              :             CALL contract_s_ra2m_ab(npgfa(iset), nshella(iset), &
     728              :                                     scon_ra2m(1:npgfa(iset), 1:m + 1, 1:nshella(iset), iset), &
     729              :                                     npgfb(jset), nshellb(jset), &
     730              :                                     sconb_shg(1:npgfb(jset), 1:nshellb(jset), jset), &
     731              :                                     swork(1:npgfa(iset), 1:npgfb(jset), 1:m + 1, 1:nds), &
     732              :                                     swork_cont(1:nds, 1:nshella(iset), 1:nshellb(jset)), &
     733         4800 :                                     m, nds)
     734              :             CALL construct_int_shg_ab(la(:, iset), first_sgfa(:, iset), nshella(iset), &
     735              :                                       lb(:, jset), first_sgfb(:, jset), nshellb(jset), &
     736         4800 :                                       swork_cont, Waux_mat, vab)
     737         5120 :             IF (calculate_forces) THEN
     738              :                !*** -rab, since Eq. in Ref. use Ra-Rb, not Rb-Ra
     739              :                CALL construct_dev_shg_ab(la(:, iset), first_sgfa(:, iset), nshella(iset), &
     740              :                                          lb(:, jset), first_sgfb(:, jset), nshellb(jset), &
     741        19200 :                                          -rab, swork_cont, Waux_mat, dWaux_mat, dvab)
     742              :             END IF
     743              :          END DO
     744              :       END DO
     745              : 
     746           32 :       DEALLOCATE (swork, swork_cont)
     747              : 
     748           32 :    END SUBROUTINE int_ra2m_ab_shg_low
     749              : ! **************************************************************************************************
     750              : !> \brief calculate integrals (a,b,fb); requires angular-dependent part as input
     751              : !> \param abbint integral (a,b,fb)
     752              : !> \param dabbint derivative of abbint
     753              : !> \param rab distance vector between A and B
     754              : !> \param oba orbital basis at center A
     755              : !> \param obb orbital basis at center B
     756              : !> \param fbb auxiliary basis set at center B
     757              : !> \param scon_oba contraction matrix for orb bas on A
     758              : !> \param sconb_mix mixed contraction matrix orb + ri basis on B
     759              : !> \param obb_index orbital basis index for sconb_mix
     760              : !> \param fbb_index ri basis index for sconb_mix
     761              : !> \param cg_coeff Clebsch-Gordon coefficients
     762              : !> \param cg_none0_list list of none-zero Clebsch-Gordon coefficients
     763              : !> \param ncg_none0 number of non-zero Clebsch-Gordon coefficients
     764              : !> \param Waux_mat W matrix that contains the angular-dependent part
     765              : !> \param dWaux_mat derivative of the W matrix
     766              : !> \param calculate_ints ...
     767              : !> \param calculate_forces ...
     768              : ! **************************************************************************************************
     769          744 :    SUBROUTINE int_overlap_abb_shg_low(abbint, dabbint, rab, oba, obb, fbb, scon_oba, sconb_mix, &
     770          744 :                                       obb_index, fbb_index, cg_coeff, cg_none0_list, &
     771          744 :                                       ncg_none0, Waux_mat, dWaux_mat, calculate_ints, calculate_forces)
     772              : 
     773              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
     774              :          INTENT(INOUT)                                   :: abbint
     775              :       REAL(KIND=dp), ALLOCATABLE, &
     776              :          DIMENSION(:, :, :, :), INTENT(INOUT)            :: dabbint
     777              :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rab
     778              :       TYPE(gto_basis_set_type), POINTER                  :: oba, obb, fbb
     779              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN)      :: scon_oba
     780              :       REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN)   :: sconb_mix
     781              :       INTEGER, DIMENSION(:, :, :), INTENT(IN)            :: obb_index, fbb_index
     782              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN)      :: cg_coeff
     783              :       INTEGER, DIMENSION(:, :, :), INTENT(IN)            :: cg_none0_list
     784              :       INTEGER, DIMENSION(:, :), INTENT(IN)               :: ncg_none0
     785              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN)      :: Waux_mat
     786              :       REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN)   :: dWaux_mat
     787              :       LOGICAL, INTENT(IN)                                :: calculate_ints, calculate_forces
     788              : 
     789              :       INTEGER :: iset, jset, kset, la_max_set, lb_max_set, lbb_max, lbb_max_set, lcb_max_set, na, &
     790              :          nb, ncb, ndev, nds, nds_max, nl, nl_set, npgfa_set, npgfb_set, npgfcb_set, nseta, nsetb, &
     791              :          nsetcb, nshella_set, nshellb_set, nshellcb_set, sgfa, sgfb, sgfcb
     792          744 :       INTEGER, DIMENSION(:), POINTER                     :: la_max, lb_max, lcb_max, npgfa, npgfb, &
     793          744 :                                                             npgfcb, nsgfa, nsgfb, nsgfcb, nshella, &
     794          744 :                                                             nshellb, nshellcb
     795          744 :       INTEGER, DIMENSION(:, :), POINTER                  :: first_sgfa, first_sgfb, first_sgfcb, la, &
     796          744 :                                                             lb, lcb
     797              :       REAL(KIND=dp)                                      :: dab, rab2
     798              :       REAL(KIND=dp), ALLOCATABLE, &
     799              :          DIMENSION(:, :, :, :, :)                        :: swork, swork_cont
     800          744 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: set_radius_a, set_radius_b, set_radius_cb
     801          744 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: zeta, zetb, zetcb
     802              : 
     803          744 :       NULLIFY (la_max, lb_max, lcb_max, npgfa, npgfb, npgfcb)
     804          744 :       NULLIFY (first_sgfa, first_sgfb, first_sgfcb, set_radius_a, set_radius_b, &
     805          744 :                set_radius_cb, zeta, zetb, zetcb)
     806              : 
     807              :       ! basis ikind
     808          744 :       first_sgfa => oba%first_sgf
     809          744 :       la_max => oba%lmax
     810          744 :       la => oba%l
     811          744 :       nsgfa => oba%nsgf_set
     812          744 :       npgfa => oba%npgf
     813          744 :       nshella => oba%nshell
     814          744 :       nseta = oba%nset
     815          744 :       set_radius_a => oba%set_radius
     816          744 :       zeta => oba%zet
     817              :       ! basis jkind
     818          744 :       first_sgfb => obb%first_sgf
     819          744 :       lb_max => obb%lmax
     820          744 :       lb => obb%l
     821          744 :       nsgfb => obb%nsgf_set
     822          744 :       npgfb => obb%npgf
     823          744 :       nshellb => obb%nshell
     824          744 :       nsetb = obb%nset
     825          744 :       set_radius_b => obb%set_radius
     826          744 :       zetb => obb%zet
     827              : 
     828              :       ! basis RI on B
     829          744 :       first_sgfcb => fbb%first_sgf
     830          744 :       lcb_max => fbb%lmax
     831          744 :       lcb => fbb%l
     832          744 :       nsgfcb => fbb%nsgf_set
     833          744 :       npgfcb => fbb%npgf
     834          744 :       nshellcb => fbb%nshell
     835          744 :       nsetcb = fbb%nset
     836          744 :       set_radius_cb => fbb%set_radius
     837          744 :       zetcb => fbb%zet
     838              : 
     839         2976 :       dab = SQRT(SUM(rab**2))
     840          744 :       rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
     841              : 
     842         1716 :       la_max_set = MAXVAL(la_max)
     843         1716 :       lb_max_set = MAXVAL(lb_max)
     844         7242 :       lcb_max_set = MAXVAL(lcb_max)
     845         1716 :       npgfa_set = MAXVAL(npgfa)
     846         1716 :       npgfb_set = MAXVAL(npgfb)
     847         7242 :       npgfcb_set = MAXVAL(npgfcb)
     848         1716 :       nshella_set = MAXVAL(nshella)
     849         1716 :       nshellb_set = MAXVAL(nshellb)
     850         7242 :       nshellcb_set = MAXVAL(nshellcb)
     851              :       !*** for forces: derivative+1 in auxiliary vector required
     852          744 :       ndev = 0
     853          744 :       IF (calculate_forces) ndev = 1
     854              : 
     855          744 :       lbb_max_set = lb_max_set + lcb_max_set
     856              : 
     857              :       ! allocate some work storage....
     858          744 :       nds_max = la_max_set + lbb_max_set + ndev + 1
     859          744 :       nl_set = INT((lbb_max_set)/2)
     860         5208 :       ALLOCATE (swork(npgfa_set, npgfb_set, npgfcb_set, nl_set + 1, nds_max))
     861         5208 :       ALLOCATE (swork_cont(nds_max, 0:nl_set, nshella_set, nshellb_set, nshellcb_set))
     862              : 
     863         1716 :       DO iset = 1, nseta
     864              : 
     865         3144 :          DO jset = 1, nsetb
     866              : 
     867         1428 :             IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
     868              : 
     869         8664 :             DO kset = 1, nsetcb
     870              : 
     871         6816 :                IF (set_radius_a(iset) + set_radius_cb(kset) < dab) CYCLE
     872              : 
     873         4004 :                lbb_max = lb_max(jset) + lcb_max(kset)
     874         4004 :                nds = la_max(iset) + lbb_max + ndev + 1
     875         4004 :                nl = INT((lbb_max)/2) + 1
     876      4749072 :                swork(1:npgfa(iset), 1:npgfb(jset), 1:npgfcb(kset), 1:nl, 1:nds) = 0.0_dp
     877              :                CALL s_overlap_abb(la_max(iset), npgfa(iset), zeta(:, iset), &
     878              :                                   lb_max(jset), npgfb(jset), zetb(:, jset), &
     879              :                                   lcb_max(kset), npgfcb(kset), zetcb(:, kset), &
     880         4004 :                                   rab, swork, calculate_forces)
     881              : 
     882              :                CALL contract_s_overlap_abb(la(:, iset), npgfa(iset), nshella(iset), &
     883              :                                            scon_oba(1:npgfa(iset), 1:nshella(iset), iset), &
     884              :                                            lb(:, jset), npgfb(jset), nshellb(jset), &
     885              :                                            lcb(:, kset), npgfcb(kset), nshellcb(kset), &
     886              :                                            obb_index(:, :, jset), fbb_index(:, :, kset), sconb_mix, nl, nds, &
     887              :                                            swork(1:npgfa(iset), 1:npgfb(jset), 1:npgfcb(kset), 1:nl, 1:nds), &
     888         4004 :                                            swork_cont, calculate_forces)
     889              : 
     890         4004 :                IF (calculate_ints) THEN
     891              :                   CALL construct_overlap_shg_abb(la(:, iset), first_sgfa(:, iset), nshella(iset), &
     892              :                                                  lb(:, jset), first_sgfb(:, jset), nshellb(jset), &
     893              :                                                  lcb(:, kset), first_sgfcb(:, kset), nshellcb(kset), &
     894              :                                                  cg_coeff, cg_none0_list, &
     895         3277 :                                                  ncg_none0, swork_cont, Waux_mat, abbint)
     896              :                END IF
     897         4004 :                IF (calculate_forces) THEN
     898              :                   !*** -rab, since Eq. in Ref. use Ra-Rb, not Rb-Ra
     899              :                   CALL dev_overlap_shg_abb(la(:, iset), first_sgfa(:, iset), nshella(iset), &
     900              :                                            lb(:, jset), first_sgfb(:, jset), nshellb(jset), &
     901              :                                            lcb(:, kset), first_sgfcb(:, kset), nshellcb(kset), &
     902              :                                            cg_coeff, cg_none0_list, ncg_none0, -rab, swork_cont, &
     903         3868 :                                            Waux_mat, dWaux_mat, dabbint)
     904              :                END IF
     905              :                ! max value of integrals in this set triple
     906         4004 :                sgfa = first_sgfa(1, iset)
     907         4004 :                na = sgfa + nsgfa(iset) - 1
     908         4004 :                sgfb = first_sgfb(1, jset)
     909         4004 :                nb = sgfb + nsgfb(jset) - 1
     910         4004 :                sgfcb = first_sgfcb(1, kset)
     911         8244 :                ncb = sgfcb + nsgfcb(kset) - 1
     912              :             END DO
     913              :          END DO
     914              :       END DO
     915              : 
     916          744 :       DEALLOCATE (swork_cont)
     917          744 :       DEALLOCATE (swork)
     918              : 
     919          744 :    END SUBROUTINE int_overlap_abb_shg_low
     920              : ! **************************************************************************************************
     921              : !> \brief obtain integrals (a,b,fb) by symmetry relations from (a,b,fa) if basis sets at a and
     922              : !>        b are of the same kind, i.e. a and b are same atom type
     923              : !> \param abbint integral (a,b,fb)
     924              : !> \param dabbint derivative of abbint
     925              : !> \param abaint integral (a,b,fa)
     926              : !> \param dabdaint derivative of abaint
     927              : !> \param rab distance vector between A and B
     928              : !> \param oba orbital basis at center A
     929              : !> \param fba auxiliary basis set at center A
     930              : !> \param calculate_ints ...
     931              : !> \param calculate_forces ...
     932              : ! **************************************************************************************************
     933        14005 :    SUBROUTINE get_abb_same_kind(abbint, dabbint, abaint, dabdaint, rab, oba, fba, &
     934              :                                 calculate_ints, calculate_forces)
     935              : 
     936              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
     937              :          INTENT(INOUT)                                   :: abbint
     938              :       REAL(KIND=dp), ALLOCATABLE, &
     939              :          DIMENSION(:, :, :, :), INTENT(INOUT)            :: dabbint
     940              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
     941              :          INTENT(INOUT)                                   :: abaint
     942              :       REAL(KIND=dp), ALLOCATABLE, &
     943              :          DIMENSION(:, :, :, :), INTENT(INOUT)            :: dabdaint
     944              :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rab
     945              :       TYPE(gto_basis_set_type), POINTER                  :: oba, fba
     946              :       LOGICAL, INTENT(IN)                                :: calculate_ints, calculate_forces
     947              : 
     948              :       INTEGER :: i, iend, iset, isgfa, ishella, istart, jend, jset, jsgfa, jshella, jstart, kend, &
     949              :          kset, ksgfa, kshella, kstart, lai, laj, lak, nseta, nsetca, nsgfa, nsgfca, sgfa_end, &
     950              :          sgfa_start
     951        14005 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: lai_set, laj_set, lak_set
     952        14005 :       INTEGER, DIMENSION(:), POINTER                     :: nsgfa_set, nsgfca_set, nshella, nshellca
     953        14005 :       INTEGER, DIMENSION(:, :), POINTER                  :: first_sgfa, first_sgfca, la, lca
     954              :       REAL(KIND=dp)                                      :: dab
     955        14005 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: set_radius_a, set_radius_ca
     956              : 
     957        14005 :       NULLIFY (nshellca, first_sgfa, first_sgfca, lca, set_radius_a, &
     958        14005 :                set_radius_ca)
     959              : 
     960              :       ! basis ikind
     961        14005 :       first_sgfa => oba%first_sgf
     962        14005 :       set_radius_a => oba%set_radius
     963        14005 :       nseta = oba%nset
     964        14005 :       nsgfa = oba%nsgf
     965        14005 :       nsgfa_set => oba%nsgf_set
     966        14005 :       nshella => oba%nshell
     967        14005 :       la => oba%l
     968              : 
     969              :       ! basis RI
     970        14005 :       first_sgfca => fba%first_sgf
     971        14005 :       set_radius_ca => fba%set_radius
     972        14005 :       nsetca = fba%nset
     973        14005 :       nshellca => fba%nshell
     974        14005 :       nsgfca = fba%nsgf
     975        14005 :       nsgfca_set => fba%nsgf_set
     976        14005 :       lca => fba%l
     977              : 
     978        42015 :       ALLOCATE (lai_set(nsgfa))
     979        28010 :       ALLOCATE (laj_set(nsgfa))
     980        42015 :       ALLOCATE (lak_set(nsgfca))
     981              : 
     982        56020 :       dab = SQRT(SUM(rab**2))
     983        28526 :       DO iset = 1, nseta
     984              : 
     985        82460 :          DO ishella = 1, nshella(iset)
     986        67939 :             sgfa_start = first_sgfa(ishella, iset)
     987        67939 :             sgfa_end = sgfa_start + 2*la(ishella, iset)
     988       259275 :             lai_set(sgfa_start:sgfa_end) = la(ishella, iset)
     989              :          END DO
     990        14521 :          istart = first_sgfa(1, iset)
     991        14521 :          iend = istart + nsgfa_set(iset) - 1
     992              : 
     993        44079 :          DO jset = 1, nseta
     994              : 
     995        15553 :             IF (set_radius_a(iset) + set_radius_a(jset) < dab) CYCLE
     996        81336 :             DO jshella = 1, nshella(jset)
     997        67177 :                sgfa_start = first_sgfa(jshella, jset)
     998        67177 :                sgfa_end = sgfa_start + 2*la(jshella, jset)
     999       252725 :                laj_set(sgfa_start:sgfa_end) = la(jshella, jset)
    1000              :             END DO
    1001        14159 :             jstart = first_sgfa(1, jset)
    1002        14159 :             jend = jstart + nsgfa_set(jset) - 1
    1003              : 
    1004       152516 :             DO kset = 1, nsetca
    1005              : 
    1006       123836 :                IF (set_radius_a(iset) + set_radius_ca(kset) < dab) CYCLE
    1007              : 
    1008       210208 :                DO kshella = 1, nshellca(kset)
    1009       150740 :                   sgfa_start = first_sgfca(kshella, kset)
    1010       150740 :                   sgfa_end = sgfa_start + 2*lca(kshella, kset)
    1011       732074 :                   lak_set(sgfa_start:sgfa_end) = lca(kshella, kset)
    1012              :                END DO
    1013        59468 :                kstart = first_sgfca(1, kset)
    1014        59468 :                kend = kstart + nsgfca_set(kset) - 1
    1015       596887 :                DO ksgfa = kstart, kend
    1016       521866 :                   lak = lak_set(ksgfa)
    1017      7050400 :                   DO jsgfa = jstart, jend
    1018      6404698 :                      laj = laj_set(jsgfa)
    1019     88237182 :                      DO isgfa = istart, iend
    1020     81310618 :                         lai = lai_set(isgfa)
    1021     87715316 :                         IF (MODULO((lai + laj + lak), 2) /= 0) THEN
    1022     40645729 :                            IF (calculate_ints) THEN
    1023              :                               abbint(isgfa, jsgfa, ksgfa) = &
    1024     22047455 :                                  -abaint(jsgfa, isgfa, ksgfa)
    1025              :                            END IF
    1026     40645729 :                            IF (calculate_forces) THEN
    1027     74393096 :                               DO i = 1, 3
    1028              :                                  dabbint(isgfa, jsgfa, ksgfa, i) = &
    1029     74393096 :                                     -dabdaint(jsgfa, isgfa, ksgfa, i)
    1030              :                               END DO
    1031              :                            END IF
    1032              :                         ELSE
    1033     40664889 :                            IF (calculate_ints) THEN
    1034              :                               abbint(isgfa, jsgfa, ksgfa) = &
    1035     22054833 :                                  abaint(jsgfa, isgfa, ksgfa)
    1036              :                            END IF
    1037     40664889 :                            IF (calculate_forces) THEN
    1038     74440224 :                               DO i = 1, 3
    1039              :                                  dabbint(isgfa, jsgfa, ksgfa, i) = &
    1040     74440224 :                                     dabdaint(jsgfa, isgfa, ksgfa, i)
    1041              :                               END DO
    1042              :                            END IF
    1043              :                         END IF
    1044              :                      END DO
    1045              :                   END DO
    1046              :                END DO
    1047              :             END DO
    1048              :          END DO
    1049              :       END DO
    1050              : 
    1051        14005 :       DEALLOCATE (lai_set, laj_set, lak_set)
    1052              : 
    1053        14005 :    END SUBROUTINE get_abb_same_kind
    1054              : 
    1055              : ! **************************************************************************************************
    1056              : !> \brief calculate integrals (a,b,fa);  requires angular-dependent part as input
    1057              : !> \param abaint integral (a,b,fa)
    1058              : !> \param dabdaint ...
    1059              : !> \param rab distance vector between A and B
    1060              : !> \param oba orbital basis at center A
    1061              : !> \param obb orbital basis at center B
    1062              : !> \param fba auxiliary basis set at center A
    1063              : !> \param scon_obb contraction matrix for orb bas on B
    1064              : !> \param scona_mix mixed contraction matrix orb + ri basis on A
    1065              : !> \param oba_index orbital basis index for scona_mix
    1066              : !> \param fba_index ri basis index for scona_mix
    1067              : !> \param cg_coeff Clebsch-Gordon coefficients
    1068              : !> \param cg_none0_list list of none-zero Clebsch-Gordon coefficients
    1069              : !> \param ncg_none0 number of non-zero Clebsch-Gordon coefficients
    1070              : !> \param Waux_mat W matrix that contains the angular-dependent part
    1071              : !> \param dWaux_mat derivative of the W matrix
    1072              : !> \param calculate_ints ...
    1073              : !> \param calculate_forces ...
    1074              : ! **************************************************************************************************
    1075        14877 :    SUBROUTINE int_overlap_aba_shg_low(abaint, dabdaint, rab, oba, obb, fba, scon_obb, scona_mix, &
    1076        14877 :                                       oba_index, fba_index, cg_coeff, cg_none0_list, &
    1077        14877 :                                       ncg_none0, Waux_mat, dWaux_mat, calculate_ints, calculate_forces)
    1078              : 
    1079              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
    1080              :          INTENT(INOUT)                                   :: abaint
    1081              :       REAL(KIND=dp), ALLOCATABLE, &
    1082              :          DIMENSION(:, :, :, :), INTENT(INOUT)            :: dabdaint
    1083              :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rab
    1084              :       TYPE(gto_basis_set_type), POINTER                  :: oba, obb, fba
    1085              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN)      :: scon_obb
    1086              :       REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN)   :: scona_mix
    1087              :       INTEGER, DIMENSION(:, :, :), INTENT(IN)            :: oba_index, fba_index
    1088              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN)      :: cg_coeff
    1089              :       INTEGER, DIMENSION(:, :, :), INTENT(IN)            :: cg_none0_list
    1090              :       INTEGER, DIMENSION(:, :), INTENT(IN)               :: ncg_none0
    1091              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN)      :: Waux_mat
    1092              :       REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN)   :: dWaux_mat
    1093              :       LOGICAL, INTENT(IN)                                :: calculate_ints, calculate_forces
    1094              : 
    1095              :       INTEGER :: iset, jset, kset, la_max_set, laa_max, laa_max_set, lb_max_set, lca_max_set, na, &
    1096              :          nb, nca, ndev, nds, nds_max, nl, nl_set, npgfa_set, npgfb_set, npgfca_set, nseta, nsetb, &
    1097              :          nsetca, nshella_set, nshellb_set, nshellca_set, sgfa, sgfb, sgfca
    1098        14877 :       INTEGER, DIMENSION(:), POINTER                     :: la_max, lb_max, lca_max, npgfa, npgfb, &
    1099        14877 :                                                             npgfca, nsgfa, nsgfb, nsgfca, nshella, &
    1100        14877 :                                                             nshellb, nshellca
    1101        14877 :       INTEGER, DIMENSION(:, :), POINTER                  :: first_sgfa, first_sgfb, first_sgfca, la, &
    1102        14877 :                                                             lb, lca
    1103              :       REAL(KIND=dp)                                      :: dab, rab2
    1104              :       REAL(KIND=dp), ALLOCATABLE, &
    1105              :          DIMENSION(:, :, :, :, :)                        :: swork, swork_cont
    1106        14877 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: set_radius_a, set_radius_b, set_radius_ca
    1107        14877 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: zeta, zetb, zetca
    1108              : 
    1109        14877 :       NULLIFY (la_max, lb_max, lca_max, npgfa, npgfb, npgfca)
    1110        14877 :       NULLIFY (first_sgfa, first_sgfb, first_sgfca, set_radius_a, set_radius_b, &
    1111        14877 :                set_radius_ca, zeta, zetb, zetca)
    1112              : 
    1113              :       ! basis ikind
    1114        14877 :       first_sgfa => oba%first_sgf
    1115        14877 :       la_max => oba%lmax
    1116        14877 :       la => oba%l
    1117        14877 :       nsgfa => oba%nsgf_set
    1118        14877 :       npgfa => oba%npgf
    1119        14877 :       nshella => oba%nshell
    1120        14877 :       nseta = oba%nset
    1121        14877 :       set_radius_a => oba%set_radius
    1122        14877 :       zeta => oba%zet
    1123              :       ! basis jkind
    1124        14877 :       first_sgfb => obb%first_sgf
    1125        14877 :       lb_max => obb%lmax
    1126        14877 :       lb => obb%l
    1127        14877 :       nsgfb => obb%nsgf_set
    1128        14877 :       npgfb => obb%npgf
    1129        14877 :       nshellb => obb%nshell
    1130        14877 :       nsetb = obb%nset
    1131        14877 :       set_radius_b => obb%set_radius
    1132        14877 :       zetb => obb%zet
    1133              : 
    1134              :       ! basis RI A
    1135        14877 :       first_sgfca => fba%first_sgf
    1136        14877 :       lca_max => fba%lmax
    1137        14877 :       lca => fba%l
    1138        14877 :       nsgfca => fba%nsgf_set
    1139        14877 :       npgfca => fba%npgf
    1140        14877 :       nshellca => fba%nshell
    1141        14877 :       nsetca = fba%nset
    1142        14877 :       set_radius_ca => fba%set_radius
    1143        14877 :       zetca => fba%zet
    1144              : 
    1145        59508 :       dab = SQRT(SUM(rab**2))
    1146        14877 :       rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
    1147              : 
    1148        30542 :       la_max_set = MAXVAL(la_max)
    1149        30542 :       lb_max_set = MAXVAL(lb_max)
    1150       146243 :       lca_max_set = MAXVAL(lca_max)
    1151        30542 :       npgfa_set = MAXVAL(npgfa)
    1152        30542 :       npgfb_set = MAXVAL(npgfb)
    1153       146243 :       npgfca_set = MAXVAL(npgfca)
    1154        30542 :       nshella_set = MAXVAL(nshella)
    1155        30542 :       nshellb_set = MAXVAL(nshellb)
    1156       146243 :       nshellca_set = MAXVAL(nshellca)
    1157              :       !*** for forces: derivative+1 in auxiliary vector required
    1158        14877 :       ndev = 0
    1159        14877 :       IF (calculate_forces) ndev = 1
    1160              : 
    1161        14877 :       laa_max_set = la_max_set + lca_max_set
    1162              : 
    1163              :       ! allocate some work storage....
    1164        14877 :       nds_max = laa_max_set + lb_max_set + ndev + 1
    1165        14877 :       nl_set = INT((laa_max_set)/2)
    1166       104139 :       ALLOCATE (swork(npgfb_set, npgfa_set, npgfca_set, nl_set + 1, nds_max))
    1167       104139 :       ALLOCATE (swork_cont(nds_max, 0:nl_set, nshella_set, nshellb_set, nshellca_set))
    1168              : 
    1169        30542 :       DO iset = 1, nseta
    1170              : 
    1171        47879 :          DO jset = 1, nsetb
    1172              : 
    1173        17337 :             IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
    1174              : 
    1175       165424 :             DO kset = 1, nsetca
    1176              : 
    1177       134368 :                IF (set_radius_b(jset) + set_radius_ca(kset) < dab) CYCLE
    1178              : 
    1179              :                !*** calculate s_baa here
    1180        67188 :                laa_max = la_max(iset) + lca_max(kset)
    1181        67188 :                nds = laa_max + lb_max(jset) + ndev + 1
    1182        67188 :                nl = INT(laa_max/2) + 1
    1183     55180887 :                swork(1:npgfb(jset), 1:npgfa(iset), 1:npgfca(kset), 1:nl, 1:nds) = 0.0_dp
    1184              :                CALL s_overlap_abb(lb_max(jset), npgfb(jset), zetb(:, jset), &
    1185              :                                   la_max(iset), npgfa(iset), zeta(:, iset), &
    1186              :                                   lca_max(kset), npgfca(kset), zetca(:, kset), &
    1187        67188 :                                   rab, swork, calculate_forces)
    1188              : 
    1189              :                CALL contract_s_overlap_aba(la(:, iset), npgfa(iset), nshella(iset), &
    1190              :                                            lb(:, jset), npgfb(jset), nshellb(jset), &
    1191              :                                            scon_obb(1:npgfb(jset), 1:nshellb(jset), jset), &
    1192              :                                            lca(:, kset), npgfca(kset), nshellca(kset), &
    1193              :                                            oba_index(:, :, iset), fba_index(:, :, kset), scona_mix, nl, nds, &
    1194              :                                            swork(1:npgfb(jset), 1:npgfa(iset), 1:npgfca(kset), 1:nl, 1:nds), &
    1195        67188 :                                            swork_cont, calculate_forces)
    1196        67188 :                IF (calculate_ints) THEN
    1197              :                   CALL construct_overlap_shg_aba(la(:, iset), first_sgfa(:, iset), nshella(iset), &
    1198              :                                                  lb(:, jset), first_sgfb(:, jset), nshellb(jset), &
    1199              :                                                  lca(:, kset), first_sgfca(:, kset), nshellca(kset), &
    1200              :                                                  cg_coeff, cg_none0_list, ncg_none0, &
    1201        40040 :                                                  swork_cont, Waux_mat, abaint)
    1202              :                END IF
    1203        67188 :                IF (calculate_forces) THEN
    1204              :                   !*** -rab, since Eq. in Ref. use Ra-Rb, not Rb-Ra
    1205              :                   CALL dev_overlap_shg_aba(la(:, iset), first_sgfa(:, iset), nshella(iset), &
    1206              :                                            lb(:, jset), first_sgfb(:, jset), nshellb(jset), &
    1207              :                                            lca(:, kset), first_sgfca(:, kset), nshellca(kset), &
    1208              :                                            cg_coeff, cg_none0_list, ncg_none0, &
    1209       109552 :                                            -rab, swork_cont, Waux_mat, dWaux_mat, dabdaint)
    1210              :                END IF
    1211              :                ! max value of integrals in this set triple
    1212        67188 :                sgfa = first_sgfa(1, iset)
    1213        67188 :                na = sgfa + nsgfa(iset) - 1
    1214        67188 :                sgfb = first_sgfb(1, jset)
    1215        67188 :                nb = sgfb + nsgfb(jset) - 1
    1216        67188 :                sgfca = first_sgfca(1, kset)
    1217       151705 :                nca = sgfca + nsgfca(kset) - 1
    1218              : 
    1219              :             END DO
    1220              :          END DO
    1221              :       END DO
    1222              : 
    1223        14877 :       DEALLOCATE (swork_cont)
    1224        14877 :       DEALLOCATE (swork)
    1225              : 
    1226        14877 :    END SUBROUTINE int_overlap_aba_shg_low
    1227              : 
    1228              : ! **************************************************************************************************
    1229              : !> \brief precalculates the angular part of the SHG integrals for the matrices
    1230              : !>        (fa,fb), (a,b), (a,b,fa) and (b,fb,a); the same Waux_mat can then be used for all
    1231              : !>        for integrals; specific for LRIGPW
    1232              : !> \param oba orbital basis on a
    1233              : !> \param obb orbital basis on b
    1234              : !> \param fba aux basis on a
    1235              : !> \param fbb aux basis on b
    1236              : !> \param rab distance vector between a and b
    1237              : !> \param Waux_mat W matrix that contains the angular-dependent part
    1238              : !> \param dWaux_mat derivative of the W matrix
    1239              : !> \param calculate_forces ...
    1240              : ! **************************************************************************************************
    1241        14733 :    SUBROUTINE lri_precalc_angular_shg_part(oba, obb, fba, fbb, rab, Waux_mat, dWaux_mat, calculate_forces)
    1242              : 
    1243              :       TYPE(gto_basis_set_type), POINTER                  :: oba, obb, fba, fbb
    1244              :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rab
    1245              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
    1246              :          INTENT(OUT)                                     :: Waux_mat
    1247              :       REAL(KIND=dp), ALLOCATABLE, &
    1248              :          DIMENSION(:, :, :, :), INTENT(OUT)              :: dWaux_mat
    1249              :       LOGICAL, INTENT(IN)                                :: calculate_forces
    1250              : 
    1251              :       INTEGER                                            :: i, isize, j, k, la_max, laa_max, lb_max, &
    1252              :                                                             lbb_max, lca_max, lcb_max, li_max, &
    1253              :                                                             lj_max, lmax, mdim(3), size_int(4, 2), &
    1254              :                                                             temp
    1255        14733 :       INTEGER, DIMENSION(:), POINTER                     :: li_max_all
    1256              :       REAL(KIND=dp)                                      :: rab2
    1257        14733 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: Rc, Rs
    1258              : 
    1259        14733 :       rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
    1260              : 
    1261              :       !*** 1 Waux_mat of size (li_max,lj_max) for elements
    1262              :       !                    i        j
    1263              :       !    [aab]    --> (laa_max, lb_max)
    1264              :       !    [bba]    --> (lbb_max, la_max) --> use for [abb]
    1265              :       !    [ab] ri  --> (lca_max, lcb_max)
    1266              :       !    [ab] orb --> (la_max , lb_max)
    1267              : 
    1268        30210 :       la_max = MAXVAL(oba%lmax)
    1269        30210 :       lb_max = MAXVAL(obb%lmax)
    1270       144483 :       lca_max = MAXVAL(fba%lmax)
    1271       144483 :       lcb_max = MAXVAL(fbb%lmax)
    1272              : 
    1273        14733 :       laa_max = la_max + lca_max
    1274        14733 :       lbb_max = lb_max + lcb_max
    1275        14733 :       li_max = MAX(laa_max, lbb_max)
    1276        14733 :       lj_max = MAX(la_max, lb_max, lcb_max)
    1277        14733 :       lmax = li_max
    1278              : 
    1279        44199 :       ALLOCATE (li_max_all(0:lj_max))
    1280        88398 :       ALLOCATE (Rc(0:lmax, -2*lmax:2*lmax), Rs(0:lmax, -2*lmax:2*lmax))
    1281      2336237 :       Rc = 0._dp
    1282      2336237 :       Rs = 0._dp
    1283        14733 :       mdim(1) = li_max + lj_max + 1
    1284        14733 :       mdim(2) = nsoset(li_max) + 1
    1285        14733 :       mdim(3) = nsoset(lj_max) + 1
    1286        73665 :       ALLOCATE (Waux_mat(mdim(1), mdim(2), mdim(3)))
    1287        73665 :       ALLOCATE (dWaux_mat(3, mdim(1), mdim(2), mdim(3)))
    1288              :       !Waux_mat = 0._dp !.. takes time
    1289              :       !dWaux_mat =0._dp !.. takes time
    1290              : 
    1291              :       !*** Waux_mat (li_max,lj_max) contains elements not needed,
    1292              :       !*** make indixing so that only required ones are computed
    1293              :       !*** li_max_all(j) --> li_max dependent on j
    1294        44199 :       size_int(1, :) = (/laa_max, lb_max/)
    1295        44199 :       size_int(2, :) = (/lbb_max, la_max/)
    1296        44199 :       size_int(3, :) = (/lca_max, lcb_max/)
    1297        44199 :       size_int(4, :) = (/la_max, lb_max/)
    1298              : 
    1299        75471 :       li_max_all(:) = 0
    1300        73665 :       DO isize = 1, 4
    1301        58932 :          i = size_int(isize, 1)
    1302        58932 :          j = size_int(isize, 2)
    1303        58932 :          k = li_max_all(j)
    1304        73665 :          IF (k < i) li_max_all(j) = i
    1305              :       END DO
    1306        14733 :       temp = li_max_all(lj_max)
    1307        75471 :       DO j = lj_max, 0, -1
    1308        75471 :          IF (li_max_all(j) < temp) THEN
    1309        30600 :             li_max_all(j) = temp
    1310              :          ELSE
    1311              :             temp = li_max_all(j)
    1312              :          END IF
    1313              :       END DO
    1314              : 
    1315              :       !*** -rab, since Eq. in Ref. use Ra-Rb, not Rb-Ra
    1316        58932 :       CALL get_real_scaled_solid_harmonic(Rc, Rs, lmax, -rab, rab2)
    1317        14733 :       CALL get_W_matrix(li_max_all, lj_max, lmax, Rc, Rs, Waux_mat)
    1318        14733 :       IF (calculate_forces) THEN
    1319         6140 :          CALL get_dW_matrix(li_max_all, lj_max, Waux_mat, dWaux_mat)
    1320              :       END IF
    1321              : 
    1322        14733 :       DEALLOCATE (Rc, Rs, li_max_all)
    1323              : 
    1324        14733 :    END SUBROUTINE lri_precalc_angular_shg_part
    1325              : 
    1326              : END MODULE generic_shg_integrals
        

Generated by: LCOV version 2.0-1