SHG integrals: performance improvement, new operators

svn-origin-rev: 16962

src/aobasis/ai_contraction_sphi.F

@@ -0,0 +1,215 @@
+!--------------------------------------------------------------------------------------------------!
+!   CP2K: A general program to perform molecular dynamics simulations                              !
+!   Copyright (C) 2000 - 2016  CP2K developers group                                               !
+!--------------------------------------------------------------------------------------------------!
+
+! **************************************************************************************************
+!> \brief Contraction of integrals over primitive Cartesian Gaussians based on the contraction
+!>        matrix sphi which is part of the gto_basis_set_type
+!> \par History
+!>      none
+!> \author Dorothea Golze (05.2016)
+! **************************************************************************************************
+MODULE ai_contraction_sphi
+
+   USE kinds,                           ONLY: dp
+#include "../base/base_uses.f90"
+
+   IMPLICIT NONE
+
+   PRIVATE
+
+   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_contraction_sphi'
+
+   PUBLIC :: ab_contract, abc_contract, abcd_contract
+
+CONTAINS
+
+! **************************************************************************************************
+!> \brief contract overlap integrals (a,b) and transfer to spherical Gaussians
+!> \param abint contracted, normalized integrals of spherical Gaussians
+!> \param sab uncontracted, unnormalized integrals of primitive Cartesian Gaussians
+!> \param sphi_a contraction matrix for center a
+!> \param sphi_b contraction matrix for center b
+!> \param ncoa number of cartesian orbitals on a
+!> \param ncob number of cartesian orbitals on b
+!> \param nsgfa number of spherical Gaussian functions on a
+!> \param nsgfb number of sperical Gaussian functions on b
+! **************************************************************************************************
+   SUBROUTINE ab_contract(abint, sab, sphi_a, sphi_b, ncoa, ncob, nsgfa, nsgfb)
+
+      REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT)      :: abint
+      REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: sab, sphi_a, sphi_b
+      INTEGER, INTENT(IN)                                :: ncoa, ncob, nsgfa, nsgfb
+
+      CHARACTER(LEN=*), PARAMETER :: routineN = 'ab_contract', routineP = moduleN//':'//routineN
+
+      INTEGER                                            :: m1, m2, msphia, msphib, nn
+      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: cpp
+
+      msphia = SIZE(sphi_a, 1)
+      msphib = SIZE(sphi_b, 1)
+
+      m1 = SIZE(sab, 1)
+      m2 = SIZE(sab, 2)
+
+      nn = SIZE(abint, 1)
+
+      ALLOCATE (cpp(nsgfa, m2))
+
+      CALL dgemm("T", "N", nsgfa, m2, ncoa, 1._dp, sphi_a, msphia, sab, m1, 0.0_dp, cpp, nsgfa)
+      CALL dgemm("N", "N", nsgfa, nsgfb, ncob, 1._dp, cpp, nsgfa, sphi_b, msphib, 0.0_dp, &
+                 abint, nn)
+
+      DEALLOCATE (cpp)
+
+   END SUBROUTINE ab_contract
+
+! **************************************************************************************************
+!> \brief contract three-center overlap integrals (a,b,c) and transfer
+!>        to spherical Gaussians
+!> \param abcint contracted, normalized integrals of spherical Gaussians
+!> \param sabc uncontracted, unnormalized integrals of primitive Cartesian Gaussians
+!> \param sphi_a contraction matrix for center a
+!> \param sphi_b contraction matrix for center b
+!> \param sphi_c contraction matrix for center c
+!> \param ncoa number of cartesian orbitals on a
+!> \param ncob number of cartesian orbitals on b
+!> \param ncoc number of cartesian orbitals on c
+!> \param nsgfa number of spherical Gaussian functions on a
+!> \param nsgfb number of spherical Gaussian functions on b
+!> \param nsgfc number of spherical Gaussian functions on c
+! **************************************************************************************************
+   SUBROUTINE abc_contract(abcint, sabc, sphi_a, sphi_b, sphi_c, ncoa, ncob, ncoc, &
+                           nsgfa, nsgfb, nsgfc)
+
+      REAL(KIND=dp), DIMENSION(:, :, :)                  :: abcint, sabc
+      REAL(KIND=dp), DIMENSION(:, :)                     :: sphi_a, sphi_b, sphi_c
+      INTEGER, INTENT(IN)                                :: ncoa, ncob, ncoc, nsgfa, nsgfb, nsgfc
+
+      CHARACTER(LEN=*), PARAMETER :: routineN = 'abc_contract', routineP = moduleN//':'//routineN
+
+      INTEGER                                            :: handle, isgfc, m1, m2, m3, msphia, &
+                                                            msphib, msphic
+      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: temp_ccc, work_cpc
+      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: cpc, cpp
+
+      CALL timeset(routineN, handle)
+
+      msphia = SIZE(sphi_a, 1)
+      msphib = SIZE(sphi_b, 1)
+      msphic = SIZE(sphi_c, 1)
+
+      m1 = SIZE(sabc, 1)
+      m2 = SIZE(sabc, 2)
+      m3 = SIZE(sabc, 3)
+
+      ALLOCATE (cpp(nsgfa, m2, m3), cpc(nsgfa, m2, nsgfc))
+      cpp = 0._dp
+      cpc = 0._dp
+      ALLOCATE (work_cpc(nsgfa, m2), temp_ccc(nsgfa, nsgfb))
+      work_cpc(:, :) = 0._dp
+      temp_ccc(:, :) = 0._dp
+
+      CALL dgemm("T", "N", nsgfa, m2*m3, ncoa, 1._dp, sphi_a, msphia, sabc, m1, 0.0_dp, cpp, nsgfa)
+      CALL dgemm("N", "N", nsgfa*m2, nsgfc, ncoc, 1._dp, cpp, nsgfa*m2, sphi_c, msphic, 0.0_dp, &
+                 cpc, nsgfa*m2)
+
+      DO isgfc = 1, nsgfc
+         work_cpc(:, :) = cpc(:, :, isgfc)
+         CALL dgemm("N", "N", nsgfa, nsgfb, ncob, 1._dp, work_cpc, nsgfa, sphi_b, msphib, &
+                    0.0_dp, temp_ccc, nsgfa)
+         abcint(:, :, isgfc) = temp_ccc(:, :)
+      END DO
+
+      DEALLOCATE (cpp, cpc, work_cpc, temp_ccc)
+
+      CALL timestop(handle)
+
+   END SUBROUTINE abc_contract
+
+! **************************************************************************************************
+!> \brief contract four-center overlap integrals (a,b,c,d) and transfer
+!>        to spherical Gaussians
+!> \param abcdint contracted, normalized integrals of spherical Gaussians
+!> \param sabcd uncontracted, unnormalized integrals of primitive Cartesian Gaussians
+!> \param sphi_a contraction matrix for center a
+!> \param sphi_b contraction matrix for center b
+!> \param sphi_c contraction matrix for center c
+!> \param sphi_d contraction matrix for center d
+!> \param ncoa number of cartesian orbitals on a
+!> \param ncob number of cartesian orbitals on b
+!> \param ncoc number of cartesian orbitals on c
+!> \param ncod number of cartesian orbitals on d
+!> \param nsgfa number of spherical Gaussian functions on a
+!> \param nsgfb number of spherical Gaussian functions on b
+!> \param nsgfc number of spherical Gaussian functions on c
+!> \param nsgfd number of spherical Gaussian functions on d
+! **************************************************************************************************
+   SUBROUTINE abcd_contract(abcdint, sabcd, sphi_a, sphi_b, sphi_c, sphi_d, ncoa, ncob, &
+                            ncoc, ncod, nsgfa, nsgfb, nsgfc, nsgfd)
+
+      REAL(KIND=dp), DIMENSION(:, :, :, :), &
+         INTENT(INOUT)                                   :: abcdint
+      REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN)   :: sabcd
+      REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: sphi_a, sphi_b, sphi_c, sphi_d
+      INTEGER, INTENT(IN)                                :: ncoa, ncob, ncoc, ncod, nsgfa, nsgfb, &
+                                                            nsgfc, nsgfd
+
+      CHARACTER(LEN=*), PARAMETER :: routineN = 'abcd_contract', routineP = moduleN//':'//routineN
+
+      INTEGER                                            :: handle, isgfc, isgfd, m1, m2, m3, m4, &
+                                                            msphia, msphib, msphic, msphid
+      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: temp_cccc, work_cpcc
+      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: temp_cpcc, work_cppc
+      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :)  :: cpcc, cppc, cppp
+
+      CALL timeset(routineN, handle)
+
+      msphia = SIZE(sphi_a, 1)
+      msphib = SIZE(sphi_b, 1)
+      msphic = SIZE(sphi_c, 1)
+      msphid = SIZE(sphi_d, 1)
+
+      m1 = SIZE(sabcd, 1)
+      m2 = SIZE(sabcd, 2)
+      m3 = SIZE(sabcd, 3)
+      m4 = SIZE(sabcd, 4)
+
+      ALLOCATE (cppp(nsgfa, m2, m3, m4), cppc(nsgfa, m2, m3, nsgfd), &
+                cpcc(nsgfa, m2, nsgfc, nsgfd))
+
+      ALLOCATE (work_cppc(nsgfa, m2, m3), temp_cpcc(nsgfa, m2, nsgfc))
+      work_cppc = 0._dp
+      temp_cpcc = 0._dp
+
+      ALLOCATE (work_cpcc(nsgfa, m2), temp_cccc(nsgfa, nsgfb))
+      work_cpcc = 0._dp
+      temp_cccc = 0._dp
+
+      CALL dgemm("T", "N", nsgfa, m2*m3*m4, ncoa, 1._dp, sphi_a, msphia, sabcd, m1, &
+                 0.0_dp, cppp, nsgfa)
+      CALL dgemm("N", "N", nsgfa*m2*m3, nsgfd, ncod, 1._dp, cppp, nsgfa*m2*m3, &
+                 sphi_d, msphid, 0.0_dp, cppc, nsgfa*m2*m3)
+
+      DO isgfd = 1, nsgfd
+         work_cppc(:, :, :) = cppc(:, :, :, isgfd)
+         CALL dgemm("N", "N", nsgfa*m2, nsgfc, ncoc, 1._dp, work_cppc, nsgfa*m2, &
+                    sphi_c, msphic, 0.0_dp, temp_cpcc, nsgfa*m2)
+         cpcc(:, :, :, isgfd) = temp_cpcc(:, :, :)
+         DO isgfc = 1, nsgfc
+            work_cpcc(:, :) = cpcc(:, :, isgfc, isgfd)
+            CALL dgemm("N", "N", nsgfa, nsgfb, ncob, 1._dp, work_cpcc, nsgfa, sphi_b, &
+                       msphib, 0.0_dp, temp_cccc, nsgfa)
+            abcdint(:, :, isgfc, isgfd) = temp_cccc(:, :)
+         END DO
+      END DO
+
+      DEALLOCATE (cpcc, cppc, cppp)
+      DEALLOCATE (work_cpcc, work_cppc, temp_cpcc, temp_cccc)
+
+      CALL timestop(handle)
+
+   END SUBROUTINE abcd_contract
+
+END MODULE ai_contraction_sphi

src/aobasis/ai_derivatives.F

@@ -41,7 +41,7 @@ MODULE ai_derivatives
    PRIVATE
 
 ! *** Public subroutines ***
-   PUBLIC :: adbdr, dabdr
+   PUBLIC :: adbdr, dabdr, dabdr_noscreen
 
 CONTAINS
 
@@ -344,6 +344,126 @@ CONTAINS
       END DO
 
    END SUBROUTINE adbdr
+! **************************************************************************************************
+!> \brief   Calculate the first derivative of an integral block.
+!>          This takes the derivative  with respect to
+!>          the atomic position Ra, i.e. the center of the primitive on the left.
+!>          Difference to routine dabdr: employs no (!!) screening, which is relevant when
+!>          calculating the derivatives of Coulomb integrals
+!> \param la_max ...
+!> \param npgfa ...
+!> \param zeta ...
+!> \param lb_max ...
+!> \param npgfb ...
+!> \param ab ...
+!> \param dabdx ...
+!> \param dabdy ...
+!> \param dabdz ...
+! **************************************************************************************************
+   SUBROUTINE dabdr_noscreen(la_max, npgfa, zeta, lb_max, npgfb, ab, dabdx, dabdy, dabdz)
+      INTEGER, INTENT(IN)                                :: la_max, npgfa
+      REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: zeta
+      INTEGER, INTENT(IN)                                :: lb_max, npgfb
+      REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: ab
+      REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT)        :: dabdx, dabdy, dabdz
+
+      INTEGER                                            :: ax, ay, az, bx, by, bz, coa, coamx, &
+                                                            coamy, coamz, coapx, coapy, coapz, &
+                                                            cob, codb, ipgf, jpgf, la, lb, na, nb, &
+                                                            nda, ndb
+      REAL(KIND=dp)                                      :: fa, fx, fy, fz
+
+!   *** Loop over all pairs of primitive Gaussian-type functions ***
+
+      na = 0
+      nda = 0
+
+      dabdx = 0.0_dp
+      dabdy = 0.0_dp
+      dabdz = 0.0_dp
+
+      DO ipgf = 1, npgfa
+
+         fa = 2.0_dp*zeta(ipgf)
+
+         nb = 0
+         ndb = 0
+
+         DO jpgf = 1, npgfb
+
+            !*** [da/dRi|O|b] = 2*zeta*[a+1i|O|b] - Ni(a)[a-1i|O|b] ***
+
+            DO la = 0, la_max !MAX(0,la_min),la_max
+
+               IF (la == 0) THEN
+
+                  coa = na+1
+                  coapx = nda+2
+                  coapy = nda+3
+                  coapz = nda+4
+
+                  DO lb = 0, lb_max !lb_min,lb_max
+                     DO bx = 0, lb
+                        DO by = 0, lb-bx
+                           bz = lb-bx-by
+                           cob = nb+coset(bx, by, bz)
+                           codb = ndb+coset(bx, by, bz)
+                           dabdx(coa, cob) = fa*ab(coapx, codb)
+                           dabdy(coa, cob) = fa*ab(coapy, codb)
+                           dabdz(coa, cob) = fa*ab(coapz, codb)
+                        END DO
+                     END DO
+                  END DO
+
+               ELSE
+
+                  DO ax = 0, la
+                     DO ay = 0, la-ax
+                        az = la-ax-ay
+
+                        coa = na+coset(ax, ay, az)
+                        coamx = nda+coset(MAX(0, ax-1), ay, az)
+                        coamy = nda+coset(ax, MAX(0, ay-1), az)
+                        coamz = nda+coset(ax, ay, MAX(0, az-1))
+                        coapx = nda+coset(ax+1, ay, az)
+                        coapy = nda+coset(ax, ay+1, az)
+                        coapz = nda+coset(ax, ay, az+1)
+
+                        fx = REAL(ax, dp)
+                        fy = REAL(ay, dp)
+                        fz = REAL(az, dp)
+
+                        DO lb = 0, lb_max !lb_min,lb_max
+                           DO bx = 0, lb
+                              DO by = 0, lb-bx
+                                 bz = lb-bx-by
+                                 cob = nb+coset(bx, by, bz)
+                                 codb = ndb+coset(bx, by, bz)
+                                 dabdx(coa, cob) = fa*ab(coapx, codb)-fx*ab(coamx, codb)
+                                 dabdy(coa, cob) = fa*ab(coapy, codb)-fy*ab(coamy, codb)
+                                 dabdz(coa, cob) = fa*ab(coapz, codb)-fz*ab(coamz, codb)
+                              END DO
+                           END DO
+                        END DO
+
+                     END DO
+                  END DO
+
+               END IF
+
+            END DO
+
+            nb = nb+ncoset(lb_max)
+            ndb = ndb+ncoset(lb_max+1)
+
+         END DO
+
+         na = na+ncoset(la_max)
+         nda = nda+ncoset(la_max+1)
+
+      END DO
+
+   END SUBROUTINE dabdr_noscreen
 
 END MODULE ai_derivatives
 

src/aobasis/ai_operators_r12.F

@@ -0,0 +1,620 @@
+!--------------------------------------------------------------------------------------------------!
+!   CP2K: A general program to perform molecular dynamics simulations                              !
+!   Copyright (C) 2000 - 2016  CP2K developers group                                               !
+!--------------------------------------------------------------------------------------------------!
+
+! **************************************************************************************************
+!> \brief Calculation of integrals over Cartesian Gaussian-type functions for different r12
+!>        operators: 1/r12, erf(omega*r12/r12), erfc(omega*r12/r12), exp(-omega*r12^2)/r12 and
+!>                   exp(-omega*r12^2)
+!> \par Literature
+!>      S. Obara and A. Saika, J. Chem. Phys. 84, 3963 (1986)
+!>      R. Ahlrichs, PCCP, 8, 3072 (2006)
+!> \par History
+!>      none
+!> \par Parameters
+!>       - ax,ay,az    : Angular momentum index numbers of orbital a.
+!>       - cx,cy,cz    : Angular momentum index numbers of orbital c.
+!>       - coset       : Cartesian orbital set pointer.
+!>       - dac         : Distance between the atomic centers a and c.
+!>       - l{a,c}      : Angular momentum quantum number of shell a or c.
+!>       - l{a,c}_max  : Maximum angular momentum quantum number of shell a or c.
+!>       - l{a,c}_min  : Minimum angular momentum quantum number of shell a or c.
+!>       - ncoset      : Number of orbitals in a Cartesian orbital set.
+!>       - npgf{a,c}   : Degree of contraction of shell a or c.
+!>       - rac         : Distance vector between the atomic centers a and c.
+!>       - rac2        : Square of the distance between the atomic centers a and c.
+!>       - zet{a,c}    : Exponents of the Gaussian-type functions a or c.
+!>       - zetp        : Reciprocal of the sum of the exponents of orbital a and b.
+!>       - zetw        : Reciprocal of the sum of the exponents of orbital a and c.
+!>       - omega       : Parameter in the operator
+!> \author Dorothea Golze (05.2016)
+! **************************************************************************************************
+MODULE ai_operators_r12
+
+   USE gamma,                           ONLY: fgamma=>fgamma_0
+   USE kinds,                           ONLY: dp
+   USE mathconstants,                   ONLY: fac,&
+                                              pi
+   USE orbital_pointers,                ONLY: coset,&
+                                              ncoset
+#include "../base/base_uses.f90"
+
+   IMPLICIT NONE
+   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_operators_r12'
+   PRIVATE
+
+   ! *** Public subroutines ***
+
+   PUBLIC :: operator2, cps_coulomb2, cps_verf2, cps_verfc2, cps_vgauss2, cps_gauss2
+
+#include "./ai_operators_r12.h"
+
+CONTAINS
+
+! **************************************************************************************************
+!> \brief Calculation of the primitive two-center integrals over Cartesian Gaussian-type
+!>        functions for different r12 operators.
+!> \param cps_operator2 procedure pointer for the respective operator. The integrals evaluation
+!>        differs only in the evaluation of the cartesian primitive s (cps) integrals [s|O(r12)|s]
+!>        and auxiliary integrals [s|O(r12)|s]^n. This pointer selects the correct routine.
+!> \param la_max ...
+!> \param npgfa ...
+!> \param zeta ...
+!> \param la_min ...
+!> \param lc_max ...
+!> \param npgfc ...
+!> \param zetc ...
+!> \param lc_min ...
+!> \param omega ...
+!> \param rac ...
+!> \param rac2 ...
+!> \param vac matrix storing the integrals
+!> \param v temporary work array
+!> \param maxder maximal derivative
+!> \param vac_plus matrix storing the integrals for highler l-quantum numbers; used to
+!>        construct the derivatives
+! **************************************************************************************************
+
+   SUBROUTINE operator2(cps_operator2, la_max, npgfa, zeta, la_min, lc_max, npgfc, zetc, lc_min, &
+                        omega, rac, rac2, vac, v, maxder, vac_plus)
+      PROCEDURE(ab_sint_os), POINTER                     :: cps_operator2
+      INTEGER, INTENT(IN)                                :: la_max, npgfa
+      REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: zeta
+      INTEGER, INTENT(IN)                                :: la_min, lc_max, npgfc
+      REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: zetc
+      INTEGER, INTENT(IN)                                :: lc_min
+      REAL(KIND=dp), INTENT(IN)                          :: omega
+      REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rac
+      REAL(KIND=dp), INTENT(IN)                          :: rac2
+      REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT)      :: vac
+      REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT)    :: v
+      INTEGER, INTENT(IN), OPTIONAL                      :: maxder
+      REAL(KIND=dp), DIMENSION(:, :), OPTIONAL           :: vac_plus
+
+      CHARACTER(len=*), PARAMETER :: routineN = 'operator2', routineP = moduleN//':'//routineN
+
+      INTEGER                                            :: ax, ay, az, coc, cocx, cocy, cocz, cx, &
+                                                            cy, cz, i, ipgf, j, jpgf, la, lc, &
+                                                            maxder_local, n, na, nap, nc, ncp, nmax
+      REAL(KIND=dp)                                      :: dac, f1, f2, f3, f4, f5, f6, fcx, &
+                                                            fcy, fcz, rho, zetp, zetq, zetw
+      REAL(KIND=dp), DIMENSION(3)                        :: raw, rcw
+
+      v = 0.0_dp
+
+      maxder_local = 0
+      IF (PRESENT(maxder)) THEN
+         maxder_local = maxder
+         vac_plus = 0.0_dp
+      END IF
+
+      nmax = la_max+lc_max+1
+
+      ! *** Calculate the distance of the centers a and c ***
+
+      dac = SQRT(rac2)
+
+      ! *** Loop over all pairs of primitive Gaussian-type functions ***
+
+      na = 0
+      nap = 0
+
+      DO ipgf = 1, npgfa
+
+         nc = 0
+         ncp = 0
+
+         DO jpgf = 1, npgfc
+
+            ! *** Calculate some prefactors ***
+
+            zetp = 1.0_dp/zeta(ipgf)
+            zetq = 1.0_dp/zetc(jpgf)
+            zetw = 1.0_dp/(zeta(ipgf)+zetc(jpgf))
+
+            rho = zeta(ipgf)*zetc(jpgf)*zetw
+
+            ! *** Calculate the basic two-center integrals [s||s]{n} ***
+
+            CALL cps_operator2(v, nmax, zetp, zetq, zetw, rho, rac2, omega)
+
+            ! *** Vertical recurrence steps: [s||s] -> [s||c] ***
+
+            IF (lc_max > 0) THEN
+
+               f1 = 0.5_dp*zetq
+               f2 = -rho*zetq
+
+               rcw(:) = -zeta(ipgf)*zetw*rac(:)
+
+               ! *** [s||p]{n} = (Wi - Ci)*[s||s]{n+1}  (i = x,y,z) ***
+
+               DO n = 1, nmax-1
+                  v(1, 2, n) = rcw(1)*v(1, 1, n+1)
+                  v(1, 3, n) = rcw(2)*v(1, 1, n+1)
+                  v(1, 4, n) = rcw(3)*v(1, 1, n+1)
+               END DO
+
+               ! ** [s||c]{n} = (Wi - Ci)*[s||c-1i]{n+1} +     ***
+               ! **             f1*Ni(c-1i)*(   [s||c-2i]{n} + ***
+               ! **                          f2*[s||c-2i]{n+1} ***
+
+               DO lc = 2, lc_max
+
+                  DO n = 1, nmax-lc
+
+                     ! **** Increase the angular momentum component z of c ***
+
+                     v(1, coset(0, 0, lc), n) = &
+                        rcw(3)*v(1, coset(0, 0, lc-1), n+1)+ &
+                        f1*REAL(lc-1, dp)*(v(1, coset(0, 0, lc-2), n)+ &
+                                           f2*v(1, coset(0, 0, lc-2), n+1))
+
+                     ! *** Increase the angular momentum component y of c ***
+
+                     cz = lc-1
+                     v(1, coset(0, 1, cz), n) = rcw(2)*v(1, coset(0, 0, cz), n+1)
+
+                     DO cy = 2, lc
+                        cz = lc-cy
+                        v(1, coset(0, cy, cz), n) = &
+                           rcw(2)*v(1, coset(0, cy-1, cz), n+1)+ &
+                           f1*REAL(cy-1, dp)*(v(1, coset(0, cy-2, cz), n)+ &
+                                              f2*v(1, coset(0, cy-2, cz), n+1))
+                     END DO
+
+                     ! *** Increase the angular momentum component x of c ***
+
+                     DO cy = 0, lc-1
+                        cz = lc-1-cy
+                        v(1, coset(1, cy, cz), n) = rcw(1)*v(1, coset(0, cy, cz), n+1)
+                     END DO
+
+                     DO cx = 2, lc
+                        f6 = f1*REAL(cx-1, dp)
+                        DO cy = 0, lc-cx
+                           cz = lc-cx-cy
+                           v(1, coset(cx, cy, cz), n) = &
+                              rcw(1)*v(1, coset(cx-1, cy, cz), n+1)+ &
+                              f6*(v(1, coset(cx-2, cy, cz), n)+ &
+                                  f2*v(1, coset(cx-2, cy, cz), n+1))
+                        END DO
+                     END DO
+
+                  END DO
+
+               END DO
+
+            END IF
+
+            ! *** Vertical recurrence steps: [s||c] -> [a||c] ***
+
+            IF (la_max > 0) THEN
+
+               f3 = 0.5_dp*zetp
+               f4 = -rho*zetp
+               f5 = 0.5_dp*zetw
+
+               raw(:) = zetc(jpgf)*zetw*rac(:)
+
+               ! *** [p||s]{n} = (Wi - Ai)*[s||s]{n+1}  (i = x,y,z) ***
+
+               DO n = 1, nmax-1
+                  v(2, 1, n) = raw(1)*v(1, 1, n+1)
+                  v(3, 1, n) = raw(2)*v(1, 1, n+1)
+                  v(4, 1, n) = raw(3)*v(1, 1, n+1)
+               END DO
+
+               ! *** [a||s]{n} = (Wi - Ai)*[a-1i||s]{n+1} +      ***
+               ! ***             f3*Ni(a-1i)*(   [a-2i||s]{n} +  ***
+               ! ***                          f4*[a-2i||s]{n+1}) ***
+
+               DO la = 2, la_max
+
+                  DO n = 1, nmax-la
+
+                     ! *** Increase the angular momentum component z of a ***
+
+                     v(coset(0, 0, la), 1, n) = &
+                        raw(3)*v(coset(0, 0, la-1), 1, n+1)+ &
+                        f3*REAL(la-1, dp)*(v(coset(0, 0, la-2), 1, n)+ &
+                                           f4*v(coset(0, 0, la-2), 1, n+1))
+
+                     ! *** Increase the angular momentum component y of a ***
+
+                     az = la-1
+                     v(coset(0, 1, az), 1, n) = raw(2)*v(coset(0, 0, az), 1, n+1)
+
+                     DO ay = 2, la
+                        az = la-ay
+                        v(coset(0, ay, az), 1, n) = &
+                           raw(2)*v(coset(0, ay-1, az), 1, n+1)+ &
+                           f3*REAL(ay-1, dp)*(v(coset(0, ay-2, az), 1, n)+ &
+                                              f4*v(coset(0, ay-2, az), 1, n+1))
+                     END DO
+
+                     ! *** Increase the angular momentum component x of a ***
+
+                     DO ay = 0, la-1
+                        az = la-1-ay
+                        v(coset(1, ay, az), 1, n) = raw(1)*v(coset(0, ay, az), 1, n+1)
+                     END DO
+
+                     DO ax = 2, la
+                        f6 = f3*REAL(ax-1, dp)
+                        DO ay = 0, la-ax
+                           az = la-ax-ay
+                           v(coset(ax, ay, az), 1, n) = &
+                              raw(1)*v(coset(ax-1, ay, az), 1, n+1)+ &
+                              f6*(v(coset(ax-2, ay, az), 1, n)+ &
+                                  f4*v(coset(ax-2, ay, az), 1, n+1))
+                        END DO
+                     END DO
+
+                  END DO
+
+               END DO
+
+               DO lc = 1, lc_max
+
+                  DO cx = 0, lc
+                     DO cy = 0, lc-cx
+                        cz = lc-cx-cy
+
+                        coc = coset(cx, cy, cz)
+                        cocx = coset(MAX(0, cx-1), cy, cz)
+                        cocy = coset(cx, MAX(0, cy-1), cz)
+                        cocz = coset(cx, cy, MAX(0, cz-1))
+
+                        fcx = f5*REAL(cx, dp)
+                        fcy = f5*REAL(cy, dp)
+                        fcz = f5*REAL(cz, dp)
+
+                        ! *** [p||c]{n} = (Wi - Ai)*[s||c]{n+1} + ***
+                        ! ***             f5*Ni(c)*[s||c-1i]{n+1} ***
+
+                        DO n = 1, nmax-1-lc
+                           v(2, coc, n) = raw(1)*v(1, coc, n+1)+fcx*v(1, cocx, n+1)
+                           v(3, coc, n) = raw(2)*v(1, coc, n+1)+fcy*v(1, cocy, n+1)
+                           v(4, coc, n) = raw(3)*v(1, coc, n+1)+fcz*v(1, cocz, n+1)
+                        END DO
+
+                        ! *** [a||c]{n} = (Wi - Ai)*[a-1i||c]{n+1} +        ***
+                        ! ***             f3*Ni(a-1i)*(   [a-2i||c]{n} +    ***
+                        ! ***                          f4*[a-2i||c]{n+1}) + ***
+                        ! ***             f5*Ni(c)*[a-1i||c-1i]{n+1}        ***
+
+                        DO la = 2, la_max
+
+                           DO n = 1, nmax-la-lc
+
+                              ! *** Increase the angular momentum component z of a ***
+
+                              v(coset(0, 0, la), coc, n) = &
+                                 raw(3)*v(coset(0, 0, la-1), coc, n+1)+ &
+                                 f3*REAL(la-1, dp)*(v(coset(0, 0, la-2), coc, n)+ &
+                                                    f4*v(coset(0, 0, la-2), coc, n+1))+ &
+                                 fcz*v(coset(0, 0, la-1), cocz, n+1)
+
+                              ! *** Increase the angular momentum component y of a ***
+
+                              az = la-1
+                              v(coset(0, 1, az), coc, n) = &
+                                 raw(2)*v(coset(0, 0, az), coc, n+1)+ &
+                                 fcy*v(coset(0, 0, az), cocy, n+1)
+
+                              DO ay = 2, la
+                                 az = la-ay
+                                 v(coset(0, ay, az), coc, n) = &
+                                    raw(2)*v(coset(0, ay-1, az), coc, n+1)+ &
+                                    f3*REAL(ay-1, dp)*(v(coset(0, ay-2, az), coc, n)+ &
+                                                       f4*v(coset(0, ay-2, az), coc, n+1))+ &
+                                    fcy*v(coset(0, ay-1, az), cocy, n+1)
+                              END DO
+
+                              ! *** Increase the angular momentum component x of a ***
+
+                              DO ay = 0, la-1
+                                 az = la-1-ay
+                                 v(coset(1, ay, az), coc, n) = &
+                                    raw(1)*v(coset(0, ay, az), coc, n+1)+ &
+                                    fcx*v(coset(0, ay, az), cocx, n+1)
+                              END DO
+
+                              DO ax = 2, la
+                                 f6 = f3*REAL(ax-1, dp)
+                                 DO ay = 0, la-ax
+                                    az = la-ax-ay
+                                    v(coset(ax, ay, az), coc, n) = &
+                                       raw(1)*v(coset(ax-1, ay, az), coc, n+1)+ &
+                                       f6*(v(coset(ax-2, ay, az), coc, n)+ &
+                                           f4*v(coset(ax-2, ay, az), coc, n+1))+ &
+                                       fcx*v(coset(ax-1, ay, az), cocx, n+1)
+                                 END DO
+                              END DO
+
+                           END DO
+
+                        END DO
+
+                     END DO
+                  END DO
+
+               END DO
+
+            END IF
+
+            DO j = ncoset(lc_min-1)+1, ncoset(lc_max-maxder_local)
+               DO i = ncoset(la_min-1)+1, ncoset(la_max-maxder_local)
+                  vac(na+i, nc+j) = v(i, j, 1)
+               END DO
+            END DO
+
+            IF (PRESENT(maxder)) THEN
+               DO j = 1, ncoset(lc_max)
+                  DO i = 1, ncoset(la_max)
+                     vac_plus(nap+i, ncp+j) = v(i, j, 1)
+                  END DO
+               END DO
+            END IF
+
+            nc = nc+ncoset(lc_max-maxder_local)
+            ncp = ncp+ncoset(lc_max)
+
+         END DO
+
+         na = na+ncoset(la_max-maxder_local)
+         nap = nap+ncoset(la_max)
+
+      END DO
+
+   END SUBROUTINE operator2
+
+! **************************************************************************************************
+!> \brief Calculation of Coulomb integrals for s-function, i.e, [s|1/r12|s], and the auxiliary
+!>        integrals [s|1/r12|s]^n
+!> \param v matrix storing the integrals
+!> \param nmax maximal n in the auxiliary integrals [s|1/r12|s]^n
+!> \param zetp = 1/zeta
+!> \param zetq = 1/zetc
+!> \param zetw = 1/(zeta+zetc)
+!> \param rho = zeta*zetc*zetw
+!> \param rac2 square distance between center A and C, |Ra-Rc|^2
+!> \param omega this parameter is actually not used, but included for the sake of the abstract
+!>        interface
+! **************************************************************************************************
+   SUBROUTINE cps_coulomb2(v, nmax, zetp, zetq, zetw, rho, rac2, omega)
+      REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT)   :: v
+      INTEGER, INTENT(IN)                                :: nmax
+      REAL(KIND=dp), INTENT(IN)                          :: zetp, zetq, zetw, rho, rac2, omega
+
+      CHARACTER(len=*), PARAMETER :: routineN = 'cps_coulomb2', routineP = moduleN//':'//routineN
+
+      INTEGER                                            :: n
+      REAL(KIND=dp)                                      :: dummy, f0, t
+      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: f
+
+      ALLOCATE (f(0:nmax))
+      dummy = omega
+      f0 = 2.0_dp*SQRT(pi**5*zetw)*zetp*zetq
+
+      ! *** Calculate the incomplete Gamma/Boys function ***
+
+      t = rho*rac2
+      CALL fgamma(nmax-1, t, f)
+
+      ! *** Calculate the basic two-center integrals [s||s]{n} ***
+
+      DO n = 1, nmax
+         v(1, 1, n) = f0*f(n-1)
+      END DO
+
+      DEALLOCATE (f)
+   END SUBROUTINE cps_coulomb2
+
+! **************************************************************************************************
+!> \brief Calculation of verf integrals for s-function, i.e, [s|erf(omega*r12)/r12|s], and the
+!>        auxiliary integrals [s|erf(omega*r12)/r12|s]^n
+!> \param v matrix storing the integrals
+!> \param nmax maximal n in the auxiliary integrals [s|erf(omega*r12)/r12|s]^n
+!> \param zetp = 1/zeta
+!> \param zetq = 1/zetc
+!> \param zetw = 1/(zeta+zetc)
+!> \param rho = zeta*zetc*zetw
+!> \param rac2 square distance between center A and C, |Ra-Rc|^2
+!> \param omega parameter in the operator
+! **************************************************************************************************
+   SUBROUTINE cps_verf2(v, nmax, zetp, zetq, zetw, rho, rac2, omega)
+      REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT)   :: v
+      INTEGER, INTENT(IN)                                :: nmax
+      REAL(KIND=dp), INTENT(IN)                          :: zetp, zetq, zetw, rho, rac2, omega
+
+      CHARACTER(len=*), PARAMETER :: routineN = 'cps_verf2', routineP = moduleN//':'//routineN
+
+      INTEGER                                            :: n
+      REAL(KIND=dp)                                      :: arg, comega, f0, t
+      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: f
+
+      ALLOCATE (f(0:nmax))
+      comega = omega**2/(omega**2+rho)
+      f0 = 2.0_dp*SQRT(pi**5*zetw*comega)*zetp*zetq
+
+      ! *** Calculate the incomplete Gamma/Boys function ***
+
+      t = rho*rac2
+      arg = comega*t
+      CALL fgamma(nmax-1, arg, f)
+
+      ! *** Calculate the basic two-center integrals [s||s]{n} ***
+
+      DO n = 1, nmax
+         v(1, 1, n) = f0*f(n-1)*comega**(n-1)
+      END DO
+
+      DEALLOCATE (f)
+
+   END SUBROUTINE cps_verf2
+
+! **************************************************************************************************
+!> \brief Calculation of verfc integrals for s-function, i.e, [s|erfc(omega*r12)/r12|s], and
+!>        the auxiliary integrals [s|erfc(omega*r12)/r12|s]^n
+!> \param v matrix storing the integrals
+!> \param nmax maximal n in the auxiliary integrals [s|erfc(omega*r12)/r12|s]^n
+!> \param zetp = 1/zeta
+!> \param zetq = 1/zetc
+!> \param zetw = 1/(zeta+zetc)
+!> \param rho = zeta*zetc*zetw
+!> \param rac2 square distance between center A and C, |Ra-Rc|^2
+!> \param omega parameter in the operator
+! **************************************************************************************************
+   SUBROUTINE cps_verfc2(v, nmax, zetp, zetq, zetw, rho, rac2, omega)
+      REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT)   :: v
+      INTEGER, INTENT(IN)                                :: nmax
+      REAL(KIND=dp), INTENT(IN)                          :: zetp, zetq, zetw, rho, rac2, omega
+
+      CHARACTER(len=*), PARAMETER :: routineN = 'cps_verfc2', routineP = moduleN//':'//routineN
+
+      INTEGER                                            :: n
+      REAL(KIND=dp)                                      :: argerf, comega, f0, t
+      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: fv, fverf
+
+      ALLOCATE (fv(0:nmax), fverf(0:nmax))
+      comega = omega**2/(omega**2+rho)
+      f0 = 2.0_dp*SQRT(pi**5*zetw)*zetp*zetq
+
+      ! *** Calculate the incomplete Gamma/Boys function ***
+
+      t = rho*rac2
+      argerf = comega*t
+
+      CALL fgamma(nmax-1, t, fv)
+      CALL fgamma(nmax-1, argerf, fverf)
+
+      ! *** Calculate the basic two-center integrals [s||s]{n} ***
+
+      DO n = 1, nmax
+         v(1, 1, n) = f0*(fv(n-1)-SQRT(comega)*comega**(n-1)*fverf(n-1))
+      END DO
+
+      DEALLOCATE (fv, fverf)
+
+   END SUBROUTINE cps_verfc2
+
+! **************************************************************************************************
+!> \brief Calculation of vgauss integrals for s-function, i.e, [s|exp(-omega*r12^2)/r12|s], and
+!>        the auxiliary integrals [s|exp(-omega*r12^2)/r12|s]
+!> \param v matrix storing the integrals
+!> \param nmax maximal n in the auxiliary integrals [s|exp(-omega*r12^2)/r12|s]
+!> \param zetp = 1/zeta
+!> \param zetq = 1/zetc
+!> \param zetw = 1/(zeta+zetc)
+!> \param rho = zeta*zetc*zetw
+!> \param rac2 square distance between center A and C, |Ra-Rc|^2
+!> \param omega parameter in the operator
+! **************************************************************************************************
+   SUBROUTINE cps_vgauss2(v, nmax, zetp, zetq, zetw, rho, rac2, omega)
+      REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT)   :: v
+      INTEGER, INTENT(IN)                                :: nmax
+      REAL(KIND=dp), INTENT(IN)                          :: zetp, zetq, zetw, rho, rac2, omega
+
+      CHARACTER(len=*), PARAMETER :: routineN = 'cps_vgauss2', routineP = moduleN//':'//routineN
+
+      INTEGER                                            :: j, n
+      REAL(KIND=dp)                                      :: arg, dummy, eta, expT, f0, fsign, t, tau
+      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: f
+
+      ALLOCATE (f(0:nmax))
+
+      dummy = zetp
+      dummy = zetq
+      eta = rho/(rho+omega)
+      tau = omega/(rho+omega)
+
+      ! *** Calculate the incomplete Gamma/Boys function ***
+
+      t = rho*rac2
+      arg = eta*t
+
+      CALL fgamma(nmax-1, arg, f)
+
+      expT = EXP(-omega/(omega+rho)*t)
+      f0 = 2.0_dp*SQRT(pi**5*zetw**3)/(rho+omega)*expT
+
+      ! *** Calculate the basic two-center integrals [s||s]{n} ***
+      v(1, 1, 1:nmax) = 0.0_dp
+      DO n = 1, nmax
+         fsign = (-1.0_dp)**(n-1)
+         DO j = 0, n-1
+            v(1, 1, n) = v(1, 1, n)+f0*fsign* &
+                         fac(n-1)/fac(n-j-1)/fac(j)*(-tau)**(n-j-1)*(-eta)**j*f(j)
+         ENDDO
+      ENDDO
+
+      DEALLOCATE (f)
+
+   END SUBROUTINE cps_vgauss2
+
+! **************************************************************************************************
+!> \brief Calculation of gauss integrals for s-function, i.e, [s|exp(-omega*r12^2)|s], and
+!>        the auxiliary integrals [s|exp(-omega*r12^2)|s]
+!> \param v matrix storing the integrals
+!> \param nmax maximal n in the auxiliary integrals [s|exp(-omega*r12^2)|s]
+!> \param zetp = 1/zeta
+!> \param zetq = 1/zetc
+!> \param zetw = 1/(zeta+zetc)
+!> \param rho = zeta*zetc*zetw
+!> \param rac2 square distance between center A and C, |Ra-Rc|^2
+!> \param omega parameter in the operator
+! **************************************************************************************************
+   SUBROUTINE cps_gauss2(v, nmax, zetp, zetq, zetw, rho, rac2, omega)
+      REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT)   :: v
+      INTEGER, INTENT(IN)                                :: nmax
+      REAL(KIND=dp), INTENT(IN)                          :: zetp, zetq, zetw, rho, rac2, omega
+
+      CHARACTER(len=*), PARAMETER :: routineN = 'cps_gauss2', routineP = moduleN//':'//routineN
+
+      INTEGER                                            :: n
+      REAL(KIND=dp)                                      :: dummy, expT, f0, t, tau
+      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: f
+
+      ALLOCATE (f(0:nmax))
+
+      dummy = zetp
+      dummy = zetq
+      tau = omega/(rho+omega)
+      t = rho*rac2
+      expT = EXP(-tau*t)
+      f0 = pi**3*SQRT(zetw**3/(rho+omega)**3)*expT
+
+      ! *** Calculate the basic two-center integrals [s||s]{n} ***
+
+      DO n = 1, nmax
+         v(1, 1, n) = f0*tau**(n-1)
+      END DO
+
+      DEALLOCATE (f)
+
+   END SUBROUTINE cps_gauss2
+
+END MODULE ai_operators_r12

src/aobasis/ai_operators_r12.h

@@ -0,0 +1,17 @@
+!--------------------------------------------------------------------------------------------------!
+!   CP2K: A general program to perform molecular dynamics simulations                              !
+!   Copyright (C) 2000 - 2016  CP2K developers group                                               !
+!--------------------------------------------------------------------------------------------------!
+! **************************************************************************************************
+!> \brief Interface for the calculation of integrals over s-functions and the s-type auxiliary
+!>        integrals using the Obara-Saika (OS) scheme
+! **************************************************************************************************
+   ABSTRACT INTERFACE
+      SUBROUTINE ab_sint_os(v, nmax, zetp, zetq, zetw, rho, rac2, omega)
+       USE kinds,                           ONLY: dp
+      REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT)   :: v
+      INTEGER, INTENT(IN)                                :: nmax
+      REAL(KIND=dp), INTENT(IN)                          :: zetp, zetq, zetw, rho, rac2, omega
+
+      END SUBROUTINE ab_sint_os
+   END INTERFACE

src/aobasis/basis_set_types.F

@@ -1608,7 +1608,7 @@ CONTAINS
       CALL reallocate(gto_basis_set%norm_cgf, 1, ncgf)
 
       basis_type = ""
-      CALL section_vals_val_get(basis_section, "_SECTION_PARAMETERS_", i_rep_val=irep, c_val=basis_type)
+      CALL section_vals_val_get(basis_section, "_SECTION_PARAMETERS_", i_rep_section=irep, c_val=basis_type)
       IF (basis_type == "Orbital") basis_type = "ORB"
 
       NULLIFY (list, val)

src/aobasis/orbital_pointers.F

@@ -22,8 +22,6 @@ MODULE orbital_pointers
 !         number l.
 ! nsoset: Number of spherical orbitals up to the angular momentum quantum
 !         number l.
-! nsoset_pm: Number of spherical orbitals up to the angular momentum quantum
-!         number l for positive m only, i.e. m >= 0
 
    USE kinds,                           ONLY: int_size
    USE termination,                     ONLY: stop_memory
@@ -39,7 +37,7 @@ MODULE orbital_pointers
 
    INTEGER, SAVE :: current_maxl = -1
 
-   INTEGER, DIMENSION(:), ALLOCATABLE     :: nco, ncoset, nso, nsoset, nsoset_pm
+   INTEGER, DIMENSION(:), ALLOCATABLE     :: nco, ncoset, nso, nsoset
    INTEGER, DIMENSION(:, :), ALLOCATABLE   :: indco, indso, indso_inv
    INTEGER, DIMENSION(:, :), ALLOCATABLE   :: so, soset
    INTEGER, DIMENSION(:, :, :), ALLOCATABLE :: co, coset
@@ -61,7 +59,6 @@ MODULE orbital_pointers
              ncoset, &
              nso, &
              nsoset, &
-             nsoset_pm, &
              soset
 
 CONTAINS
@@ -192,18 +189,10 @@ CONTAINS
          CALL stop_memory(routineN, moduleN, __LINE__, &
                           "nsoset", (maxl+2)*int_size)
       END IF
-      ALLOCATE (nsoset_pm(-1:maxl), STAT=istat)
-      IF (istat /= 0) THEN
-         CALL stop_memory(routineN, moduleN, __LINE__, &
-                          "nsoset_pm", (maxl+2)*int_size)
-      END IF
-
       nsoset(-1) = 0
-      nsoset_pm(-1) = 0
 
       DO l = 0, maxl
          nsoset(l) = nsoset(l-1)+nso(l)
-         nsoset_pm(l) = nsoset_pm(l-1)+l+1
       END DO
 
       ALLOCATE (indso(2, nsoset(maxl)), STAT=istat)
@@ -281,8 +270,6 @@ CONTAINS
 
          DEALLOCATE (nsoset)
 
-         DEALLOCATE (nsoset_pm)
-
          DEALLOCATE (so)
 
          DEALLOCATE (soset)

src/debug_os_integrals.F

@@ -4,13 +4,12 @@
 !--------------------------------------------------------------------------------------------------!
 
 ! **************************************************************************************************
-!> \brief Debugs the  integral matrices for the LRIGPW method
-!>        lri : local resolution of the identity
+!> \brief Debugs Obara-Saika integral matrices
 !> \par History
 !>      created [07.2014]
 !> \authors Dorothea Golze
 ! **************************************************************************************************
-MODULE lri_debug_integrals
+MODULE debug_os_integrals
 
    USE ai_overlap,                      ONLY: overlap
    USE ai_overlap3,                     ONLY: overlap3
@@ -31,7 +30,7 @@ MODULE lri_debug_integrals
 
 ! **************************************************************************************************
 
-   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'lri_debug_integrals'
+   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'debug_os_integrals'
 
    PUBLIC :: overlap_ab_test, overlap_abc_test, overlap_aabb_test
 
@@ -247,11 +246,6 @@ CONTAINS
                     lb_max_set=lb_max, npgfb=1, zetb=xb_work, rpgfb=rpgfb, lb_min_set=lb_min, &
                     lc_max_set=lc_max, npgfc=1, zetc=xc_work, rpgfc=rpgfc, lc_min_set=lc_min, &
                     rab=rab, dab=dab, rac=rac, dac=dac, rbc=rbc, dbc=dbc, sabc=sabc)
-      !sabc(:,:,:)=0._dp
-      !CALL overlap3_old(la_max=la_max,npgfa=1,zeta=xa_work,rpgfa=rpgfa,la_min=la_min,&
-      !              lb_max=lb_max,npgfb=1,zetb=xb_work,rpgfb=rpgfb,lb_min=lb_min,&
-      !              lc_max=lc_max,npgfc=1,zetc=xc,rpgfc=20._dp,lc_min=lc_min,&
-      !              rab=rab,dab=dab,rac=rac,dac=dac,rbc=rbc,dbc=dbc,sabc=sabc)
 
       !---------------------------------------
 
@@ -624,4 +618,4 @@ CONTAINS
 
    END SUBROUTINE overlap_aabb_test
 
-END MODULE lri_debug_integrals
+END MODULE debug_os_integrals

src/generic_os_integrals.F

@@ -4,23 +4,39 @@
 !--------------------------------------------------------------------------------------------------!
 
 ! **************************************************************************************************
-!> \brief Calculates integral matrices for LRIGPW method using the Obara-Saika
-!>        (OS) integral scheme
-!>        lri : local resolution of the identity
+!> \brief Calculation of contracted, spherical Gaussian integrals using the (OS) integral
+!>        scheme. Routines for the following two-center integrals:
+!>        i)  (a|O(r12)|b) where O(r12) is the overlap, coulomb operator etc.
+!>        ii) (aba) and (abb) s-overlaps
 !> \par History
 !>      created [06.2015]
 !> \author Dorothea Golze
 ! **************************************************************************************************
-MODULE lri_os_integrals
+MODULE generic_os_integrals
+   USE ai_contraction_sphi,             ONLY: ab_contract,&
+                                              abc_contract,&
+                                              abcd_contract
+   USE ai_derivatives,                  ONLY: dabdr_noscreen
+   USE ai_operators_r12,                ONLY: cps_coulomb2,&
+                                              cps_gauss2,&
+                                              cps_verf2,&
+                                              cps_verfc2,&
+                                              cps_vgauss2,&
+                                              operator2
    USE ai_overlap,                      ONLY: overlap
    USE ai_overlap3,                     ONLY: overlap3
    USE ai_overlap_aabb,                 ONLY: overlap_aabb
    USE basis_set_types,                 ONLY: get_gto_basis_set,&
                                               gto_basis_set_type
-   USE kinds,                           ONLY: dp
-   USE lri_debug_integrals,             ONLY: overlap_aabb_test,&
+   USE debug_os_integrals,              ONLY: overlap_aabb_test,&
                                               overlap_ab_test,&
                                               overlap_abc_test
+   USE input_constants,                 ONLY: operator_coulomb,&
+                                              operator_gauss,&
+                                              operator_verf,&
+                                              operator_verfc,&
+                                              operator_vgauss
+   USE kinds,                           ONLY: dp
    USE orbital_pointers,                ONLY: ncoset
 #include "./base/base_uses.f90"
 
@@ -30,29 +46,220 @@ MODULE lri_os_integrals
 
 ! **************************************************************************************************
 
-   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'lri_os_integrals'
+   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'generic_os_integrals'
 
-   PUBLIC :: lri_int_ab, lri_int_aba, lri_int_abb, lri_int_aabb
+   PUBLIC :: two_center_os_integrals, int_overlap_ab_os, int_overlap_aba_os, &
+             int_overlap_abb_os, int_overlap_aabb_os
+
+#include "./aobasis/ai_operators_r12.h"
 
 ! **************************************************************************************************
 
 CONTAINS
 
+! **************************************************************************************************
+!> \brief Calcululates the two-center integrals of the type (a|O(r12)|b) using the OS scheme
+!> \param r12_operator the integral operator, which depends on r12=|r1-r2|
+!> \param vab integral matrix of spherical contracted Gaussian functions
+!> \param dvab derivative of the integrals
+!> \param rab distance vector between center A and B
+!> \param fba basis at center A
+!> \param fbb basis at center B
+!> \param omega parameter in the operator
+!> \param calculate_forces ...
+! **************************************************************************************************
+   SUBROUTINE two_center_os_integrals(r12_operator, vab, dvab, rab, fba, fbb, omega, calculate_forces)
+
+      INTEGER, INTENT(IN)                                :: r12_operator
+      REAL(KIND=dp), DIMENSION(:, :), POINTER            :: vab
+      REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, &
+         POINTER                                         :: dvab
+      REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rab
+      TYPE(gto_basis_set_type), POINTER                  :: fba, fbb
+      REAL(KIND=dp), INTENT(IN), OPTIONAL                :: omega
+      LOGICAL, INTENT(IN)                                :: calculate_forces
+
+      CHARACTER(LEN=*), PARAMETER :: routineN = 'two_center_os_integrals', &
+         routineP = moduleN//':'//routineN
+
+      INTEGER                                            :: handle
+      REAL(KIND=dp)                                      :: my_omega
+
+      PROCEDURE(ab_sint_os), POINTER                     :: cps_operator2 => null()
+
+      CALL timeset(routineN, handle)
+
+      my_omega = 1.0_dp
+
+      SELECT CASE (r12_operator)
+      CASE (operator_coulomb)
+         cps_operator2 => cps_coulomb2
+      CASE (operator_verf)
+         cps_operator2 => cps_verf2
+         IF (PRESENT(omega)) my_omega = omega
+      CASE (operator_verfc)
+         cps_operator2 => cps_verfc2
+         IF (PRESENT(omega)) my_omega = omega
+      CASE (operator_vgauss)
+         cps_operator2 => cps_vgauss2
+         IF (PRESENT(omega)) my_omega = omega
+      CASE (operator_gauss)
+         cps_operator2 => cps_gauss2
+         IF (PRESENT(omega)) my_omega = omega
+      CASE DEFAULT
+         CPABORT("Operator not available")
+      END SELECT
+
+      CALL int_operator_ab_os(cps_operator2, vab, dvab, rab, fba, fbb, my_omega, calculate_forces)
+
+      CALL timestop(handle)
+
+   END SUBROUTINE two_center_os_integrals
+
+! **************************************************************************************************
+!> \brief calculate integrals (a|O(r12)|b)
+!> \param cps_operator2 procedure pointer for the respective operator.
+!> \param vab integral matrix of spherical contracted Gaussian functions
+!> \param dvab derivative of the integrals
+!> \param rab distance vector between center A and B
+!> \param fba basis at center A
+!> \param fbb basis at center B
+!> \param omega parameter in the operator
+!> \param calculate_forces ...
+! **************************************************************************************************
+   SUBROUTINE int_operator_ab_os(cps_operator2, vab, dvab, rab, fba, fbb, omega, calculate_forces)
+
+      PROCEDURE(ab_sint_os), POINTER                     :: cps_operator2
+      REAL(KIND=dp), DIMENSION(:, :), POINTER            :: vab
+      REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, &
+         POINTER                                         :: dvab
+      REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rab
+      TYPE(gto_basis_set_type), POINTER                  :: fba, fbb
+      REAL(KIND=dp), INTENT(IN)                          :: omega
+      LOGICAL, INTENT(IN)                                :: calculate_forces
+
+      CHARACTER(LEN=*), PARAMETER :: routineN = 'int_operator_ab_os', routineP = moduleN//':'//routineN
+
+      INTEGER                                            :: handle, i, iset, jset, lds, m1, m2, &
+                                                            maxco, maxcoa, maxcob, maxl, maxla, &
+                                                            maxlb, ncoa, ncoap, ncob, ncobp, &
+                                                            nseta, nsetb, sgfa, sgfb
+      INTEGER, DIMENSION(:), POINTER                     :: la_max, la_min, lb_max, lb_min, npgfa, &
+                                                            npgfb, nsgfa, nsgfb
+      INTEGER, DIMENSION(:, :), POINTER                  :: first_sgfa, first_sgfb
+      REAL(KIND=dp)                                      :: dab, rab2
+      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: vac, vac_plus
+      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: devab, vwork
+      REAL(KIND=dp), DIMENSION(:), POINTER               :: set_radius_a, set_radius_b
+      REAL(KIND=dp), DIMENSION(:, :), POINTER            :: rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb
+
+      CALL timeset(routineN, handle)
+      NULLIFY (la_max, la_min, lb_max, lb_min, npgfa, npgfb, nsgfa, nsgfb, &
+               first_sgfa, first_sgfb, set_radius_a, set_radius_b, rpgfa, rpgfb, &
+               sphi_a, sphi_b, zeta, zetb)
+
+      ! basis ikind
+      first_sgfa => fba%first_sgf
+      la_max => fba%lmax
+      la_min => fba%lmin
+      npgfa => fba%npgf
+      nseta = fba%nset
+      nsgfa => fba%nsgf_set
+      rpgfa => fba%pgf_radius
+      set_radius_a => fba%set_radius
+      sphi_a => fba%sphi
+      zeta => fba%zet
+      ! basis jkind
+      first_sgfb => fbb%first_sgf
+      lb_max => fbb%lmax
+      lb_min => fbb%lmin
+      npgfb => fbb%npgf
+      nsetb = fbb%nset
+      nsgfb => fbb%nsgf_set
+      rpgfb => fbb%pgf_radius
+      set_radius_b => fbb%set_radius
+      sphi_b => fbb%sphi
+      zetb => fbb%zet
+
+      CALL get_gto_basis_set(fba, maxco=maxcoa, maxl=maxla)
+      CALL get_gto_basis_set(fbb, maxco=maxcob, maxl=maxlb)
+      maxco = MAX(maxcoa, maxcob)
+      IF (calculate_forces) THEN
+         maxl = MAX(maxla+1, maxlb)
+      ELSE
+         maxl = MAX(maxla, maxlb)
+      ENDIF
+      lds = ncoset(maxl)
+
+      rab2 = SUM(rab*rab)
+      dab = SQRT(rab2)
+
+      DO iset = 1, nseta
+
+         ncoa = npgfa(iset)*ncoset(la_max(iset))
+         ncoap = npgfa(iset)*ncoset(la_max(iset)+1)
+         sgfa = first_sgfa(1, iset)
+
+         DO jset = 1, nsetb
+
+            ncob = npgfb(jset)*ncoset(lb_max(jset))
+            ncobp = npgfb(jset)*ncoset(lb_max(jset)+1)
+            sgfb = first_sgfb(1, jset)
+            m1 = sgfa+nsgfa(iset)-1
+            m2 = sgfb+nsgfb(jset)-1
+
+            ! calculate integrals
+            IF (calculate_forces) THEN
+               ALLOCATE (vwork(ncoap, ncobp, la_max(iset)+lb_max(jset)+3), &
+                         vac(ncoa, ncob), vac_plus(ncoap, ncobp), devab(ncoa, ncob, 3))
+               devab = 0._dp
+               vwork = 0.0_dp
+               vac = 0.0_dp
+               CALL operator2(cps_operator2, la_max(iset)+1, npgfa(iset), zeta(:, iset), la_min(iset), &
+                              lb_max(jset)+1, npgfb(jset), zetb(:, jset), lb_min(jset), &
+                              omega, rab, rab2, vac, vwork, maxder=1, vac_plus=vac_plus)
+               CALL dabdr_noscreen(la_max(iset), npgfa(iset), zeta(:, iset), lb_max(jset), npgfb(jset), &
+                                   vac_plus, devab(:, :, 1), devab(:, :, 2), devab(:, :, 3))
+               DO i = 1, 3
+                  CALL ab_contract(dvab(sgfa:m1, sgfb:m2, i), devab(:, :, i), sphi_a(:, sgfa:), &
+                                   sphi_b(:, sgfb:), ncoa, ncob, nsgfa(iset), nsgfb(jset))
+               ENDDO
+
+            ELSE
+               ALLOCATE (vwork(ncoa, ncob, la_max(iset)+lb_max(jset)+1), &
+                         vac(ncoa, ncob), vac_plus(ncoap, ncobp), devab(ncoa, ncob, 3))
+               vwork = 0.0_dp
+               vac = 0.0_dp
+               CALL operator2(cps_operator2, la_max(iset), npgfa(iset), zeta(:, iset), la_min(iset), &
+                              lb_max(jset), npgfb(jset), zetb(:, jset), lb_min(jset), &
+                              omega, rab, rab2, vac, vwork)
+            ENDIF
+
+            CALL ab_contract(vab(sgfa:m1, sgfb:m2), vac(1:ncoa, 1:ncob), sphi_a(:, sgfa:), sphi_b(:, sgfb:), &
+                             ncoa, ncob, nsgfa(iset), nsgfb(jset))
+            DEALLOCATE (vwork, vac, vac_plus, devab)
+         END DO
+      END DO
+
+      CALL timestop(handle)
+
+   END SUBROUTINE int_operator_ab_os
+
 ! **************************************************************************************************
 !> \brief calculate overlap integrals (a,b)
 !> \param sab integral (a,b)
 !> \param dsab derivative of sab with respect to A
-!> \param ra ...
-!> \param rb ...
-!> \param rab ...
+!> \param ra position center A, for debug only
+!> \param rb position center B, for debug only
+!> \param rab distance vector between center A and B
 !> \param fba basis at center A
 !> \param fbb basis at center B
 !> \param calculate_forces ...
 !> \param debug integrals are debugged by recursive routines if requested
 !> \param dmax maximal deviation between integrals when debugging
 ! **************************************************************************************************
-   SUBROUTINE lri_int_ab(sab, dsab, ra, rb, rab, fba, fbb, calculate_forces, debug, &
-                         dmax)
+   SUBROUTINE int_overlap_ab_os(sab, dsab, ra, rb, rab, fba, fbb, calculate_forces, debug, &
+                                dmax)
 
       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: sab
       REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, &
@@ -62,7 +269,8 @@ CONTAINS
       LOGICAL, INTENT(IN)                                :: calculate_forces, debug
       REAL(KIND=dp), INTENT(INOUT)                       :: dmax
 
-      CHARACTER(LEN=*), PARAMETER :: routineN = 'lri_int_ab', routineP = moduleN//':'//routineN
+      CHARACTER(LEN=*), PARAMETER :: routineN = 'int_overlap_ab_os', &
+         routineP = moduleN//':'//routineN
 
       INTEGER                                            :: handle, i, iset, jset, lds, m1, m2, &
                                                             maxco, maxcoa, maxcob, maxl, maxla, &
@@ -181,15 +389,15 @@ CONTAINS
 
       CALL timestop(handle)
 
-   END SUBROUTINE lri_int_ab
+   END SUBROUTINE int_overlap_ab_os
 
 ! **************************************************************************************************
 !> \brief calculate integrals (a,b,fa)
 !> \param abaint integral (a,b,fa)
 !> \param dabdaint derivative of abaint with respect to A
-!> \param ra ...
-!> \param rb ...
-!> \param rab ...
+!> \param ra position center A, for debug only
+!> \param rb position center B, for debug only
+!> \param rab distance vector between center A and B
 !> \param oba orbital basis at center A
 !> \param obb orbital basis at center B
 !> \param fba auxiliary basis set at center A
@@ -197,8 +405,8 @@ CONTAINS
 !> \param debug integrals are debugged by recursive routines if requested
 !> \param dmax maximal deviation between integrals when debugging
 ! **************************************************************************************************
-   SUBROUTINE lri_int_aba(abaint, dabdaint, ra, rb, rab, oba, obb, fba, &
-                          calculate_forces, debug, dmax)
+   SUBROUTINE int_overlap_aba_os(abaint, dabdaint, ra, rb, rab, oba, obb, fba, &
+                                 calculate_forces, debug, dmax)
 
       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: abaint
       REAL(KIND=dp), DIMENSION(:, :, :, :), OPTIONAL, &
@@ -208,7 +416,8 @@ CONTAINS
       LOGICAL, INTENT(IN)                                :: calculate_forces, debug
       REAL(KIND=dp), INTENT(INOUT)                       :: dmax
 
-      CHARACTER(LEN=*), PARAMETER :: routineN = 'lri_int_aba', routineP = moduleN//':'//routineN
+      CHARACTER(LEN=*), PARAMETER :: routineN = 'int_overlap_aba_os', &
+         routineP = moduleN//':'//routineN
 
       INTEGER                                            :: handle, i, iset, jset, kaset, m1, m2, &
                                                             m3, ncoa, ncob, ncoc, nseta, nsetb, &
@@ -345,15 +554,15 @@ CONTAINS
 
       CALL timestop(handle)
 
-   END SUBROUTINE lri_int_aba
+   END SUBROUTINE int_overlap_aba_os
 
 ! **************************************************************************************************
 !> \brief calculate integrals (a,b,fb)
 !> \param abbint integral (a,b,fb)
 !> \param dabbint derivative of abbint with respect to A
-!> \param ra ...
-!> \param rb ...
-!> \param rab ...
+!> \param ra position center A, for debug only
+!> \param rb position center B, for debug only
+!> \param rab distance vector between center A and B
 !> \param oba orbital basis at center A
 !> \param obb orbital basis at center B
 !> \param fbb auxiliary basis set at center B
@@ -361,8 +570,8 @@ CONTAINS
 !> \param debug integrals are debugged by recursive routines if requested
 !> \param dmax maximal deviation between integrals when debugging
 ! **************************************************************************************************
-   SUBROUTINE lri_int_abb(abbint, dabbint, ra, rb, rab, oba, obb, fbb, calculate_forces, &
-                          debug, dmax)
+   SUBROUTINE int_overlap_abb_os(abbint, dabbint, ra, rb, rab, oba, obb, fbb, calculate_forces, &
+                                 debug, dmax)
 
       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: abbint
       REAL(KIND=dp), DIMENSION(:, :, :, :), OPTIONAL, &
@@ -372,7 +581,8 @@ CONTAINS
       LOGICAL, INTENT(IN)                                :: calculate_forces, debug
       REAL(KIND=dp), INTENT(INOUT)                       :: dmax
 
-      CHARACTER(LEN=*), PARAMETER :: routineN = 'lri_int_abb', routineP = moduleN//':'//routineN
+      CHARACTER(LEN=*), PARAMETER :: routineN = 'int_overlap_abb_os', &
+         routineP = moduleN//':'//routineN
 
       INTEGER                                            :: handle, i, iset, jset, kbset, m1, m2, &
                                                             m3, ncoa, ncob, ncoc, nseta, nsetb, &
@@ -501,7 +711,7 @@ CONTAINS
 
       CALL timestop(handle)
 
-   END SUBROUTINE lri_int_abb
+   END SUBROUTINE int_overlap_abb_os
 
 ! **************************************************************************************************
 !> \brief calculate overlap integrals (aa,bb)
@@ -514,7 +724,7 @@ CONTAINS
 !> \param debug integrals are debugged by recursive routines if requested
 !> \param dmax maximal deviation between integrals when debugging
 ! **************************************************************************************************
-   SUBROUTINE lri_int_aabb(saabb, oba, obb, rab, ra, rb, debug, dmax)
+   SUBROUTINE int_overlap_aabb_os(saabb, oba, obb, rab, ra, rb, debug, dmax)
 
       REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER      :: saabb
       TYPE(gto_basis_set_type), POINTER                  :: oba, obb
@@ -522,7 +732,8 @@ CONTAINS
       LOGICAL, INTENT(IN)                                :: debug
       REAL(KIND=dp), INTENT(INOUT)                       :: dmax
 
-      CHARACTER(LEN=*), PARAMETER :: routineN = 'lri_int_aabb', routineP = moduleN//':'//routineN
+      CHARACTER(LEN=*), PARAMETER :: routineN = 'int_overlap_aabb_os', &
+         routineP = moduleN//':'//routineN
 
       INTEGER :: handle, iset, isgfa1, jset, jsgfa2, kset, ksgfb1, lds, lset, lsgfb2, m1, m2, m3, &
          m4, maxco, maxcoa, maxcob, maxl, maxla, maxlb, ncoa1, ncoa2, ncob1, ncob2, nseta, nsetb, &
@@ -649,193 +860,6 @@ CONTAINS
 
       CALL timestop(handle)
 
-   END SUBROUTINE lri_int_aabb
+   END SUBROUTINE int_overlap_aabb_os
 
-! **************************************************************************************************
-!> \brief contract overlap integrals (a,b) and transfer to spherical Gaussians
-!> \param abint ...
-!> \param sab ...
-!> \param sphi_a ...
-!> \param sphi_b ...
-!> \param ncoa ...
-!> \param ncob ...
-!> \param nsgfa ...
-!> \param nsgfb ...
-! **************************************************************************************************
-   SUBROUTINE ab_contract(abint, sab, sphi_a, sphi_b, ncoa, ncob, nsgfa, nsgfb)
-
-      REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT)      :: abint
-      REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: sab, sphi_a, sphi_b
-      INTEGER, INTENT(IN)                                :: ncoa, ncob, nsgfa, nsgfb
-
-      CHARACTER(LEN=*), PARAMETER :: routineN = 'ab_contract', routineP = moduleN//':'//routineN
-
-      INTEGER                                            :: m1, m2, msphia, msphib, nn
-      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: cpp
-
-      msphia = SIZE(sphi_a, 1)
-      msphib = SIZE(sphi_b, 1)
-
-      m1 = SIZE(sab, 1)
-      m2 = SIZE(sab, 2)
-
-      nn = SIZE(abint, 1)
-
-      ALLOCATE (cpp(nsgfa, m2))
-
-      CALL dgemm("T", "N", nsgfa, m2, ncoa, 1._dp, sphi_a, msphia, sab, m1, 0.0_dp, cpp, nsgfa)
-      CALL dgemm("N", "N", nsgfa, nsgfb, ncob, 1._dp, cpp, nsgfa, sphi_b, msphib, 0.0_dp, &
-                 abint, nn)
-
-      DEALLOCATE (cpp)
-
-   END SUBROUTINE ab_contract
-
-! **************************************************************************************************
-!> \brief contract three-center overlap integrals (a,b,c) and transfer
-!>        to spherical Gaussians
-!> \param abcint ...
-!> \param sabc ...
-!> \param sphi_a ...
-!> \param sphi_b ...
-!> \param sphi_c ...
-!> \param ncoa ...
-!> \param ncob ...
-!> \param ncoc ...
-!> \param nsgfa ...
-!> \param nsgfb ...
-!> \param nsgfc ...
-! **************************************************************************************************
-   SUBROUTINE abc_contract(abcint, sabc, sphi_a, sphi_b, sphi_c, ncoa, ncob, ncoc, &
-                           nsgfa, nsgfb, nsgfc)
-
-      REAL(KIND=dp), DIMENSION(:, :, :)                  :: abcint, sabc
-      REAL(KIND=dp), DIMENSION(:, :)                     :: sphi_a, sphi_b, sphi_c
-      INTEGER, INTENT(IN)                                :: ncoa, ncob, ncoc, nsgfa, nsgfb, nsgfc
-
-      CHARACTER(LEN=*), PARAMETER :: routineN = 'abc_contract', routineP = moduleN//':'//routineN
-
-      INTEGER                                            :: handle, isgfc, m1, m2, m3, msphia, &
-                                                            msphib, msphic
-      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: temp_ccc, work_cpc
-      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: cpc, cpp
-
-      CALL timeset(routineN, handle)
-
-      msphia = SIZE(sphi_a, 1)
-      msphib = SIZE(sphi_b, 1)
-      msphic = SIZE(sphi_c, 1)
-
-      m1 = SIZE(sabc, 1)
-      m2 = SIZE(sabc, 2)
-      m3 = SIZE(sabc, 3)
-
-      ALLOCATE (cpp(nsgfa, m2, m3), cpc(nsgfa, m2, nsgfc))
-      cpp = 0._dp
-      cpc = 0._dp
-      ALLOCATE (work_cpc(nsgfa, m2), temp_ccc(nsgfa, nsgfb))
-      work_cpc(:, :) = 0._dp
-      temp_ccc(:, :) = 0._dp
-
-      CALL dgemm("T", "N", nsgfa, m2*m3, ncoa, 1._dp, sphi_a, msphia, sabc, m1, 0.0_dp, cpp, nsgfa)
-      CALL dgemm("N", "N", nsgfa*m2, nsgfc, ncoc, 1._dp, cpp, nsgfa*m2, sphi_c, msphic, 0.0_dp, &
-                 cpc, nsgfa*m2)
-
-      DO isgfc = 1, nsgfc
-         work_cpc(:, :) = cpc(:, :, isgfc)
-         CALL dgemm("N", "N", nsgfa, nsgfb, ncob, 1._dp, work_cpc, nsgfa, sphi_b, msphib, &
-                    0.0_dp, temp_ccc, nsgfa)
-         abcint(:, :, isgfc) = temp_ccc(:, :)
-      END DO
-
-      DEALLOCATE (cpp, cpc, work_cpc, temp_ccc)
-
-      CALL timestop(handle)
-
-   END SUBROUTINE abc_contract
-
-! **************************************************************************************************
-!> \brief contract four-center overlap integrals (a,b,c,d) and transfer
-!>        to spherical Gaussians
-!> \param abcdint ...
-!> \param sabcd ...
-!> \param sphi_a ...
-!> \param sphi_b ...
-!> \param sphi_c ...
-!> \param sphi_d ...
-!> \param ncoa ...
-!> \param ncob ...
-!> \param ncoc ...
-!> \param ncod ...
-!> \param nsgfa ...
-!> \param nsgfb ...
-!> \param nsgfc ...
-!> \param nsgfd ...
-! **************************************************************************************************
-   SUBROUTINE abcd_contract(abcdint, sabcd, sphi_a, sphi_b, sphi_c, sphi_d, ncoa, ncob, &
-                            ncoc, ncod, nsgfa, nsgfb, nsgfc, nsgfd)
-
-      REAL(KIND=dp), DIMENSION(:, :, :, :), &
-         INTENT(INOUT)                                   :: abcdint
-      REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN)   :: sabcd
-      REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: sphi_a, sphi_b, sphi_c, sphi_d
-      INTEGER, INTENT(IN)                                :: ncoa, ncob, ncoc, ncod, nsgfa, nsgfb, &
-                                                            nsgfc, nsgfd
-
-      CHARACTER(LEN=*), PARAMETER :: routineN = 'abcd_contract', routineP = moduleN//':'//routineN
-
-      INTEGER                                            :: handle, isgfc, isgfd, m1, m2, m3, m4, &
-                                                            msphia, msphib, msphic, msphid
-      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: temp_cccc, work_cpcc
-      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: temp_cpcc, work_cppc
-      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :)  :: cpcc, cppc, cppp
-
-      CALL timeset(routineN, handle)
-
-      msphia = SIZE(sphi_a, 1)
-      msphib = SIZE(sphi_b, 1)
-      msphic = SIZE(sphi_c, 1)
-      msphid = SIZE(sphi_d, 1)
-
-      m1 = SIZE(sabcd, 1)
-      m2 = SIZE(sabcd, 2)
-      m3 = SIZE(sabcd, 3)
-      m4 = SIZE(sabcd, 4)
-
-      ALLOCATE (cppp(nsgfa, m2, m3, m4), cppc(nsgfa, m2, m3, nsgfd), &
-                cpcc(nsgfa, m2, nsgfc, nsgfd))
-
-      ALLOCATE (work_cppc(nsgfa, m2, m3), temp_cpcc(nsgfa, m2, nsgfc))
-      work_cppc = 0._dp
-      temp_cpcc = 0._dp
-
-      ALLOCATE (work_cpcc(nsgfa, m2), temp_cccc(nsgfa, nsgfb))
-      work_cpcc = 0._dp
-      temp_cccc = 0._dp
-
-      CALL dgemm("T", "N", nsgfa, m2*m3*m4, ncoa, 1._dp, sphi_a, msphia, sabcd, m1, &
-                 0.0_dp, cppp, nsgfa)
-      CALL dgemm("N", "N", nsgfa*m2*m3, nsgfd, ncod, 1._dp, cppp, nsgfa*m2*m3, &
-                 sphi_d, msphid, 0.0_dp, cppc, nsgfa*m2*m3)
-
-      DO isgfd = 1, nsgfd
-         work_cppc(:, :, :) = cppc(:, :, :, isgfd)
-         CALL dgemm("N", "N", nsgfa*m2, nsgfc, ncoc, 1._dp, work_cppc, nsgfa*m2, &
-                    sphi_c, msphic, 0.0_dp, temp_cpcc, nsgfa*m2)
-         cpcc(:, :, :, isgfd) = temp_cpcc(:, :, :)
-         DO isgfc = 1, nsgfc
-            work_cpcc(:, :) = cpcc(:, :, isgfc, isgfd)
-            CALL dgemm("N", "N", nsgfa, nsgfb, ncob, 1._dp, work_cpcc, nsgfa, sphi_b, &
-                       msphib, 0.0_dp, temp_cccc, nsgfa)
-            abcdint(:, :, isgfc, isgfd) = temp_cccc(:, :)
-         END DO
-      END DO
-
-      DEALLOCATE (cpcc, cppc, cppp)
-      DEALLOCATE (work_cpcc, work_cppc, temp_cpcc, temp_cccc)
-
-      CALL timestop(handle)
-
-   END SUBROUTINE abcd_contract
-
-END MODULE lri_os_integrals
+END MODULE generic_os_integrals

src/generic_shg_integrals.F

@@ -0,0 +1,849 @@
+!--------------------------------------------------------------------------------------------------!
+!   CP2K: A general program to perform molecular dynamics simulations                              !
+!   Copyright (C) 2000 - 2016  CP2K developers group                                               !
+!--------------------------------------------------------------------------------------------------!
+
+! **************************************************************************************************
+!> \brief Calculation of contracted, spherical Gaussian integrals using the solid harmonic
+!>        Gaussian (SHG) integral scheme. Routines for the following two-center integrals:
+!>        i)  (a|O(r12)|b) where O(r12) is the overlap, coulomb operator etc.
+!>        ii) (aba) and (abb) s-overlaps
+!> \par Literature
+!>      T.J. Giese and D. M. York, J. Chem. Phys, 128, 064104 (2008)
+!>      T. Helgaker, P Joergensen, J. Olsen, Molecular Electronic-Structure
+!>                                           Theory, Wiley
+!> \par History
+!>      created [05.2016]
+!> \author Dorothea Golze
+! **************************************************************************************************
+MODULE generic_shg_integrals
+   USE ai_construct_shg,                ONLY: &
+        construct_dev_shg_ab, construct_int_shg_ab, construct_overlap_shg_aba, &
+        construct_overlap_shg_abb, dev_overlap_shg_aba, dev_overlap_shg_abb, get_W_matrix, &
+        get_dW_matrix, get_real_scaled_solid_harmonic
+   USE ai_s_contract_shg,               ONLY: &
+        contract_s_overlap_aba, contract_s_overlap_abb, contract_sint_ab_chigh, &
+        contract_sint_ab_clow, s_coulomb_ab, s_gauss_ab, s_overlap_ab, s_overlap_abx, s_verf_ab, &
+        s_verfc_ab, s_vgauss_ab
+   USE basis_set_types,                 ONLY: gto_basis_set_type
+   USE input_constants,                 ONLY: operator_coulomb,&
+                                              operator_gauss,&
+                                              operator_verf,&
+                                              operator_verfc,&
+                                              operator_vgauss
+   USE kinds,                           ONLY: dp
+   USE orbital_pointers,                ONLY: nsoset
+#include "base/base_uses.f90"
+
+   IMPLICIT NONE
+
+   PRIVATE
+
+   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'generic_shg_integrals'
+
+   PUBLIC :: two_center_shg_integrals, lri_precalc_angular_shg_part, &
+             int_overlap_ab_shg, int_overlap_aba_shg, int_overlap_abb_shg
+
+   ABSTRACT INTERFACE
+! **************************************************************************************************
+!> \brief Interface for the calculation of integrals over s-functions and their scalar derivatives
+!>        with respect to rab2
+!> \param la_max ...
+!> \param npgfa ...
+!> \param zeta ...
+!> \param lb_max ...
+!> \param npgfb ...
+!> \param zetb ...
+!> \param omega ...
+!> \param rab ...
+!> \param v matrix storing the integrals and scalar derivatives
+!> \param calculate_forces ...
+! **************************************************************************************************
+      SUBROUTINE ab_sint_shg(la_max, npgfa, zeta, lb_max, npgfb, zetb, omega, rab, v, calculate_forces)
+      USE kinds,                        ONLY: dp
+      INTEGER, INTENT(IN)                                :: la_max, npgfa
+      REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: zeta
+      INTEGER, INTENT(IN)                                :: lb_max, npgfb
+      REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: zetb
+      REAL(KIND=dp), INTENT(IN)                          :: omega
+      REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rab
+      REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT)   :: v
+      LOGICAL, INTENT(IN)                                :: calculate_forces
+
+      END SUBROUTINE ab_sint_shg
+   END INTERFACE
+
+! **************************************************************************************************
+
+CONTAINS
+
+! **************************************************************************************************
+!> \brief Calcululates the two-center integrals of the type (a|O(r12)|b) using the SHG scheme
+!> \param r12_operator the integral operator, which depends on r12=|r1-r2|
+!> \param vab integral matrix of spherical contracted Gaussian functions
+!> \param dvab derivative of the integrals
+!> \param rab distance vector between center A and B
+!> \param fba basis at center A
+!> \param fbb basis at center B
+!> \param scona_shg SHG contraction matrix for A
+!> \param sconb_shg SHG contraction matrix for B
+!> \param omega parameter in the operator
+!> \param calculate_forces ...
+! **************************************************************************************************
+   SUBROUTINE two_center_shg_integrals(r12_operator, vab, dvab, rab, fba, fbb, scona_shg, sconb_shg, &
+                                       omega, calculate_forces)
+
+      INTEGER, INTENT(IN)                                :: r12_operator
+      REAL(KIND=dp), DIMENSION(:, :), POINTER            :: vab
+      REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: dvab
+      REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rab
+      TYPE(gto_basis_set_type), POINTER                  :: fba, fbb
+      REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: scona_shg, sconb_shg
+      REAL(KIND=dp), INTENT(IN), OPTIONAL                :: omega
+      LOGICAL, INTENT(IN)                                :: calculate_forces
+
+      CHARACTER(LEN=*), PARAMETER :: routineN = 'two_center_shg_integrals', &
+         routineP = moduleN//':'//routineN
+
+      INTEGER                                            :: handle
+      REAL(KIND=dp)                                      :: my_omega
+      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: Waux_mat
+      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :)  :: dWaux_mat
+
+      PROCEDURE(ab_sint_shg), POINTER                    :: s_operator_ab => null()
+
+      CALL timeset(routineN, handle)
+
+      CALL precalc_angular_shg_part(fba, fbb, rab, Waux_mat, dWaux_mat, calculate_forces)
+      my_omega = 1.0_dp
+
+      SELECT CASE (r12_operator)
+      CASE (operator_coulomb)
+         s_operator_ab => s_coulomb_ab
+      CASE (operator_verf)
+         s_operator_ab => s_verf_ab
+         IF (PRESENT(omega)) my_omega = omega
+      CASE (operator_verfc)
+         s_operator_ab => s_verfc_ab
+         IF (PRESENT(omega)) my_omega = omega
+      CASE (operator_vgauss)
+         s_operator_ab => s_vgauss_ab
+         IF (PRESENT(omega)) my_omega = omega
+      CASE (operator_gauss)
+         s_operator_ab => s_gauss_ab
+         IF (PRESENT(omega)) my_omega = omega
+      CASE DEFAULT
+         CPABORT("Operator not available")
+      END SELECT
+
+      CALL int_operator_ab_shg(s_operator_ab, vab, dvab, rab, fba, fbb, scona_shg, sconb_shg, &
+                               my_omega, Waux_mat, dWaux_mat, calculate_forces)
+
+      DEALLOCATE (Waux_mat, dWaux_mat)
+
+      CALL timestop(handle)
+
+   END SUBROUTINE two_center_shg_integrals
+
+! **************************************************************************************************
+!> \brief precalculates the angular part of the SHG integrals
+!> \param fba basis at center A
+!> \param fbb basis at center B
+!> \param rab distance vector between a and b
+!> \param Waux_mat W matrix that contains the angular-dependent part
+!> \param dWaux_mat derivative of the W matrix
+!> \param calculate_forces ...
+! **************************************************************************************************
+   SUBROUTINE precalc_angular_shg_part(fba, fbb, rab, Waux_mat, dWaux_mat, calculate_forces)
+
+      TYPE(gto_basis_set_type), POINTER                  :: fba, fbb
+      REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rab
+      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
+         INTENT(OUT)                                     :: Waux_mat
+      REAL(KIND=dp), ALLOCATABLE, &
+         DIMENSION(:, :, :, :), INTENT(OUT)              :: dWaux_mat
+      LOGICAL, INTENT(IN)                                :: calculate_forces
+
+      CHARACTER(len=*), PARAMETER :: routineN = 'precalc_angular_shg_part', &
+         routineP = moduleN//':'//routineN
+
+      INTEGER                                            :: handle, la_max, lb_max, lmax, mdim(3)
+      INTEGER, DIMENSION(:), POINTER                     :: la_max_all
+      REAL(KIND=dp)                                      :: rab2
+      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: Rc, Rs
+
+      CALL timeset(routineN, handle)
+
+      NULLIFY (la_max_all)
+      rab2 = rab(1)*rab(1)+rab(2)*rab(2)+rab(3)*rab(3)
+
+      la_max = MAXVAL(fba%lmax)
+      lb_max = MAXVAL(fbb%lmax)
+      lmax = MAX(la_max, lb_max)
+
+      ALLOCATE (la_max_all(0:lb_max))
+      ALLOCATE (Rc(0:lmax, -2*lmax:2*lmax), Rs(0:lmax, -2*lmax:2*lmax))
+      Rc = 0._dp
+      Rs = 0._dp
+      mdim(1) = MIN(la_max, lb_max)+1
+      mdim(2) = nsoset(la_max)+1
+      mdim(3) = nsoset(lb_max)+1
+      ALLOCATE (Waux_mat(mdim(1), mdim(2), mdim(3)))
+      ALLOCATE (dWaux_mat(3, mdim(1), mdim(2), mdim(3)))
+
+      la_max_all(0:lb_max) = la_max
+      !*** -rab, since Eq. in Ref. use Ra-Rb, not Rb-Ra
+      CALL get_real_scaled_solid_harmonic(Rc, Rs, lmax, -rab, rab2)
+      CALL get_W_matrix(la_max_all, lb_max, lmax, Rc, Rs, Waux_mat)
+      IF (calculate_forces) THEN
+         CALL get_dW_matrix(la_max_all, lb_max, Waux_mat, dWaux_mat)
+      ENDIF
+
+      DEALLOCATE (Rc, Rs, la_max_all)
+      CALL timestop(handle)
+
+   END SUBROUTINE precalc_angular_shg_part
+
+! **************************************************************************************************
+!> \brief calculate integrals (a|O(r12)|b)
+!> \param s_operator_ab procedure pointer for the respective operator. The integral evaluation
+!>        differs only in the calculation of the [s|O(r12)|s] integrals and their scalar
+!>        derivatives.
+!> \param vab integral matrix of spherical contracted Gaussian functions
+!> \param dvab derivative of the integrals
+!> \param rab distance vector between center A and B
+!> \param fba basis at center A
+!> \param fbb basis at center B
+!> \param scona_shg SHG contraction matrix for A
+!> \param sconb_shg SHG contraction matrix for B
+!> \param omega parameter in the operator
+!> \param Waux_mat W matrix that contains the angular-dependent part
+!> \param dWaux_mat derivative of the W matrix
+!> \param calculate_forces ...
+! **************************************************************************************************
+   SUBROUTINE int_operator_ab_shg(s_operator_ab, vab, dvab, rab, fba, fbb, scona_shg, sconb_shg, &
+                                  omega, Waux_mat, dWaux_mat, calculate_forces)
+
+      PROCEDURE(ab_sint_shg), POINTER                    :: s_operator_ab
+      REAL(KIND=dp), DIMENSION(:, :), POINTER            :: vab
+      REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: dvab
+      REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rab
+      TYPE(gto_basis_set_type), POINTER                  :: fba, fbb
+      REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: scona_shg, sconb_shg
+      REAL(KIND=dp), INTENT(IN)                          :: omega
+      REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN)      :: Waux_mat
+      REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN)   :: dWaux_mat
+      LOGICAL, INTENT(IN)                                :: calculate_forces
+
+      CHARACTER(LEN=*), PARAMETER :: routineN = 'int_operator_ab_shg', routineP = moduleN//':'//routineN
+
+      INTEGER :: handle, iset, jset, la_max_set, lb_max_set, ndev, nds, nds_max, npgfa_set, &
+                 npgfb_set, nseta, nsetb, nsgfa_set, nsgfb_set, nshella_set, nshellb_set, sgfa, sgfb
+      INTEGER, DIMENSION(:), POINTER                     :: la_max, lb_max, npgfa, npgfb, nsgfa, &
+                                                            nsgfb, nshella, nshellb
+      INTEGER, DIMENSION(:, :), POINTER                  :: first_sgfa, first_sgfb, la, lb
+      REAL(KIND=dp)                                      :: dab
+      REAL(KIND=dp), DIMENSION(:), POINTER               :: set_radius_a, set_radius_b
+      REAL(KIND=dp), DIMENSION(:, :), POINTER            :: zeta, zetb
+      REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: swork, swork_cont
+
+      CALL timeset(routineN, handle)
+      NULLIFY (la_max, lb_max, npgfa, npgfb, first_sgfa, first_sgfb, set_radius_a, &
+               set_radius_b, zeta, zetb)
+
+      ! basis ikind
+      first_sgfa => fba%first_sgf
+      la_max => fba%lmax
+      la => fba%l
+      npgfa => fba%npgf
+      nsgfa => fba%nsgf_set
+      nseta = fba%nset
+      set_radius_a => fba%set_radius
+      zeta => fba%zet
+      nshella => fba%nshell
+      ! basis jkind
+      first_sgfb => fbb%first_sgf
+      lb_max => fbb%lmax
+      lb => fbb%l
+      npgfb => fbb%npgf
+      nsgfb => fbb%nsgf_set
+      nsetb = fbb%nset
+      set_radius_b => fbb%set_radius
+      zetb => fbb%zet
+      nshellb => fbb%nshell
+
+      dab = SQRT(SUM(rab**2))
+
+      la_max_set = MAXVAL(la_max)
+      lb_max_set = MAXVAL(lb_max)
+
+      ! allocate some work matrices
+      npgfa_set = MAXVAL(npgfa)
+      npgfb_set = MAXVAL(npgfb)
+      nshella_set = MAXVAL(nshella)
+      nshellb_set = MAXVAL(nshellb)
+      nsgfa_set = MAXVAL(nsgfa)
+      nsgfb_set = MAXVAL(nsgfb)
+      ndev = 0
+      IF (calculate_forces) ndev = 1
+      nds_max = la_max_set+lb_max_set+ndev+1
+      ALLOCATE (swork(npgfa_set, npgfb_set, nds_max))
+      ALLOCATE (swork_cont(nds_max, nshella_set, nshellb_set))
+
+      vab = 0.0_dp
+      dvab = 0.0_dp
+
+      DO iset = 1, nseta
+         sgfa = first_sgfa(1, iset)
+
+         DO jset = 1, nsetb
+
+            sgfb = first_sgfb(1, jset)
+            CALL s_operator_ab(la_max(iset), npgfa(iset), zeta(:, iset), &
+                               lb_max(jset), npgfb(jset), zetb(:, jset), &
+                               omega, rab, swork, calculate_forces)
+            nds = la_max(iset)+lb_max(jset)+ndev+1
+            CALL contract_sint_ab_chigh(npgfa(iset), nshella(iset), scona_shg(1:npgfa(iset), 1:nshella(iset), iset), &
+                                        npgfb(jset), nshellb(jset), sconb_shg(1:npgfb(jset), 1:nshellb(jset), jset), &
+                                        nds, swork(1:npgfa(iset), 1:npgfb(jset), 1:nds), &
+                                        swork_cont(1:nds, 1:nshella(iset), 1:nshellb(jset)))
+            CALL construct_int_shg_ab(la(:, iset), first_sgfa(:, iset), nshella(iset), &
+                                      lb(:, jset), first_sgfb(:, jset), nshellb(jset), &
+                                      swork_cont, Waux_mat, vab)
+            IF (calculate_forces) THEN
+               !*** -rab, since Eq. in Ref. use Ra-Rb, not Rb-Ra
+               CALL construct_dev_shg_ab(la(:, iset), first_sgfa(:, iset), nshella(iset), &
+                                         lb(:, jset), first_sgfb(:, jset), nshellb(jset), &
+                                         -rab, swork_cont, Waux_mat, dWaux_mat, dvab)
+            ENDIF
+         END DO
+      END DO
+
+      DEALLOCATE (swork, swork_cont)
+
+      CALL timestop(handle)
+   END SUBROUTINE int_operator_ab_shg
+
+! **************************************************************************************************
+!> \brief calculate overlap integrals (a,b)
+!> \param sab integral (a,b)
+!> \param dsab derivative of sab
+!> \param rab distance vector
+!> \param fba basis at center A
+!> \param fbb basis at center B
+!> \param scona_shg contraction matrix A
+!> \param sconb_shg contraxtion matrix B
+!> \param Waux_mat W matrix that contains the angular-dependent part
+!> \param dWaux_mat derivative of the W matrix
+!> \param calculate_forces ...
+! **************************************************************************************************
+   SUBROUTINE int_overlap_ab_shg(sab, dsab, rab, fba, fbb, scona_shg, sconb_shg, Waux_mat, dWaux_mat, &
+                                 calculate_forces)
+
+      REAL(KIND=dp), DIMENSION(:, :), POINTER            :: sab
+      REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: dsab
+      REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rab
+      TYPE(gto_basis_set_type), POINTER                  :: fba, fbb
+      REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: scona_shg, sconb_shg
+      REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN)      :: Waux_mat
+      REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN)   :: dWaux_mat
+      LOGICAL, INTENT(IN)                                :: calculate_forces
+
+      CHARACTER(LEN=*), PARAMETER :: routineN = 'int_overlap_ab_shg', &
+         routineP = moduleN//':'//routineN
+
+      INTEGER                                            :: handle, iset, jset, la_max_set, &
+                                                            lb_max_set, ndev, nds_max, npgfa_set, &
+                                                            npgfb_set, nseta, nsetb, nshella_set, &
+                                                            nshellb_set
+      INTEGER, DIMENSION(:), POINTER                     :: la_max, lb_max, npgfa, npgfb, nshella, &
+                                                            nshellb
+      INTEGER, DIMENSION(:, :), POINTER                  :: first_sgfa, first_sgfb, la, lb
+      REAL(KIND=dp)                                      :: dab
+      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: swork, swork_cont
+      REAL(KIND=dp), DIMENSION(:), POINTER               :: set_radius_a, set_radius_b
+      REAL(KIND=dp), DIMENSION(:, :), POINTER            :: zeta, zetb
+
+      CALL timeset(routineN, handle)
+      NULLIFY (la_max, lb_max, npgfa, npgfb, first_sgfa, first_sgfb, set_radius_a, &
+               set_radius_b, zeta, zetb)
+
+      ! basis ikind
+      first_sgfa => fba%first_sgf
+      la_max => fba%lmax
+      la => fba%l
+      npgfa => fba%npgf
+      nseta = fba%nset
+      set_radius_a => fba%set_radius
+      zeta => fba%zet
+      nshella => fba%nshell
+      ! basis jkind
+      first_sgfb => fbb%first_sgf
+      lb_max => fbb%lmax
+      lb => fbb%l
+      npgfb => fbb%npgf
+      nsetb = fbb%nset
+      set_radius_b => fbb%set_radius
+      zetb => fbb%zet
+      nshellb => fbb%nshell
+
+      dab = SQRT(SUM(rab**2))
+
+      la_max_set = MAXVAL(la_max)
+      lb_max_set = MAXVAL(lb_max)
+
+      ! allocate some work matrices
+      npgfa_set = MAXVAL(npgfa)
+      npgfb_set = MAXVAL(npgfb)
+      nshella_set = MAXVAL(nshella)
+      nshellb_set = MAXVAL(nshellb)
+      ndev = 0
+      IF (calculate_forces) ndev = 1
+      nds_max = la_max_set+lb_max_set+ndev+1
+      ALLOCATE (swork(npgfa_set, npgfb_set, nds_max))
+      ALLOCATE (swork_cont(nds_max, nshella_set, nshellb_set))
+
+      DO iset = 1, nseta
+
+         DO jset = 1, nsetb
+
+            IF (set_radius_a(iset)+set_radius_b(jset) < dab) CYCLE
+
+            CALL s_overlap_ab(la_max(iset), npgfa(iset), zeta(:, iset), &
+                              lb_max(jset), npgfb(jset), zetb(:, jset), &
+                              rab, swork, calculate_forces)
+            CALL contract_sint_ab_clow(la(:, iset), npgfa(iset), nshella(iset), &
+                                       zeta(:, iset), scona_shg(:, :, iset), &
+                                       lb(:, jset), npgfb(jset), nshellb(jset), &
+                                       zetb(:, jset), sconb_shg(:, :, jset), &
+                                       swork, swork_cont, calculate_forces)
+            CALL construct_int_shg_ab(la(:, iset), first_sgfa(:, iset), nshella(iset), &
+                                      lb(:, jset), first_sgfb(:, jset), nshellb(jset), &
+                                      swork_cont, Waux_mat, sab)
+            IF (calculate_forces) THEN
+               !*** -rab, since Eq. in Ref. use Ra-Rb, not Rb-Ra
+               CALL construct_dev_shg_ab(la(:, iset), first_sgfa(:, iset), nshella(iset), &
+                                         lb(:, jset), first_sgfb(:, jset), nshellb(jset), &
+                                         -rab, swork_cont, Waux_mat, dWaux_mat, dsab)
+            ENDIF
+         END DO
+      END DO
+
+      DEALLOCATE (swork, swork_cont)
+
+      CALL timestop(handle)
+
+   END SUBROUTINE int_overlap_ab_shg
+
+! **************************************************************************************************
+!> \brief calculate integrals (a,b,fa)
+!> \param abaint integral (a,b,fa)
+!> \param dabdaint ...
+!> \param rab distance vector between A and B
+!> \param oba orbital basis at center A
+!> \param obb orbital basis at center B
+!> \param fba auxiliary basis set at center A
+!> \param scona_shg contraction matrix for orb bas on A
+!> \param sconb_shg contraction matrix for orb bas on B
+!> \param sconca_shg contraction matrix for aux bas on A
+!> \param cg_coeff Clebsch-Gordon coefficients
+!> \param cg_none0_list list of none-zero Clebsch-Gordon coefficients
+!> \param ncg_none0 number of non-zero Clebsch-Gordon coefficients
+!> \param Waux_mat W matrix that contains the angular-dependent part
+!> \param dWaux_mat derivative of the W matrix
+!> \param calculate_forces ...
+! **************************************************************************************************
+   SUBROUTINE int_overlap_aba_shg(abaint, dabdaint, rab, oba, obb, fba, scona_shg, &
+                                  sconb_shg, sconca_shg, cg_coeff, cg_none0_list, ncg_none0, &
+                                  Waux_mat, dWaux_mat, calculate_forces)
+
+      REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: abaint
+      REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER      :: dabdaint
+      REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rab
+      TYPE(gto_basis_set_type), POINTER                  :: oba, obb, fba
+      REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN)      :: scona_shg, sconb_shg, sconca_shg, &
+                                                            cg_coeff
+      INTEGER, DIMENSION(:, :, :), INTENT(IN)            :: cg_none0_list
+      INTEGER, DIMENSION(:, :), INTENT(IN)               :: ncg_none0
+      REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN)      :: Waux_mat
+      REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN)   :: dWaux_mat
+      LOGICAL, INTENT(IN)                                :: calculate_forces
+
+      CHARACTER(LEN=*), PARAMETER :: routineN = 'int_overlap_aba_shg', &
+         routineP = moduleN//':'//routineN
+
+      INTEGER :: handle, iset, jset, kset, la_max_set, laa_max_set, lb_max_set, lca_max_set, ndev, &
+         nds_max, nl_set, npgfa_set, npgfb_set, npgfca_set, nseta, nsetb, nsetca, nshella_set, &
+         nshellb_set, nshellca_set
+      INTEGER, DIMENSION(:), POINTER                     :: la_max, lb_max, lca_max, npgfa, npgfb, &
+                                                            npgfca, nshella, nshellb, nshellca
+      INTEGER, DIMENSION(:, :), POINTER                  :: first_sgfa, first_sgfb, first_sgfca, la, &
+                                                            lb, lca
+      REAL(KIND=dp)                                      :: dab, rab2
+      REAL(KIND=dp), ALLOCATABLE, &
+         DIMENSION(:, :, :, :, :)                        :: swork, swork_cont
+      REAL(KIND=dp), DIMENSION(:), POINTER               :: set_radius_a, set_radius_b, set_radius_ca
+      REAL(KIND=dp), DIMENSION(:, :), POINTER            :: zeta, zetb, zetca
+
+      CALL timeset(routineN, handle)
+      NULLIFY (la_max, lb_max, lca_max, npgfa, npgfb, npgfca)
+      NULLIFY (first_sgfa, first_sgfb, first_sgfca, set_radius_a, set_radius_b, &
+               set_radius_ca, zeta, zetb, zetca)
+
+      ! basis ikind
+      first_sgfa => oba%first_sgf
+      la_max => oba%lmax
+      la => oba%l
+      npgfa => oba%npgf
+      nshella => oba%nshell
+      nseta = oba%nset
+      set_radius_a => oba%set_radius
+      zeta => oba%zet
+      ! basis jkind
+      first_sgfb => obb%first_sgf
+      lb_max => obb%lmax
+      lb => obb%l
+      npgfb => obb%npgf
+      nshellb => obb%nshell
+      nsetb = obb%nset
+      set_radius_b => obb%set_radius
+      zetb => obb%zet
+
+      ! basis RI A
+      first_sgfca => fba%first_sgf
+      lca_max => fba%lmax
+      lca => fba%l
+      npgfca => fba%npgf
+      nshellca => fba%nshell
+      nsetca = fba%nset
+      set_radius_ca => fba%set_radius
+      zetca => fba%zet
+
+      dab = SQRT(SUM(rab**2))
+      rab2 = rab(1)*rab(1)+rab(2)*rab(2)+rab(3)*rab(3)
+
+      la_max_set = MAXVAL(la_max)
+      lb_max_set = MAXVAL(lb_max)
+      lca_max_set = MAXVAL(lca_max)
+      npgfa_set = MAXVAL(npgfa)
+      npgfb_set = MAXVAL(npgfb)
+      npgfca_set = MAXVAL(npgfca)
+      nshella_set = MAXVAL(nshella)
+      nshellb_set = MAXVAL(nshellb)
+      nshellca_set = MAXVAL(nshellca)
+      !*** for forces: derivative+1 in auxiliary vector required
+      ndev = 0
+      IF (calculate_forces) ndev = 1
+
+      laa_max_set = la_max_set+lca_max_set
+
+      ! allocate some work storage....
+      nds_max = laa_max_set+lb_max_set+ndev+1
+      nl_set = INT((laa_max_set)/2)
+      ALLOCATE (swork(npgfa_set, npgfb_set, npgfca_set, 0:nl_set, nds_max))
+      ALLOCATE (swork_cont(nds_max, 0:nl_set, nshella_set, nshellb_set, nshellca_set))
+
+      DO iset = 1, nseta
+
+         DO jset = 1, nsetb
+
+            IF (set_radius_a(iset)+set_radius_b(jset) < dab) CYCLE
+
+            DO kset = 1, nsetca
+
+               IF (set_radius_b(jset)+set_radius_ca(kset) < dab) CYCLE
+
+               CALL s_overlap_abx(la_max(iset), npgfa(iset), zeta(:, iset), &
+                                  lb_max(jset), npgfb(jset), zetb(:, jset), &
+                                  lca_max(kset), npgfca(kset), zetca(:, kset), &
+                                  rab, swork, calculate_forces, calc_aba=.TRUE.)
+
+               CALL contract_s_overlap_aba(la(:, iset), npgfa(iset), nshella(iset), zeta(:, iset), scona_shg(:, :, iset), &
+                                           lb(:, jset), npgfb(jset), nshellb(jset), zetb(:, jset), sconb_shg(:, :, jset), &
+                                           lca(:, kset), npgfca(kset), nshellca(kset), zetca(:, kset), sconca_shg(:, :, kset), &
+                                           swork, swork_cont, calculate_forces)
+
+               CALL construct_overlap_shg_aba(la(:, iset), first_sgfa(:, iset), nshella(iset), &
+                                              lb(:, jset), first_sgfb(:, jset), nshellb(jset), &
+                                              lca(:, kset), first_sgfca(:, kset), nshellca(kset), &
+                                              cg_coeff, cg_none0_list, ncg_none0, &
+                                              swork_cont, Waux_mat, abaint)
+               IF (calculate_forces) THEN
+                  !*** -rab, since Eq. in Ref. use Ra-Rb, not Rb-Ra
+                  CALL dev_overlap_shg_aba(la(:, iset), first_sgfa(:, iset), nshella(iset), &
+                                           lb(:, jset), first_sgfb(:, jset), nshellb(jset), &
+                                           lca(:, kset), first_sgfca(:, kset), nshellca(kset), &
+                                           cg_coeff, cg_none0_list, ncg_none0, &
+                                           -rab, swork_cont, Waux_mat, dWaux_mat, dabdaint)
+               ENDIF
+
+            END DO
+         END DO
+      END DO
+
+      DEALLOCATE (swork_cont)
+      DEALLOCATE (swork)
+
+      CALL timestop(handle)
+
+   END SUBROUTINE int_overlap_aba_shg
+
+! **************************************************************************************************
+!> \brief calculate integrals (a,b,fb)
+!> \param abbint integral (a,b,fb)
+!> \param dabbint derivative of abbint
+!> \param rab distance vector between A and B
+!> \param oba orbital basis at center A
+!> \param obb orbital basis at center B
+!> \param fbb auxiliary basis set at center B
+!> \param scona_shg contraction matrix for orb bas on A
+!> \param sconb_shg contraction matrix for orb bas on B
+!> \param sconcb_shg contraction matrix for aux bas on B
+!> \param cg_coeff Clebsch-Gordon coefficients
+!> \param cg_none0_list list of none-zero Clebsch-Gordon coefficients
+!> \param ncg_none0 number of non-zero Clebsch-Gordon coefficients
+!> \param Waux_mat W matrix that contains the angular-dependent part
+!> \param dWaux_mat derivative of the W matrix
+!> \param calculate_forces ...
+! **************************************************************************************************
+   SUBROUTINE int_overlap_abb_shg(abbint, dabbint, rab, oba, obb, fbb, scona_shg, &
+                                  sconb_shg, sconcb_shg, cg_coeff, cg_none0_list, ncg_none0, &
+                                  Waux_mat, dWaux_mat, calculate_forces)
+
+      REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: abbint
+      REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER      :: dabbint
+      REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rab
+      TYPE(gto_basis_set_type), POINTER                  :: oba, obb, fbb
+      REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN)      :: scona_shg, sconb_shg, sconcb_shg, &
+                                                            cg_coeff
+      INTEGER, DIMENSION(:, :, :), INTENT(IN)            :: cg_none0_list
+      INTEGER, DIMENSION(:, :), INTENT(IN)               :: ncg_none0
+      REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN)      :: Waux_mat
+      REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN)   :: dWaux_mat
+      LOGICAL, INTENT(IN)                                :: calculate_forces
+
+      CHARACTER(LEN=*), PARAMETER :: routineN = 'int_overlap_abb_shg', &
+         routineP = moduleN//':'//routineN
+
+      INTEGER :: handle, iset, jset, kset, la_max_set, lb_max_set, lbb_max_set, lcb_max_set, ndev, &
+         nds_max, nl_set, npgfa_set, npgfb_set, npgfcb_set, nseta, nsetb, nsetcb, nshella_set, &
+         nshellb_set, nshellcb_set
+      INTEGER, DIMENSION(:), POINTER                     :: la_max, lb_max, lcb_max, npgfa, npgfb, &
+                                                            npgfcb, nshella, nshellb, nshellcb
+      INTEGER, DIMENSION(:, :), POINTER                  :: first_sgfa, first_sgfb, first_sgfcb, la, &
+                                                            lb, lcb
+      REAL(KIND=dp)                                      :: dab, rab2
+      REAL(KIND=dp), ALLOCATABLE, &
+         DIMENSION(:, :, :, :, :)                        :: swork, swork_cont
+      REAL(KIND=dp), DIMENSION(:), POINTER               :: set_radius_a, set_radius_b, set_radius_cb
+      REAL(KIND=dp), DIMENSION(:, :), POINTER            :: zeta, zetb, zetcb
+
+      CALL timeset(routineN, handle)
+      NULLIFY (la_max, lb_max, lcb_max, npgfa, npgfb, npgfcb)
+      NULLIFY (first_sgfa, first_sgfb, first_sgfcb, set_radius_a, set_radius_b, &
+               set_radius_cb, zeta, zetb, zetcb)
+
+      ! basis ikind
+      first_sgfa => oba%first_sgf
+      la_max => oba%lmax
+      la => oba%l
+      npgfa => oba%npgf
+      nshella => oba%nshell
+      nseta = oba%nset
+      set_radius_a => oba%set_radius
+      zeta => oba%zet
+      ! basis jkind
+      first_sgfb => obb%first_sgf
+      lb_max => obb%lmax
+      lb => obb%l
+      npgfb => obb%npgf
+      nshellb => obb%nshell
+      nsetb = obb%nset
+      set_radius_b => obb%set_radius
+      zetb => obb%zet
+
+      ! basis RI on B
+      first_sgfcb => fbb%first_sgf
+      lcb_max => fbb%lmax
+      lcb => fbb%l
+      npgfcb => fbb%npgf
+      nshellcb => fbb%nshell
+      nsetcb = fbb%nset
+      set_radius_cb => fbb%set_radius
+      zetcb => fbb%zet
+
+      dab = SQRT(SUM(rab**2))
+      rab2 = rab(1)*rab(1)+rab(2)*rab(2)+rab(3)*rab(3)
+
+      la_max_set = MAXVAL(la_max)
+      lb_max_set = MAXVAL(lb_max)
+      lcb_max_set = MAXVAL(lcb_max)
+      npgfa_set = MAXVAL(npgfa)
+      npgfb_set = MAXVAL(npgfb)
+      npgfcb_set = MAXVAL(npgfcb)
+      nshella_set = MAXVAL(nshella)
+      nshellb_set = MAXVAL(nshellb)
+      nshellcb_set = MAXVAL(nshellcb)
+      !*** for forces: derivative+1 in auxiliary vector required
+      ndev = 0
+      IF (calculate_forces) ndev = 1
+
+      lbb_max_set = lb_max_set+lcb_max_set
+
+      ! allocate some work storage....
+      nds_max = la_max_set+lbb_max_set+ndev+1
+      nl_set = INT((lbb_max_set)/2)
+      ALLOCATE (swork(1:npgfa_set, 1:npgfb_set, 1:npgfcb_set, 0:nl_set, 1:nds_max))
+      ALLOCATE (swork_cont(nds_max, 0:nl_set, nshella_set, nshellb_set, nshellcb_set))
+
+      DO iset = 1, nseta
+
+         DO jset = 1, nsetb
+
+            IF (set_radius_a(iset)+set_radius_b(jset) < dab) CYCLE
+
+            DO kset = 1, nsetcb
+
+               IF (set_radius_a(iset)+set_radius_cb(kset) < dab) CYCLE
+
+               CALL s_overlap_abx(la_max(iset), npgfa(iset), zeta(:, iset), &
+                                  lb_max(jset), npgfb(jset), zetb(:, jset), &
+                                  lcb_max(kset), npgfcb(kset), zetcb(:, kset), &
+                                  rab, swork, calculate_forces, calc_aba=.FALSE.)
+
+               CALL contract_s_overlap_abb(la(:, iset), npgfa(iset), nshella(iset), zeta(:, iset), scona_shg(:, :, iset), &
+                                           lb(:, jset), npgfb(jset), nshellb(jset), zetb(:, jset), sconb_shg(:, :, jset), &
+                                           lcb(:, kset), npgfcb(kset), nshellcb(kset), zetcb(:, kset), sconcb_shg(:, :, kset), &
+                                           swork, swork_cont, calculate_forces)
+
+               CALL construct_overlap_shg_abb(la(:, iset), first_sgfa(:, iset), nshella(iset), &
+                                              lb(:, jset), first_sgfb(:, jset), nshellb(jset), &
+                                              lcb(:, kset), first_sgfcb(:, kset), nshellcb(kset), &
+                                              cg_coeff, cg_none0_list, &
+                                              ncg_none0, swork_cont, Waux_mat, abbint)
+               IF (calculate_forces) THEN
+                  !*** -rab, since Eq. in Ref. use Ra-Rb, not Rb-Ra
+                  CALL dev_overlap_shg_abb(la(:, iset), first_sgfa(:, iset), nshella(iset), &
+                                           lb(:, jset), first_sgfb(:, jset), nshellb(jset), &
+                                           lcb(:, kset), first_sgfcb(:, kset), nshellcb(kset), &
+                                           cg_coeff, cg_none0_list, ncg_none0, -rab, swork_cont, &
+                                           Waux_mat, dWaux_mat, dabbint)
+               ENDIF
+
+            END DO
+         END DO
+      END DO
+
+      DEALLOCATE (swork_cont)
+      DEALLOCATE (swork)
+
+      CALL timestop(handle)
+
+   END SUBROUTINE int_overlap_abb_shg
+
+! **************************************************************************************************
+!> \brief precalculates the angular part of the SHG integrals for the matrices
+!>        (fa,fb), (a,b), (a,b,fa) and (b,fb,a); the same Waux_mat can then be used for all
+!>        for integrals; specific for LRIGPW
+!> \param oba orbital basis on a
+!> \param obb orbital basis on b
+!> \param fba aux basis on a
+!> \param fbb aux basis on b
+!> \param rab distance vector between a and b
+!> \param Waux_mat W matrix that contains the angular-dependent part
+!> \param dWaux_mat derivative of the W matrix
+!> \param calculate_forces ...
+! **************************************************************************************************
+   SUBROUTINE lri_precalc_angular_shg_part(oba, obb, fba, fbb, rab, Waux_mat, dWaux_mat, calculate_forces)
+
+      TYPE(gto_basis_set_type), POINTER                  :: oba, obb, fba, fbb
+      REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rab
+      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
+         INTENT(OUT)                                     :: Waux_mat
+      REAL(KIND=dp), ALLOCATABLE, &
+         DIMENSION(:, :, :, :), INTENT(OUT)              :: dWaux_mat
+      LOGICAL, INTENT(IN)                                :: calculate_forces
+
+      CHARACTER(len=*), PARAMETER :: routineN = 'lri_precalc_angular_shg_part', &
+         routineP = moduleN//':'//routineN
+
+      INTEGER                                            :: handle, i, isize, j, k, la_max, laa_max, &
+                                                            lb_max, lbb_max, lca_max, lcb_max, &
+                                                            li_max, lj_max, lmax, mdim(3), &
+                                                            size_int(4, 2), temp
+      INTEGER, DIMENSION(:), POINTER                     :: li_max_all
+      REAL(KIND=dp)                                      :: rab2
+      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: Rc, Rs
+
+      CALL timeset(routineN, handle)
+
+      rab2 = rab(1)*rab(1)+rab(2)*rab(2)+rab(3)*rab(3)
+
+      !*** 1 Waux_mat of size (li_max,lj_max) for elements
+      !                    i        j
+      !    [aab]    --> (laa_max, lb_max)
+      !    [bba]    --> (lbb_max, la_max) --> use for [abb]
+      !    [ab] ri  --> (lca_max, lcb_max)
+      !    [ab] orb --> (la_max , lb_max)
+
+      la_max = MAXVAL(oba%lmax)
+      lb_max = MAXVAL(obb%lmax)
+      lca_max = MAXVAL(fba%lmax)
+      lcb_max = MAXVAL(fbb%lmax)
+
+      laa_max = la_max+lca_max
+      lbb_max = lb_max+lcb_max
+      li_max = MAX(laa_max, lbb_max)
+      lj_max = MAX(la_max, lb_max, lcb_max)
+      lmax = li_max
+
+      ALLOCATE (li_max_all(0:lj_max))
+      ALLOCATE (Rc(0:lmax, -2*lmax:2*lmax), Rs(0:lmax, -2*lmax:2*lmax))
+      Rc = 0._dp
+      Rs = 0._dp
+      mdim(1) = li_max+lj_max+1
+      mdim(2) = nsoset(li_max)+1
+      mdim(3) = nsoset(lj_max)+1
+      ALLOCATE (Waux_mat(mdim(1), mdim(2), mdim(3)))
+      ALLOCATE (dWaux_mat(3, mdim(1), mdim(2), mdim(3)))
+      !Waux_mat = 0._dp !.. takes time
+      !dWaux_mat =0._dp !.. takes time
+
+      !*** Waux_mat (li_max,lj_max) contains elements not needed,
+      !*** make indixing so that only required ones are computed
+      !*** li_max_all(j) --> li_max dependent on j
+      size_int(1, :) = (/laa_max, lb_max/)
+      size_int(2, :) = (/lbb_max, la_max/)
+      size_int(3, :) = (/lca_max, lcb_max/)
+      size_int(4, :) = (/la_max, lb_max/)
+
+      li_max_all(:) = 0
+      DO isize = 1, 4
+         i = size_int(isize, 1)
+         j = size_int(isize, 2)
+         k = li_max_all(j)
+         IF (k < i) li_max_all(j) = i
+      ENDDO
+      temp = li_max_all(lj_max)
+      DO j = lj_max, 0, -1
+         IF (li_max_all(j) < temp) THEN
+            li_max_all(j) = temp
+         ELSE
+            temp = li_max_all(j)
+         ENDIF
+      ENDDO
+
+      !*** -rab, since Eq. in Ref. use Ra-Rb, not Rb-Ra
+      CALL get_real_scaled_solid_harmonic(Rc, Rs, lmax, -rab, rab2)
+      CALL get_W_matrix(li_max_all, lj_max, lmax, Rc, Rs, Waux_mat)
+      IF (calculate_forces) THEN
+         CALL get_dW_matrix(li_max_all, lj_max, Waux_mat, dWaux_mat)
+      ENDIF
+
+      DEALLOCATE (Rc, Rs, li_max_all)
+
+      CALL timestop(handle)
+
+   END SUBROUTINE lri_precalc_angular_shg_part
+
+END MODULE generic_shg_integrals

src/input_constants.F

@@ -989,4 +989,10 @@ MODULE input_constants
    INTEGER, PARAMETER, PUBLIC               :: do_eri_gpw = 0, &
                                                do_eri_mme = 1
 
+   ! operators for integrals
+   INTEGER, PARAMETER, PUBLIC                :: operator_coulomb = 0, &
+                                                operator_verf = 1, &
+                                                operator_verfc = 2, &
+                                                operator_vgauss = 3, &
+                                                operator_gauss = 4
 END MODULE input_constants

src/input_cp2k_subsys.F

@@ -62,7 +62,8 @@ MODULE input_cp2k_subsys
    PUBLIC :: create_subsys_section, &
              create_cell_section, &
              create_structure_data_section, &
-             create_rng_section
+             create_rng_section, &
+             create_basis_section
 !***
 CONTAINS
 

src/library_tests.F

@@ -105,6 +105,7 @@ MODULE library_tests
         pw2rs, realspace_grid_desc_type, realspace_grid_input_type, realspace_grid_type, rs2pw, &
         rs_grid_create, rs_grid_create_descriptor, rs_grid_print, rs_grid_release, &
         rs_grid_release_descriptor, rs_grid_zero, rs_pw_transfer
+   USE shg_integrals_test,              ONLY: shg_integrals_perf_acc_test
    USE termination,                     ONLY: stop_memory
 #include "./base/base_uses.f90"
 
@@ -141,7 +142,8 @@ CONTAINS
       LOGICAL                                            :: explicit
       TYPE(cp_logger_type), POINTER                      :: logger
       TYPE(section_vals_type), POINTER :: cp_dbcsr_test_section, cp_fm_gemm_test_section, &
-         eigensolver_section, eri_mme_test_section, pw_transfer_section, rs_pw_transfer_section
+         eigensolver_section, eri_mme_test_section, pw_transfer_section, rs_pw_transfer_section, &
+         shg_integrals_test_section
 
       CALL timeset(routineN, handle)
 
@@ -206,6 +208,12 @@ CONTAINS
          CALL eri_mme_perf_acc_test(para_env, iw, eri_mme_test_section)
       ENDIF
 
+      shg_integrals_test_section => section_vals_get_subs_vals(root_section, "TEST%SHG_INTEGRALS_TEST")
+      CALL section_vals_get(shg_integrals_test_section, explicit=explicit)
+      IF (explicit) THEN
+         CALL shg_integrals_perf_acc_test(iw, shg_integrals_test_section)
+      ENDIF
+
       ! DBCSR tests
       ! matrix_multiplication
       cp_dbcsr_test_section => section_vals_get_subs_vals(root_section, &

src/lri_environment_init.F

@@ -410,7 +410,7 @@ CONTAINS
 
    END SUBROUTINE basis_ovlp
 
-!*****************************************************************************
+! **************************************************************************************************
 !> \brief contraction matrix for SHG integrals
 !> \param basis ...
 !> \param scon_shg contraction matrix
@@ -438,6 +438,7 @@ CONTAINS
       maxshell = SIZE(basis%gcc, 2)
       ALLOCATE (norm(basis%nset, maxshell))
       ALLOCATE (scon_shg(maxpgf, maxshell, nset))
+      scon_shg = 0.0_dp
 
       CALL basis_norm_shg(basis, norm)
 

src/lri_environment_methods.F

@@ -29,6 +29,14 @@ MODULE lri_environment_methods
    USE cp_output_handling,              ONLY: cp_print_key_finished_output,&
                                               cp_print_key_unit_nr
    USE cp_para_types,                   ONLY: cp_para_env_type
+   USE generic_os_integrals,            ONLY: int_overlap_aabb_os,&
+                                              int_overlap_ab_os,&
+                                              int_overlap_aba_os,&
+                                              int_overlap_abb_os
+   USE generic_shg_integrals,           ONLY: int_overlap_ab_shg,&
+                                              int_overlap_aba_shg,&
+                                              int_overlap_abb_shg,&
+                                              lri_precalc_angular_shg_part
    USE input_constants,                 ONLY: do_lri_inv,&
                                               do_lri_inv_auto,&
                                               do_lri_pseudoinv_diag,&
@@ -40,14 +48,6 @@ MODULE lri_environment_methods
         deallocate_lri_ints, deallocate_lri_ints_rho, lri_density_create, lri_density_release, &
         lri_density_type, lri_environment_type, lri_int_rho_type, lri_int_type, lri_kind_type, &
         lri_list_type, lri_rhoab_type
-   USE lri_os_integrals,                ONLY: lri_int_aabb,&
-                                              lri_int_ab,&
-                                              lri_int_aba,&
-                                              lri_int_abb
-   USE lri_shg_integrals,               ONLY: lri_int_ab_shg,&
-                                              lri_int_aba_shg,&
-                                              lri_int_abb_shg,&
-                                              precalc_angular_shg_part
    USE mathlib,                         ONLY: get_pseudo_inverse_diag,&
                                               get_pseudo_inverse_svd,&
                                               invmat_symm
@@ -165,7 +165,6 @@ CONTAINS
             CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
                                    nlist=nlist, ilist=ilist, nnode=nneighbor, inode=jneighbor, &
                                    iatom=iatom, jatom=jatom, r=rab)
-
             iac = ikind+nkind*(jkind-1)
             dab = SQRT(SUM(rab*rab))
 
@@ -200,7 +199,6 @@ CONTAINS
                CALL lri_int_os(lri_env, lrii, ra, rb, rab, obasa, obasb, fbasa, fbasb, &
                                iatom, jatom, ikind)
             ENDIF
-
             ! construct and invert S matrix
             lrii%sinv(1:nfa, 1:nfa) = lri_env%bas_prop(ikind)%ri_ovlp(1:nfa, 1:nfa)
             IF (dab > lri_env%delta) THEN
@@ -311,8 +309,8 @@ CONTAINS
             rb(:) = pbc(particle_set(jatom)%r, cell)
 
             ! calculate integrals (aa,bb)
-            CALL lri_int_aabb(lriir%soaabb, obasa, obasb, rab, ra, rb, lri_env%debug, &
-                              lriir%dmax_aabb)
+            CALL int_overlap_aabb_os(lriir%soaabb, obasa, obasb, rab, ra, rb, lri_env%debug, &
+                                     lriir%dmax_aabb)
 
          END DO
 
@@ -678,26 +676,26 @@ CONTAINS
          lrii%soo(1:nba, 1:nba) = lri_env%bas_prop(ikind)%orb_ovlp(1:nba, 1:nba)
          lrii%dsoo = 0._dp
          !*** (a,a,fa)
-         CALL lri_int_aba(lrii%abaint, ra=ra, rb=rb, rab=rab, oba=obasa, obb=obasb, &
-                          fba=fbasa, calculate_forces=.FALSE., debug=lri_env%debug, &
-                          dmax=lrii%dmax_aba)
+         CALL int_overlap_aba_os(lrii%abaint, ra=ra, rb=rb, rab=rab, oba=obasa, obb=obasb, &
+                                 fba=fbasa, calculate_forces=.FALSE., debug=lri_env%debug, &
+                                 dmax=lrii%dmax_aba)
          lrii%dabdaint = 0.0_dp
          lrii%dabbint = 0.0_dp
       ELSE
          !*** (fa,fb)
-         CALL lri_int_ab(lrii%sab, lrii%dsab, ra, rb, rab, fbasa, fbasb, &
-                         lrii%calc_force_pair, lri_env%debug, &
-                         lrii%dmax_ab)
+         CALL int_overlap_ab_os(lrii%sab, lrii%dsab, ra, rb, rab, fbasa, fbasb, &
+                                lrii%calc_force_pair, lri_env%debug, &
+                                lrii%dmax_ab)
          !*** (a,b)
-         CALL lri_int_ab(lrii%soo, lrii%dsoo, ra, rb, rab, obasa, obasb, &
-                         lrii%calc_force_pair, lri_env%debug, &
-                         lrii%dmax_oo)
+         CALL int_overlap_ab_os(lrii%soo, lrii%dsoo, ra, rb, rab, obasa, obasb, &
+                                lrii%calc_force_pair, lri_env%debug, &
+                                lrii%dmax_oo)
          !*** (a,b,fa)
-         CALL lri_int_aba(lrii%abaint, lrii%dabdaint, ra, rb, rab, obasa, obasb, fbasa, &
-                          lrii%calc_force_pair, lri_env%debug, lrii%dmax_aba)
+         CALL int_overlap_aba_os(lrii%abaint, lrii%dabdaint, ra, rb, rab, obasa, obasb, fbasa, &
+                                 lrii%calc_force_pair, lri_env%debug, lrii%dmax_aba)
          !*** (a,b,fb)
-         CALL lri_int_abb(lrii%abbint, lrii%dabbint, ra, rb, rab, obasa, obasb, fbasb, &
-                          lrii%calc_force_pair, lri_env%debug, lrii%dmax_abb)
+         CALL int_overlap_abb_os(lrii%abbint, lrii%dabbint, ra, rb, rab, obasa, obasb, fbasb, &
+                                 lrii%calc_force_pair, lri_env%debug, lrii%dmax_abb)
       ENDIF
 
       CALL timestop(handle)
@@ -733,13 +731,12 @@ CONTAINS
 
       INTEGER                                            :: handle, nba, nbb, nfa, nfb
       REAL(KIND=dp)                                      :: dab
+      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: Waux_mat
+      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :)  :: dWaux_mat
       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: scon_orba, scon_orbb, scon_ria, scon_rib
-      REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER      :: Waux_mat
-      REAL(KIND=dp), DIMENSION(:, :, :, :, :), POINTER   :: dWaux_mat
 
       CALL timeset(routineN, handle)
-      NULLIFY (scon_orba, scon_orbb, scon_ria, scon_rib, Waux_mat, dWaux_mat)
-
+      NULLIFY (scon_orba, scon_orbb, scon_ria, scon_rib)
       dab = SQRT(SUM(rab*rab))
       nba = obasa%nsgf
       nbb = obasb%nsgf
@@ -755,9 +752,9 @@ CONTAINS
          lrii%soo(1:nba, 1:nba) = lri_env%bas_prop(ikind)%orb_ovlp(1:nba, 1:nba)
          lrii%dsoo = 0._dp
          !*** (a,a,fa)
-         CALL lri_int_aba(lrii%abaint, ra=ra, rb=rb, rab=rab, oba=obasa, obb=obasb, &
-                          fba=fbasa, calculate_forces=.FALSE., debug=lri_env%debug, &
-                          dmax=lrii%dmax_aba)
+         CALL int_overlap_aba_os(lrii%abaint, ra=ra, rb=rb, rab=rab, oba=obasa, obb=obasb, &
+                                 fba=fbasa, calculate_forces=.FALSE., debug=lri_env%debug, &
+                                 dmax=lrii%dmax_aba)
          lrii%dabdaint = 0.0_dp
          lrii%dabbint = 0.0_dp
       ELSE
@@ -765,22 +762,23 @@ CONTAINS
          scon_orbb => lri_env%bas_prop(jkind)%scon_orb
          scon_ria => lri_env%bas_prop(ikind)%scon_ri
          scon_rib => lri_env%bas_prop(jkind)%scon_ri
-         CALL precalc_angular_shg_part(obasa, obasb, fbasa, fbasb, rab, Waux_mat, dWaux_mat, &
-                                       lrii%calc_force_pair)
+         CALL lri_precalc_angular_shg_part(obasa, obasb, fbasa, fbasb, rab, Waux_mat, dWaux_mat, lrii%calc_force_pair)
          !*** (fa,fb)
-         CALL lri_int_ab_shg(lrii%sab, lrii%dsab, rab, fbasa, fbasb, scon_ria, scon_rib, &
-                             Waux_mat, dWaux_mat, lrii%calc_force_pair)
+         CALL int_overlap_ab_shg(lrii%sab, lrii%dsab, rab, fbasa, fbasb, scon_ria, scon_rib, &
+                                 Waux_mat, dWaux_mat, lrii%calc_force_pair)
          !*** (a,b)
-         CALL lri_int_ab_shg(lrii%soo, lrii%dsoo, rab, obasa, obasb, scon_orba, scon_orbb, &
-                             Waux_mat, dWaux_mat, lrii%calc_force_pair)
+         CALL int_overlap_ab_shg(lrii%soo, lrii%dsoo, rab, obasa, obasb, scon_orba, scon_orbb, &
+                                 Waux_mat, dWaux_mat, lrii%calc_force_pair)
          !*** (a,b,fa)
-         CALL lri_int_aba_shg(lri_env, lrii%abaint, lrii%dabdaint, rab, obasa, obasb, fbasa, &
-                              scon_orba, scon_orbb, scon_ria, Waux_mat, dWaux_mat, &
-                              lrii%calc_force_pair)
+         CALL int_overlap_aba_shg(lrii%abaint, lrii%dabdaint, rab, obasa, obasb, fbasa, &
+                                  scon_orba, scon_orbb, scon_ria, &
+                                  lri_env%cg_shg%cg_coeff, lri_env%cg_shg%cg_none0_list, &
+                                  lri_env%cg_shg%ncg_none0, Waux_mat, dWaux_mat, lrii%calc_force_pair)
          !*** (a,b,fb)
-         CALL lri_int_abb_shg(lri_env, lrii%abbint, lrii%dabbint, rab, obasa, obasb, fbasb, &
-                              scon_orba, scon_orbb, scon_rib, Waux_mat, dWaux_mat, &
-                              lrii%calc_force_pair)
+         CALL int_overlap_abb_shg(lrii%abbint, lrii%dabbint, rab, obasa, obasb, fbasb, &
+                                  scon_orba, scon_orbb, scon_rib, &
+                                  lri_env%cg_shg%cg_coeff, lri_env%cg_shg%cg_none0_list, &
+                                  lri_env%cg_shg%ncg_none0, Waux_mat, dWaux_mat, lrii%calc_force_pair)
          DEALLOCATE (Waux_mat, dWaux_mat)
       ENDIF
 

src/lri_environment_types.F

@@ -72,25 +72,25 @@ MODULE lri_environment_types
       ! condition number of overlap matrix
       REAL(KIND=dp)                                           :: cond_num
       ! integrals (a,b,ai)
-      REAL(KIND=dp), DIMENSION(:, :, :), POINTER                :: abaint
-      ! integrals (a,b,bi)
-      REAL(KIND=dp), DIMENSION(:, :, :), POINTER                :: abbint
+      REAL(KIND=dp), DIMENSION(:, :, :), POINTER              :: abaint
+      ! integrals (a,b,b)
+      REAL(KIND=dp), DIMENSION(:, :, :), POINTER              :: abbint
       ! integrals (da/dA,b,dai/dA)
-      REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER              :: dabdaint
+      REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER           :: dabdaint
       ! integrals (da/dA,b,bi)
-      REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER              :: dabbint
+      REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER           :: dabbint
       ! integrals (a,b)
-      REAL(KIND=dp), DIMENSION(:, :), POINTER                  :: soo
+      REAL(KIND=dp), DIMENSION(:, :), POINTER                 :: soo
       ! derivative d(a,b)/dA
-      REAL(KIND=dp), DIMENSION(:, :, :), POINTER                :: dsoo
+      REAL(KIND=dp), DIMENSION(:, :, :), POINTER              :: dsoo
       ! integrals (ai,bi)
-      REAL(KIND=dp), DIMENSION(:, :), POINTER                  :: sab
+      REAL(KIND=dp), DIMENSION(:, :), POINTER                 :: sab
       ! derivative d(ai,bi)/dA
-      REAL(KIND=dp), DIMENSION(:, :, :), POINTER                :: dsab
+      REAL(KIND=dp), DIMENSION(:, :, :), POINTER              :: dsab
       ! derivative of fit coeff dacoef/dpmatrix
-      REAL(KIND=dp), DIMENSION(:, :, :), POINTER                :: dacoef
+      REAL(KIND=dp), DIMENSION(:, :, :), POINTER              :: dacoef
       ! inverse of integrals (ai,bi)
-      REAL(KIND=dp), DIMENSION(:, :), POINTER                  :: sinv
+      REAL(KIND=dp), DIMENSION(:, :), POINTER                 :: sinv
       ! integral (ai) / (bi), dim(1..nfa,nfa+1..nfa+nfb)
       REAL(KIND=dp), DIMENSION(:), POINTER                    :: n
       ! sinv * (ai)
@@ -142,13 +142,13 @@ MODULE lri_environment_types
       ! integral of ri basis fbas
       REAL(KIND=dp), DIMENSION(:), POINTER                    :: int_fbas
       ! self overlap ri basis
-      REAL(KIND=dp), DIMENSION(:, :), POINTER                  :: ri_ovlp
+      REAL(KIND=dp), DIMENSION(:, :), POINTER                 :: ri_ovlp
       ! self overlap orb basis
-      REAL(KIND=dp), DIMENSION(:, :), POINTER                  :: orb_ovlp
+      REAL(KIND=dp), DIMENSION(:, :), POINTER                 :: orb_ovlp
       ! contraction matrix for SHG integrals ri basis
-      REAL(KIND=dp), DIMENSION(:, :, :), POINTER                :: scon_ri
+      REAL(KIND=dp), DIMENSION(:, :, :), POINTER              :: scon_ri
       ! contraction matrix for SHG integrals orb basis
-      REAL(KIND=dp), DIMENSION(:, :, :), POINTER                :: scon_orb
+      REAL(KIND=dp), DIMENSION(:, :, :), POINTER              :: scon_orb
    END TYPE lri_bas_type
 
 ! **************************************************************************************************
@@ -159,7 +159,7 @@ MODULE lri_environment_types
       ! list of non-zero CG coefficients
       INTEGER, DIMENSION(:, :, :), POINTER                      :: cg_none0_list
       ! number of non-zero CG coefficients
-      INTEGER, DIMENSION(:, :), POINTER                        :: ncg_none0
+      INTEGER, DIMENSION(:, :), POINTER                         :: ncg_none0
    END TYPE lri_clebsch_gordon_type
 
 ! **************************************************************************************************

src/lri_shg_integrals.F

@@ -1,732 +0,0 @@
-!--------------------------------------------------------------------------------------------------!
-!   CP2K: A general program to perform molecular dynamics simulations                              !
-!   Copyright (C) 2000 - 2016  CP2K developers group                                               !
-!--------------------------------------------------------------------------------------------------!
-
-! **************************************************************************************************
-!> \brief Calculates integral matrices for LRIGPW method using solid harmonic
-!>        Gaussians (SHG)
-!>        lri : local resolution of the identity
-!> \par History
-!>      created [06.2015]
-!> \author Dorothea Golze
-! **************************************************************************************************
-MODULE lri_shg_integrals
-   USE ai_shg_overlap,                  ONLY: &
-        dev_overlap_shg_ab, dev_overlap_shg_aba, dev_overlap_shg_abb, get_W_matrix, get_dW_matrix, &
-        get_real_scaled_solid_harmonic, overlap_shg_ab, overlap_shg_aba, overlap_shg_abb, &
-        s_overlap_ab, s_overlap_abx
-   USE basis_set_types,                 ONLY: gto_basis_set_type
-   USE kinds,                           ONLY: dp
-   USE lri_environment_types,           ONLY: lri_environment_type
-   USE mathconstants,                   ONLY: fac,&
-                                              ifac,&
-                                              pi
-   USE orbital_pointers,                ONLY: nsoset
-#include "./base/base_uses.f90"
-
-   IMPLICIT NONE
-
-   PRIVATE
-
-! **************************************************************************************************
-
-   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'lri_shg_integrals'
-
-   PUBLIC :: precalc_angular_shg_part, lri_int_ab_shg, lri_int_aba_shg, &
-             lri_int_abb_shg
-
-! **************************************************************************************************
-
-CONTAINS
-
-! **************************************************************************************************
-!> \brief precalculates the angular part of the SHG integrals for the matrices
-!>        (fa,fb), (a,b), (a,b,fa) and (b,fb,a)
-!> \param oba orbital basis on a
-!> \param obb orbital basis on b
-!> \param fba aux basis on a
-!> \param fbb aux basis on b
-!> \param rab distance vector between a and b
-!> \param Waux_mat matrix storing angular-dependent part
-!> \param dWaux_mat ...
-!> \param calculate_forces ...
-! **************************************************************************************************
-   SUBROUTINE precalc_angular_shg_part(oba, obb, fba, fbb, rab, Waux_mat, dWaux_mat, &
-                                       calculate_forces)
-
-      TYPE(gto_basis_set_type), POINTER                  :: oba, obb, fba, fbb
-      REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rab
-      REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER      :: Waux_mat
-      REAL(KIND=dp), DIMENSION(:, :, :, :, :), POINTER   :: dWaux_mat
-      LOGICAL, INTENT(IN)                                :: calculate_forces
-
-      CHARACTER(len=*), PARAMETER :: routineN = 'precalc_angular_shg_part', &
-         routineP = moduleN//':'//routineN
-
-      INTEGER                                            :: handle, i, isize, j, k, la_max, laa_max, &
-                                                            lb_max, lbb_max, lca_max, lcb_max, &
-                                                            li_max, lj_max, lmax, mdim(4), &
-                                                            size_int(4, 2), temp
-      INTEGER, DIMENSION(:), POINTER                     :: li_max_all
-      REAL(KIND=dp)                                      :: rab2
-      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: Rc, Rs
-
-      CALL timeset(routineN, handle)
-      NULLIFY (li_max_all)
-
-      rab2 = rab(1)*rab(1)+rab(2)*rab(2)+rab(3)*rab(3)
-
-      !*** 1 Waux_mat of size (li_max,lj_max) for elements
-      !                    i        j
-      !    [aab]    --> (laa_max, lb_max)
-      !    [bba]    --> (lbb_max, la_max) --> use for [abb]
-      !    [ab] ri  --> (lca_max, lcb_max)
-      !    [ab] orb --> (la_max , lb_max)
-
-      la_max = MAXVAL(oba%lmax)
-      lb_max = MAXVAL(obb%lmax)
-      lca_max = MAXVAL(fba%lmax)
-      lcb_max = MAXVAL(fbb%lmax)
-
-      laa_max = la_max+lca_max
-      lbb_max = lb_max+lcb_max
-      li_max = MAX(laa_max, lbb_max)
-      lj_max = MAX(la_max, lb_max, lcb_max)
-      lmax = li_max
-
-      ALLOCATE (li_max_all(0:lj_max))
-      ALLOCATE (Rc(0:lmax, -2*lmax:2*lmax), Rs(0:lmax, -2*lmax:2*lmax))
-      Rc = 0._dp
-      Rs = 0._dp
-      mdim(1) = nsoset(li_max)+1
-      mdim(2) = nsoset(lj_max)+1
-      mdim(3) = li_max+lj_max+1
-      mdim(4) = 4
-      ALLOCATE (Waux_mat(mdim(1), mdim(2), mdim(3), mdim(4)))
-      ALLOCATE (dWaux_mat(3, mdim(1), mdim(2), mdim(3), mdim(4)))
-      !Waux_mat = 0._dp !.. takes time
-
-      !*** Waux_mat (li_max,lj_max) contains elements not needed,
-      !*** make indixing so that only required ones are computed
-      !*** li_max_all(j) --> li_max dependent on j
-      size_int(1, :) = (/laa_max, lb_max/)
-      size_int(2, :) = (/lbb_max, la_max/)
-      size_int(3, :) = (/lca_max, lcb_max/)
-      size_int(4, :) = (/la_max, lb_max/)
-
-      li_max_all(:) = 0
-      DO isize = 1, 4
-         i = size_int(isize, 1)
-         j = size_int(isize, 2)
-         k = li_max_all(j)
-         IF (k < i) li_max_all(j) = i
-      ENDDO
-      temp = li_max_all(lj_max)
-      DO j = lj_max, 0, -1
-         IF (li_max_all(j) < temp) THEN
-            li_max_all(j) = temp
-         ELSE
-            temp = li_max_all(j)
-         ENDIF
-      ENDDO
-
-      !*** -rab, since Eq. in Ref. use Ra-Rb, not Rb-Ra
-      CALL get_real_scaled_solid_harmonic(Rc, Rs, lmax, -rab, rab2)
-      CALL get_W_matrix(li_max_all, lj_max, lmax, Rc, Rs, Waux_mat)
-      IF (calculate_forces) THEN
-         CALL get_dW_matrix(li_max_all, lj_max, Waux_mat, dWaux_mat)
-      ENDIF
-
-      DEALLOCATE (Rc, Rs, li_max_all)
-
-      CALL timestop(handle)
-
-   END SUBROUTINE precalc_angular_shg_part
-
-! **************************************************************************************************
-!> \brief calculate integrals (a,b,fa)
-!> \param lri_env ...
-!> \param abaint integral (a,b,fa)
-!> \param dabdaint ...
-!> \param rab distance vector between A and B
-!> \param oba orbital basis at center A
-!> \param obb orbital basis at center B
-!> \param fba auxiliary basis set at center A
-!> \param scona_shg contraction matrix for orb bas on A
-!> \param sconb_shg contraction matrix for orb bas on B
-!> \param sconca_shg contraction matrix for aux bas on A
-!> \param Waux_mat angular-depedent part
-!> \param dWaux_mat derivative of Waux_mat
-!> \param calculate_forces ...
-! **************************************************************************************************
-   SUBROUTINE lri_int_aba_shg(lri_env, abaint, dabdaint, rab, oba, obb, fba, scona_shg, &
-                              sconb_shg, sconca_shg, Waux_mat, dWaux_mat, &
-                              calculate_forces)
-
-      TYPE(lri_environment_type), POINTER                :: lri_env
-      REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: abaint
-      REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER      :: dabdaint
-      REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rab
-      TYPE(gto_basis_set_type), POINTER                  :: oba, obb, fba
-      REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: scona_shg, sconb_shg, sconca_shg
-      REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER      :: Waux_mat
-      REAL(KIND=dp), DIMENSION(:, :, :, :, :), POINTER   :: dWaux_mat
-      LOGICAL, INTENT(IN)                                :: calculate_forces
-
-      CHARACTER(LEN=*), PARAMETER :: routineN = 'lri_int_aba_shg', &
-         routineP = moduleN//':'//routineN
-
-      INTEGER :: handle, ids, ids_start, iil, il, ipgfa, iset, ishella, jpgfb, jset, jshellb, &
-         kpgfa, kset, kshella, la_max_set, laa, laa_max_set, lai, lak, lb_max_set, lbj, &
-         lca_max_set, ndev, nds, nds_max, nl, nl_set, npgfa_set, npgfb_set, npgfca_set, nseta, &
-         nsetb, nsetca, nshella_set, nshellb_set, nshellca_set
-      INTEGER, DIMENSION(:), POINTER                     :: la_max, lb_max, lca_max, npgfa, npgfb, &
-                                                            npgfca, nshella, nshellb, nshellca
-      INTEGER, DIMENSION(:, :), POINTER                  :: first_sgfa, first_sgfb, first_sgfca, la, &
-                                                            lb, lca
-      REAL(KIND=dp)                                      :: a, aif, bjf, const, consta, constca, &
-                                                            dab, prefac, rab2, shg_fac(5)
-      REAL(KIND=dp), DIMENSION(:), POINTER               :: set_radius_a, set_radius_b, set_radius_ca
-      REAL(KIND=dp), DIMENSION(:, :), POINTER            :: zeta, zetb, zetca
-      REAL(KIND=dp), DIMENSION(:, :, :, :, :), POINTER   :: swork, swork_cont
-
-      CALL timeset(routineN, handle)
-      NULLIFY (la_max, lb_max, lca_max, npgfa, npgfb, npgfca)
-      NULLIFY (first_sgfa, first_sgfb, first_sgfca, set_radius_a, set_radius_b, &
-               set_radius_ca, zeta, zetb, zetca)
-
-      ! basis ikind
-      first_sgfa => oba%first_sgf
-      la_max => oba%lmax
-      la => oba%l
-      npgfa => oba%npgf
-      nshella => oba%nshell
-      nseta = oba%nset
-      set_radius_a => oba%set_radius
-      zeta => oba%zet
-      ! basis jkind
-      first_sgfb => obb%first_sgf
-      lb_max => obb%lmax
-      lb => obb%l
-      npgfb => obb%npgf
-      nshellb => obb%nshell
-      nsetb = obb%nset
-      set_radius_b => obb%set_radius
-      zetb => obb%zet
-
-      ! basis RI A
-      first_sgfca => fba%first_sgf
-      lca_max => fba%lmax
-      lca => fba%l
-      npgfca => fba%npgf
-      nshellca => fba%nshell
-      nsetca = fba%nset
-      set_radius_ca => fba%set_radius
-      zetca => fba%zet
-
-      dab = SQRT(SUM(rab**2))
-      rab2 = rab(1)*rab(1)+rab(2)*rab(2)+rab(3)*rab(3)
-
-      la_max_set = MAXVAL(la_max)
-      lb_max_set = MAXVAL(lb_max)
-      lca_max_set = MAXVAL(lca_max)
-      npgfa_set = MAXVAL(npgfa)
-      npgfb_set = MAXVAL(npgfb)
-      npgfca_set = MAXVAL(npgfca)
-      nshella_set = MAXVAL(nshella)
-      nshellb_set = MAXVAL(nshellb)
-      nshellca_set = MAXVAL(nshellca)
-      !*** for forces: derivative+1 in auxiliary vector required
-      ndev = 0
-      IF (calculate_forces) ndev = 1
-
-      laa_max_set = la_max_set+lca_max_set
-
-      ! allocate some work storage....
-      nds_max = laa_max_set+lb_max_set+ndev+1
-      nl_set = INT((laa_max_set)/2)
-      ALLOCATE (swork(npgfa_set, npgfb_set, npgfca_set, 0:nl_set, nds_max))
-      ALLOCATE (swork_cont(nshella_set, nshellb_set, nshellca_set, 0:nl_set, nds_max))
-
-      DO iset = 1, nseta
-
-         DO jset = 1, nsetb
-
-            IF (set_radius_a(iset)+set_radius_b(jset) < dab) CYCLE
-
-            DO kset = 1, nsetca
-
-               IF (set_radius_b(jset)+set_radius_ca(kset) < dab) CYCLE
-
-               CALL s_overlap_abx(la_max(iset), npgfa(iset), zeta(:, iset), &
-                                  lb_max(jset), npgfb(jset), zetb(:, jset), &
-                                  lca_max(kset), npgfca(kset), zetca(:, kset), &
-                                  rab, swork, calculate_forces, calc_aba=.TRUE.)
-
-               swork_cont = 0.0_dp
-               DO ishella = 1, nshella(iset)
-                  lai = la(ishella, iset)
-                  consta = SQRT(1.0_dp/REAL(2*lai+1, dp))
-                  DO jshellb = 1, nshellb(jset)
-                     lbj = lb(jshellb, jset)
-                     DO kshella = 1, nshellca(kset)
-                        lak = lca(kshella, kset)
-                        constca = SQRT(1.0_dp/REAL(2*lak+1, dp))
-                        nl = INT((lai+lak)/2)
-                        IF (lai == 0 .OR. lak == 0) nl = 0
-                        DO il = 0, nl
-                           laa = lai+lak-2*il
-                           nds = laa+lbj+1
-                           ids_start = nds-MIN(laa, lbj)
-                           const = consta*constca*2.0_dp*SQRT(pi*REAL(2*laa+1, dp))
-                           DO iil = 1, il
-                              shg_fac(iil) = fac(laa+iil-1)*ifac(laa)*REAL(laa, dp) &
-                                             *fac(il)/fac(il-iil)/fac(iil)
-                           ENDDO
-                           DO ipgfa = 1, npgfa(iset)
-                              DO jpgfb = 1, npgfb(jset)
-                                 bjf = 1.0_dp/((2._dp*zetb(jpgfb, jset))**lbj)
-                                 DO kpgfa = 1, npgfca(kset)
-                                    a = zeta(ipgfa, iset)+zetca(kpgfa, kset)
-                                    aif = 1.0_dp/((2._dp*a)**laa)
-                                    prefac = const*aif*bjf &
-                                             *scona_shg(ipgfa, ishella, iset) &
-                                             *sconb_shg(jpgfb, jshellb, jset) &
-                                             *sconca_shg(kpgfa, kshella, kset)
-                                    DO ids = ids_start, nds+ndev
-                                       IF (il == 0) THEN
-                                          swork_cont(ishella, jshellb, kshella, il, ids) = &
-                                             swork_cont(ishella, jshellb, kshella, il, ids) &
-                                             +prefac*swork(ipgfa, jpgfb, kpgfa, il, ids)
-                                       ELSEIF (il == 1) THEN
-                                          swork_cont(ishella, jshellb, kshella, il, ids) = &
-                                             swork_cont(ishella, jshellb, kshella, il, ids) &
-                                             +prefac*(swork(ipgfa, jpgfb, kpgfa, il, ids) &
-                                                      +shg_fac(il)*swork(ipgfa, jpgfb, kpgfa, il-1, ids)/a)
-                                       ELSEIF (il == 2) THEN
-                                          swork_cont(ishella, jshellb, kshella, il, ids) = &
-                                             swork_cont(ishella, jshellb, kshella, il, ids) &
-                                             +prefac*(swork(ipgfa, jpgfb, kpgfa, il, ids) &
-                                                      +shg_fac(il-1)*swork(ipgfa, jpgfb, kpgfa, il-1, ids)/a &
-                                                      +shg_fac(il)*swork(ipgfa, jpgfb, kpgfa, il-2, ids)/(a**2))
-                                       ELSEIF (il == 3) THEN
-                                          swork_cont(ishella, jshellb, kshella, il, ids) = &
-                                             swork_cont(ishella, jshellb, kshella, il, ids) &
-                                             +prefac*(swork(ipgfa, jpgfb, kpgfa, il, ids) &
-                                                      +shg_fac(il-2)*swork(ipgfa, jpgfb, kpgfa, il-1, ids)/a &
-                                                      +shg_fac(il-1)*swork(ipgfa, jpgfb, kpgfa, il-2, ids)/(a**2) &
-                                                      +shg_fac(il)*swork(ipgfa, jpgfb, kpgfa, il-3, ids)/(a**3))
-                                       ELSEIF (il == 4) THEN
-                                          swork_cont(ishella, jshellb, kshella, il, ids) = &
-                                             swork_cont(ishella, jshellb, kshella, il, ids) &
-                                             +prefac*(swork(ipgfa, jpgfb, kpgfa, il, ids) &
-                                                      +shg_fac(il-3)*swork(ipgfa, jpgfb, kpgfa, il-1, ids)/a &
-                                                      +shg_fac(il-2)*swork(ipgfa, jpgfb, kpgfa, il-2, ids)/(a**2) &
-                                                      +shg_fac(il-1)*swork(ipgfa, jpgfb, kpgfa, il-3, ids)/(a**3) &
-                                                      +shg_fac(il)*swork(ipgfa, jpgfb, kpgfa, il-4, ids)/(a**4))
-                                       ELSEIF (il == 5) THEN
-                                          swork_cont(ishella, jshellb, kshella, il, ids) = &
-                                             swork_cont(ishella, jshellb, kshella, il, ids) &
-                                             +prefac*(swork(ipgfa, jpgfb, kpgfa, il, ids) &
-                                                      +shg_fac(il-4)*swork(ipgfa, jpgfb, kpgfa, il-1, ids)/a &
-                                                      +shg_fac(il-3)*swork(ipgfa, jpgfb, kpgfa, il-2, ids)/a**2 &
-                                                      +shg_fac(il-2)*swork(ipgfa, jpgfb, kpgfa, il-3, ids)/(a**3) &
-                                                      +shg_fac(il-1)*swork(ipgfa, jpgfb, kpgfa, il-4, ids)/(a**4) &
-                                                      +shg_fac(il)*swork(ipgfa, jpgfb, kpgfa, il-5, ids)/(a**5))
-                                       ELSE
-                                          CALL cp_abort(__LOCATION__, &
-                                                        "SHG integrals not implemented when l quantum number"// &
-                                                        " of orbital and ri basis larger than 11")
-                                       ENDIF
-                                    ENDDO
-                                 ENDDO
-                              ENDDO
-                           ENDDO
-                        ENDDO
-                     ENDDO
-                  ENDDO
-               ENDDO
-               CALL overlap_shg_aba(la(:, iset), first_sgfa(:, iset), nshella(iset), &
-                                    lb(:, jset), first_sgfb(:, jset), nshellb(jset), &
-                                    lca(:, kset), first_sgfca(:, kset), nshellca(kset), &
-                                    lri_env%cg_shg%cg_coeff, lri_env%cg_shg%cg_none0_list, &
-                                    lri_env%cg_shg%ncg_none0, swork_cont, Waux_mat, abaint)
-               !
-               IF (calculate_forces) THEN
-                  !*** -rab, since Eq. in Ref. use Ra-Rb, not Rb-Ra
-                  CALL dev_overlap_shg_aba(la(:, iset), first_sgfa(:, iset), nshella(iset), &
-                                           lb(:, jset), first_sgfb(:, jset), nshellb(jset), &
-                                           lca(:, kset), first_sgfca(:, kset), nshellca(kset), &
-                                           lri_env%cg_shg%cg_coeff, lri_env%cg_shg%cg_none0_list, &
-                                           lri_env%cg_shg%ncg_none0, -rab, swork_cont, &
-                                           Waux_mat, dWaux_mat, dabdaint)
-               ENDIF
-
-            END DO
-         END DO
-      END DO
-
-      DEALLOCATE (swork_cont)
-      DEALLOCATE (swork)
-
-      CALL timestop(handle)
-
-   END SUBROUTINE lri_int_aba_shg
-
-! **************************************************************************************************
-!> \brief calculate integrals (a,b,fb)
-!> \param lri_env ...
-!> \param abbint integral (a,b,fb)
-!> \param dabbint derivative of abbint
-!> \param rab distance vector between A and B
-!> \param oba orbital basis at center A
-!> \param obb orbital basis at center B
-!> \param fbb auxiliary basis set at center B
-!> \param scona_shg contraction matrix for orb bas on A
-!> \param sconb_shg contraction matrix for orb bas on B
-!> \param sconcb_shg contraction matrix for aux bas on B
-!> \param Waux_mat angular-dependent part
-!> \param dWaux_mat derivative of Waux_mat
-!> \param calculate_forces ...
-! **************************************************************************************************
-   SUBROUTINE lri_int_abb_shg(lri_env, abbint, dabbint, rab, oba, obb, fbb, scona_shg, &
-                              sconb_shg, sconcb_shg, Waux_mat, dWaux_mat, &
-                              calculate_forces)
-
-      TYPE(lri_environment_type), POINTER                :: lri_env
-      REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: abbint
-      REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER      :: dabbint
-      REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rab
-      TYPE(gto_basis_set_type), POINTER                  :: oba, obb, fbb
-      REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: scona_shg, sconb_shg, sconcb_shg
-      REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER      :: Waux_mat
-      REAL(KIND=dp), DIMENSION(:, :, :, :, :), POINTER   :: dWaux_mat
-      LOGICAL, INTENT(IN)                                :: calculate_forces
-
-      CHARACTER(LEN=*), PARAMETER :: routineN = 'lri_int_abb_shg', &
-         routineP = moduleN//':'//routineN
-
-      INTEGER :: handle, ids, ids_start, iil, il, ipgfa, iset, ishella, jpgfb, jset, jshellb, &
-         kpgfb, kset, kshellb, la_max_set, lai, lb_max_set, lbb, lbb_max_set, lbj, lbk, &
-         lcb_max_set, ndev, nds, nds_max, nl, nl_set, npgfa_set, npgfb_set, npgfcb_set, nseta, &
-         nsetb, nsetcb, nshella_set, nshellb_set, nshellcb_set
-      INTEGER, DIMENSION(:), POINTER                     :: la_max, lb_max, lcb_max, npgfa, npgfb, &
-                                                            npgfcb, nshella, nshellb, nshellcb
-      INTEGER, DIMENSION(:, :), POINTER                  :: first_sgfa, first_sgfb, first_sgfcb, la, &
-                                                            lb, lcb
-      REAL(KIND=dp)                                      :: aif, b, bjf, const, constb, constcb, &
-                                                            dab, prefac, rab2, shg_fac(5)
-      REAL(KIND=dp), DIMENSION(:), POINTER               :: set_radius_a, set_radius_b, set_radius_cb
-      REAL(KIND=dp), DIMENSION(:, :), POINTER            :: zeta, zetb, zetcb
-      REAL(KIND=dp), DIMENSION(:, :, :, :, :), POINTER   :: swork, swork_cont
-
-      CALL timeset(routineN, handle)
-      NULLIFY (la_max, lb_max, lcb_max, npgfa, npgfb, npgfcb)
-      NULLIFY (first_sgfa, first_sgfb, first_sgfcb, set_radius_a, set_radius_b, &
-               set_radius_cb, zeta, zetb, zetcb)
-
-      ! basis ikind
-      first_sgfa => oba%first_sgf
-      la_max => oba%lmax
-      la => oba%l
-      npgfa => oba%npgf
-      nshella => oba%nshell
-      nseta = oba%nset
-      set_radius_a => oba%set_radius
-      zeta => oba%zet
-      ! basis jkind
-      first_sgfb => obb%first_sgf
-      lb_max => obb%lmax
-      lb => obb%l
-      npgfb => obb%npgf
-      nshellb => obb%nshell
-      nsetb = obb%nset
-      set_radius_b => obb%set_radius
-      zetb => obb%zet
-
-      ! basis RI on B
-      first_sgfcb => fbb%first_sgf
-      lcb_max => fbb%lmax
-      lcb => fbb%l
-      npgfcb => fbb%npgf
-      nshellcb => fbb%nshell
-      nsetcb = fbb%nset
-      set_radius_cb => fbb%set_radius
-      zetcb => fbb%zet
-
-      dab = SQRT(SUM(rab**2))
-      rab2 = rab(1)*rab(1)+rab(2)*rab(2)+rab(3)*rab(3)
-
-      la_max_set = MAXVAL(la_max)
-      lb_max_set = MAXVAL(lb_max)
-      lcb_max_set = MAXVAL(lcb_max)
-      npgfa_set = MAXVAL(npgfa)
-      npgfb_set = MAXVAL(npgfb)
-      npgfcb_set = MAXVAL(npgfcb)
-      nshella_set = MAXVAL(nshella)
-      nshellb_set = MAXVAL(nshellb)
-      nshellcb_set = MAXVAL(nshellcb)
-      !*** for forces: derivative+1 in auxiliary vector required
-      ndev = 0
-      IF (calculate_forces) ndev = 1
-
-      lbb_max_set = lb_max_set+lcb_max_set
-
-      ! allocate some work storage....
-      nds_max = la_max_set+lbb_max_set+ndev+1
-      nl_set = INT((lbb_max_set)/2)
-      ALLOCATE (swork(1:npgfa_set, 1:npgfb_set, 1:npgfcb_set, 0:nl_set, 1:nds_max))
-      ALLOCATE (swork_cont(nshella_set, nshellb_set, nshellcb_set, 0:nl_set, nds_max))
-
-      DO iset = 1, nseta
-
-         DO jset = 1, nsetb
-
-            IF (set_radius_a(iset)+set_radius_b(jset) < dab) CYCLE
-
-            DO kset = 1, nsetcb
-
-               IF (set_radius_a(iset)+set_radius_cb(kset) < dab) CYCLE
-
-               CALL s_overlap_abx(la_max(iset), npgfa(iset), zeta(:, iset), &
-                                  lb_max(jset), npgfb(jset), zetb(:, jset), &
-                                  lcb_max(kset), npgfcb(kset), zetcb(:, kset), &
-                                  rab, swork, calculate_forces, calc_aba=.FALSE.)
-
-               swork_cont = 0.0_dp
-               DO ishella = 1, nshella(iset)
-                  lai = la(ishella, iset)
-                  DO jshellb = 1, nshellb(jset)
-                     lbj = lb(jshellb, jset)
-                     constb = SQRT(1.0_dp/REAL(2*lbj+1, dp))
-                     DO kshellb = 1, nshellcb(kset)
-                        lbk = lcb(kshellb, kset)
-                        constcb = SQRT(1.0_dp/REAL(2*lbk+1, dp))
-                        nl = INT((lbj+lbk)/2)
-                        IF (lbj == 0 .OR. lbk == 0) nl = 0
-                        DO il = 0, nl
-                           lbb = lbj+lbk-2*il
-                           nds = lai+lbb+1
-                           ids_start = nds-MIN(lai, lbb)
-                           const = constb*constcb*2.0_dp*SQRT(pi*REAL(2*lbb+1, dp))
-                           DO iil = 1, il
-                              shg_fac(iil) = fac(lbb+iil-1)*ifac(lbb)*REAL(lbb, dp) &
-                                             *fac(il)/fac(il-iil)/fac(iil)
-                           ENDDO
-                           DO ipgfa = 1, npgfa(iset)
-                              aif = 1.0_dp/((2._dp*zeta(ipgfa, iset))**lai)
-                              DO jpgfb = 1, npgfb(jset)
-                                 DO kpgfb = 1, npgfcb(kset)
-                                    b = zetb(jpgfb, jset)+zetcb(kpgfb, kset)
-                                    bjf = 1.0_dp/((2._dp*b)**lbb)
-                                    prefac = const*aif*bjf &
-                                             *scona_shg(ipgfa, ishella, iset) &
-                                             *sconb_shg(jpgfb, jshellb, jset) &
-                                             *sconcb_shg(kpgfb, kshellb, kset)
-                                    DO ids = ids_start, nds+ndev
-                                       IF (il == 0) THEN
-                                          swork_cont(ishella, jshellb, kshellb, il, ids) = &
-                                             swork_cont(ishella, jshellb, kshellb, il, ids) &
-                                             +prefac*swork(ipgfa, jpgfb, kpgfb, il, ids)
-                                       ELSEIF (il == 1) THEN
-                                          swork_cont(ishella, jshellb, kshellb, il, ids) = &
-                                             swork_cont(ishella, jshellb, kshellb, il, ids) &
-                                             +prefac*(swork(ipgfa, jpgfb, kpgfb, il, ids) &
-                                                      +shg_fac(il)*swork(ipgfa, jpgfb, kpgfb, il-1, ids)/b)
-                                       ELSEIF (il == 2) THEN
-                                          swork_cont(ishella, jshellb, kshellb, il, ids) = &
-                                             swork_cont(ishella, jshellb, kshellb, il, ids) &
-                                             +prefac*(swork(ipgfa, jpgfb, kpgfb, il, ids) &
-                                                      +shg_fac(il-1)*swork(ipgfa, jpgfb, kpgfb, il-1, ids)/b &
-                                                      +shg_fac(il)*swork(ipgfa, jpgfb, kpgfb, il-2, ids)/(b**2))
-                                       ELSEIF (il == 3) THEN
-                                          swork_cont(ishella, jshellb, kshellb, il, ids) = &
-                                             swork_cont(ishella, jshellb, kshellb, il, ids) &
-                                             +prefac*(swork(ipgfa, jpgfb, kpgfb, il, ids) &
-                                                      +shg_fac(il-2)*swork(ipgfa, jpgfb, kpgfb, il-1, ids)/b &
-                                                      +shg_fac(il-1)*swork(ipgfa, jpgfb, kpgfb, il-2, ids)/(b**2) &
-                                                      +shg_fac(il)*swork(ipgfa, jpgfb, kpgfb, il-3, ids)/(b**3))
-                                       ELSEIF (il == 4) THEN
-                                          swork_cont(ishella, jshellb, kshellb, il, ids) = &
-                                             swork_cont(ishella, jshellb, kshellb, il, ids) &
-                                             +prefac*(swork(ipgfa, jpgfb, kpgfb, il, ids) &
-                                                      +shg_fac(il-3)*swork(ipgfa, jpgfb, kpgfb, il-1, ids)/b &
-                                                      +shg_fac(il-2)*swork(ipgfa, jpgfb, kpgfb, il-2, ids)/(b**2) &
-                                                      +shg_fac(il-1)*swork(ipgfa, jpgfb, kpgfb, il-3, ids)/(b**3) &
-                                                      +shg_fac(il)*swork(ipgfa, jpgfb, kpgfb, il-4, ids)/(b**4))
-                                       ELSEIF (il == 5) THEN
-                                          swork_cont(ishella, jshellb, kshellb, il, ids) = &
-                                             swork_cont(ishella, jshellb, kshellb, il, ids) &
-                                             +prefac*(swork(ipgfa, jpgfb, kpgfb, il, ids) &
-                                                      +shg_fac(il-4)*swork(ipgfa, jpgfb, kpgfb, il-1, ids)/b &
-                                                      +shg_fac(il-3)*swork(ipgfa, jpgfb, kpgfb, il-2, ids)/b**2 &
-                                                      +shg_fac(il-2)*swork(ipgfa, jpgfb, kpgfb, il-3, ids)/(b**3) &
-                                                      +shg_fac(il-1)*swork(ipgfa, jpgfb, kpgfb, il-4, ids)/(b**4) &
-                                                      +shg_fac(il)*swork(ipgfa, jpgfb, kpgfb, il-5, ids)/(b**5))
-                                       ELSE
-                                          CALL cp_abort(__LOCATION__, &
-                                                        "SHG integrals not implemented when l quantum number"// &
-                                                        " of orbital and ri basis larger than 11")
-                                       ENDIF
-                                    ENDDO
-                                 ENDDO
-                              ENDDO
-                           ENDDO
-                        ENDDO
-                     ENDDO
-                  ENDDO
-               ENDDO
-               CALL overlap_shg_abb(la(:, iset), first_sgfa(:, iset), nshella(iset), &
-                                    lb(:, jset), first_sgfb(:, jset), nshellb(jset), &
-                                    lcb(:, kset), first_sgfcb(:, kset), nshellcb(kset), &
-                                    lri_env%cg_shg%cg_coeff, lri_env%cg_shg%cg_none0_list, &
-                                    lri_env%cg_shg%ncg_none0, swork_cont, Waux_mat, abbint)
-               IF (calculate_forces) THEN
-                  !*** -rab, since Eq. in Ref. use Ra-Rb, not Rb-Ra
-                  CALL dev_overlap_shg_abb(la(:, iset), first_sgfa(:, iset), nshella(iset), &
-                                           lb(:, jset), first_sgfb(:, jset), nshellb(jset), &
-                                           lcb(:, kset), first_sgfcb(:, kset), nshellcb(kset), &
-                                           lri_env%cg_shg%cg_coeff, lri_env%cg_shg%cg_none0_list, &
-                                           lri_env%cg_shg%ncg_none0, -rab, swork_cont, &
-                                           Waux_mat, dWaux_mat, dabbint)
-               ENDIF
-
-            END DO
-         END DO
-      END DO
-
-      DEALLOCATE (swork_cont)
-      DEALLOCATE (swork)
-
-      CALL timestop(handle)
-
-   END SUBROUTINE lri_int_abb_shg
-
-! **************************************************************************************************
-!> \brief calculate overlap integrals (a,b)
-!> \param sab integral (a,b)
-!> \param dsab derivative of sab
-!> \param rab distance vector
-!> \param fba basis at center A
-!> \param fbb basis at center B
-!> \param scona_shg contraction matrix A
-!> \param sconb_shg contraxtion matrix B
-!> \param Waux_mat angular-dependent part
-!> \param dWaux_mat derivative of Waux_mat
-!> \param calculate_forces ...
-! **************************************************************************************************
-   SUBROUTINE lri_int_ab_shg(sab, dsab, rab, fba, fbb, scona_shg, sconb_shg, Waux_mat, &
-                             dWaux_mat, calculate_forces)
-
-      REAL(KIND=dp), DIMENSION(:, :), POINTER            :: sab
-      REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: dsab
-      REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rab
-      TYPE(gto_basis_set_type), POINTER                  :: fba, fbb
-      REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: scona_shg, sconb_shg
-      REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER      :: Waux_mat
-      REAL(KIND=dp), DIMENSION(:, :, :, :, :), POINTER   :: dWaux_mat
-      LOGICAL, INTENT(IN)                                :: calculate_forces
-
-      CHARACTER(LEN=*), PARAMETER :: routineN = 'lri_int_ab_shg', routineP = moduleN//':'//routineN
-
-      INTEGER :: handle, ids, ids_start, ipgfa, iset, ishella, jpgfb, jset, jshellb, la_max_set, &
-         lai, lb_max_set, lbj, ndev, nds, nds_max, npgfa_set, npgfb_set, nseta, nsetb, &
-         nshella_set, nshellb_set
-      INTEGER, DIMENSION(:), POINTER                     :: la_max, lb_max, npgfa, npgfb, nshella, &
-                                                            nshellb
-      INTEGER, DIMENSION(:, :), POINTER                  :: first_sgfa, first_sgfb, la, lb
-      REAL(KIND=dp)                                      :: aif, bjf, dab
-      REAL(KIND=dp), DIMENSION(:), POINTER               :: set_radius_a, set_radius_b
-      REAL(KIND=dp), DIMENSION(:, :), POINTER            :: zeta, zetb
-      REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: swork, swork_cont
-
-      CALL timeset(routineN, handle)
-      NULLIFY (la_max, lb_max, npgfa, npgfb, first_sgfa, first_sgfb, set_radius_a, &
-               set_radius_b, zeta, zetb)
-
-      ! basis ikind
-      first_sgfa => fba%first_sgf
-      la_max => fba%lmax
-      la => fba%l
-      npgfa => fba%npgf
-      nseta = fba%nset
-      set_radius_a => fba%set_radius
-      zeta => fba%zet
-      nshella => fba%nshell
-      ! basis jkind
-      first_sgfb => fbb%first_sgf
-      lb_max => fbb%lmax
-      lb => fbb%l
-      npgfb => fbb%npgf
-      nsetb = fbb%nset
-      set_radius_b => fbb%set_radius
-      zetb => fbb%zet
-      nshellb => fbb%nshell
-
-      dab = SQRT(SUM(rab**2))
-
-      la_max_set = MAXVAL(la_max)
-      lb_max_set = MAXVAL(lb_max)
-
-      ! allocate some work matrices
-      npgfa_set = MAXVAL(npgfa)
-      npgfb_set = MAXVAL(npgfb)
-      nshella_set = MAXVAL(nshella)
-      nshellb_set = MAXVAL(nshellb)
-      ndev = 0
-      IF (calculate_forces) ndev = 1
-      nds_max = la_max_set+lb_max_set+ndev+1
-      ALLOCATE (swork(npgfa_set, npgfb_set, nds_max))
-      ALLOCATE (swork_cont(nshella_set, nshellb_set, nds_max))
-
-      DO iset = 1, nseta
-
-         DO jset = 1, nsetb
-
-            IF (set_radius_a(iset)+set_radius_b(jset) < dab) CYCLE
-
-            CALL s_overlap_ab(la_max(iset), npgfa(iset), zeta(:, iset), &
-                              lb_max(jset), npgfb(jset), zetb(:, jset), &
-                              rab, swork, calculate_forces)
-
-            swork_cont = 0.0_dp
-            DO ishella = 1, nshella(iset)
-               lai = la(ishella, iset)
-               DO jshellb = 1, nshellb(jset)
-                  lbj = lb(jshellb, jset)
-                  nds = lai+lbj+1
-                  ids_start = nds-MIN(lai, lbj)
-                  DO ipgfa = 1, npgfa(iset)
-                     aif = 1.0_dp/((2._dp*zeta(ipgfa, iset))**lai)
-                     DO jpgfb = 1, npgfb(jset)
-                        bjf = 1.0_dp/((2._dp*zetb(jpgfb, jset))**lbj)
-                        DO ids = ids_start, nds+ndev
-                           swork_cont(ishella, jshellb, ids) = swork_cont(ishella, jshellb, ids) &
-                                                               +scona_shg(ipgfa, ishella, iset) &
-                                                               *sconb_shg(jpgfb, jshellb, jset) &
-                                                               *aif*bjf*swork(ipgfa, jpgfb, ids)
-                        ENDDO
-                     ENDDO
-                  ENDDO
-               ENDDO
-            ENDDO
-            CALL overlap_shg_ab(la(:, iset), first_sgfa(:, iset), nshella(iset), &
-                                lb(:, jset), first_sgfb(:, jset), nshellb(jset), &
-                                swork_cont, Waux_mat, sab)
-            IF (calculate_forces) THEN
-               !*** -rab, since Eq. in Ref. use Ra-Rb, not Rb-Ra
-               CALL dev_overlap_shg_ab(la(:, iset), first_sgfa(:, iset), nshella(iset), &
-                                       lb(:, jset), first_sgfb(:, jset), nshellb(jset), &
-                                       -rab, swork_cont, Waux_mat, dWaux_mat, dsab)
-            ENDIF
-         END DO
-      END DO
-
-      DEALLOCATE (swork, swork_cont)
-
-      CALL timestop(handle)
-
-   END SUBROUTINE lri_int_ab_shg
-
-END MODULE lri_shg_integrals

src/shg_integrals_test.F

@@ -0,0 +1,423 @@
+!--------------------------------------------------------------------------------------------------!
+!   CP2K: A general program to perform molecular dynamics simulations                              !
+!   Copyright (C) 2000 - 2016  CP2K developers group                                               !
+!--------------------------------------------------------------------------------------------------!
+
+! **************************************************************************************************
+!> \brief Calculates 2-center integrals for different r12 operators comparing the Solid harmonic
+!>        Gaussian integral scheme to the Obara-Saika (OS) scheme
+!> \author  Dorothea Golze [05.2016]
+! **************************************************************************************************
+MODULE shg_integrals_test
+
+   USE basis_set_types,                 ONLY: allocate_gto_basis_set,&
+                                              deallocate_gto_basis_set,&
+                                              gto_basis_set_type,&
+                                              init_orb_basis_set,&
+                                              read_gto_basis_set
+   USE generic_os_integrals,            ONLY: two_center_os_integrals
+   USE generic_shg_integrals,           ONLY: two_center_shg_integrals
+   USE input_constants,                 ONLY: operator_coulomb,&
+                                              operator_gauss,&
+                                              operator_verf,&
+                                              operator_verfc,&
+                                              operator_vgauss
+   USE input_cp2k_subsys,               ONLY: create_basis_section
+   USE input_keyword_types,             ONLY: keyword_create,&
+                                              keyword_release,&
+                                              keyword_type
+   USE input_section_types,             ONLY: &
+        section_add_keyword, section_add_subsection, section_create, section_release, &
+        section_type, section_vals_get, section_vals_get_subs_vals, section_vals_type, &
+        section_vals_val_get
+   USE input_val_types,                 ONLY: real_t
+   USE kinds,                           ONLY: default_string_length,&
+                                              dp
+   USE mathconstants,                   ONLY: fac,&
+                                              pi
+   USE orbital_pointers,                ONLY: init_orbital_pointers
+   USE orbital_transformation_matrices, ONLY: init_spherical_harmonics
+#include "./base/base_uses.f90"
+
+   IMPLICIT NONE
+
+   PRIVATE
+
+! **************************************************************************************************
+
+   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'shg_integrals_test'
+
+   PUBLIC :: create_shg_integrals_test_section, shg_integrals_perf_acc_test
+
+CONTAINS
+
+! **************************************************************************************************
+!> \brief Create input section for unit testing
+!> \param section ...
+! **************************************************************************************************
+   SUBROUTINE create_shg_integrals_test_section(section)
+      TYPE(section_type), INTENT(INOUT), POINTER         :: section
+
+      CHARACTER(len=*), PARAMETER :: routineN = 'create_shg_integrals_test_section', &
+         routineP = moduleN//':'//routineN
+
+      TYPE(keyword_type), POINTER                        :: keyword
+      TYPE(section_type), POINTER                        :: subsection
+
+      NULLIFY (keyword, subsection)
+
+      CPASSERT(.NOT. ASSOCIATED(section))
+      CALL section_create(section, "SHG_INTEGRALS_TEST", &
+                          description="Parameters for testing the SHG 2-center integrals for "// &
+                          "different r12 operators. Test w.r.t. performance and accurarcy.", &
+                          n_keywords=4, n_subsections=1)
+
+      CALL create_basis_section(subsection)
+      CALL section_add_subsection(section, subsection)
+      CALL section_release(subsection)
+
+      CALL keyword_create(keyword, &
+                          name="_SECTION_PARAMETERS_", &
+                          description="Controls the activation the SHG integral test. ", &
+                          default_l_val=.FALSE., &
+                          lone_keyword_l_val=.TRUE.)
+      CALL section_add_keyword(section, keyword)
+      CALL keyword_release(keyword)
+
+      CALL keyword_create(keyword, name="ABC", &
+                          description="Specify the lengths of the cell vectors A, B, and C. ", &
+                          usage="ABC 10.000 10.000 10.000", unit_str="angstrom", &
+                          n_var=3, type_of_var=real_t)
+      CALL section_add_keyword(section, keyword)
+      CALL keyword_release(keyword)
+
+      CALL keyword_create(keyword, name="NAB_MIN", &
+                          description="Minimum number of atomic distances to consider. ", &
+                          default_i_val=8)
+      CALL section_add_keyword(section, keyword)
+      CALL keyword_release(keyword)
+
+      CALL keyword_create(keyword, name="NREP", &
+                          description="Number of repeated calculation of each integral. ", &
+                          default_i_val=1)
+      CALL section_add_keyword(section, keyword)
+      CALL keyword_release(keyword)
+
+      CALL keyword_create(keyword, name="CHECK_ACCURACY", &
+                          description="Causes abortion when SHG and OS integrals differ "// &
+                          "more than 1.0E-08.", &
+                          default_l_val=.TRUE., lone_keyword_l_val=.TRUE.)
+      CALL section_add_keyword(section, keyword)
+      CALL keyword_release(keyword)
+
+   END SUBROUTINE create_shg_integrals_test_section
+
+! **************************************************************************************************
+!> \brief Unit test for performance and accuracy of the SHG integrals
+!> \param iw output unit
+!> \param shg_integrals_test_section ...
+! **************************************************************************************************
+   SUBROUTINE shg_integrals_perf_acc_test(iw, shg_integrals_test_section)
+      INTEGER, INTENT(IN)                                :: iw
+      TYPE(section_vals_type), INTENT(INOUT), POINTER    :: shg_integrals_test_section
+
+      CHARACTER(len=*), PARAMETER :: routineN = 'shg_integrals_perf_acc_test', &
+         routineP = moduleN//':'//routineN
+
+      CHARACTER(LEN=default_string_length)               :: basis_type
+      INTEGER                                            :: count_ab, iab, irep, jab, kab, lamax, &
+                                                            lbmax, lmax, nab, nab_min, nab_xyz, &
+                                                            nfa, nfb, nrep, nrep_bas
+      LOGICAL                                            :: acc_check
+      REAL(KIND=dp) :: acc_param, ddmax_coulomb, ddmax_gauss, ddmax_verf, ddmax_verfc, &
+         ddmax_vgauss, ddtemp, dmax_coulomb, dmax_gauss, dmax_verf, dmax_verfc, dmax_vgauss, &
+         dtemp, omega
+      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: rab
+      REAL(KIND=dp), DIMENSION(:), POINTER               :: cell_par
+      REAL(KIND=dp), DIMENSION(:, :), POINTER            :: vab_os, vab_shg
+      REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: dvab_os, dvab_shg, scona_shg, sconb_shg
+      TYPE(gto_basis_set_type), POINTER                  :: fba, fbb
+      TYPE(section_vals_type), POINTER                   :: basis_section
+
+      NULLIFY (fba, fbb, basis_section, cell_par)
+      CALL section_vals_val_get(shg_integrals_test_section, "ABC", r_vals=cell_par)
+      CALL section_vals_val_get(shg_integrals_test_section, "NAB_MIN", i_val=nab_min)
+      CALL section_vals_val_get(shg_integrals_test_section, "NREP", i_val=nrep)
+      CALL section_vals_val_get(shg_integrals_test_section, "CHECK_ACCURACY", l_val=acc_check)
+
+      !*** Read the basis set information
+      basis_section => section_vals_get_subs_vals(shg_integrals_test_section, "BASIS")
+      CALL section_vals_get(basis_section, n_repetition=nrep_bas)
+      IF (nrep_bas /= 2) RETURN
+      CALL allocate_gto_basis_set(fba)
+      CALL read_gto_basis_set(TRIM("A"), basis_type, fba, basis_section, irep=1)
+      lamax = MAXVAL(fba%lmax)
+      CALL allocate_gto_basis_set(fbb)
+      CALL read_gto_basis_set(TRIM("B"), basis_type, fbb, basis_section, irep=2)
+      lbmax = MAXVAL(fbb%lmax)
+      lmax = MAX(lamax, lbmax)
+
+      !*** Initialize basis set information
+      CALL init_orbital_pointers(lmax+1)
+      CALL init_spherical_harmonics(lmax, output_unit=-100)
+      fba%norm_type = 2
+      CALL init_orb_basis_set(fba)
+      fbb%norm_type = 2
+      CALL init_orb_basis_set(fbb)
+      ! if shg integrals are later actually used in the code, contraction_matrix_shg should be
+      ! moved to init_orb_basis_set and scon_shg should become an element of gto_basis_set_type
+      CALL contraction_matrix_shg(fba, scona_shg)
+      CALL contraction_matrix_shg(fbb, sconb_shg)
+
+      !*** Create range of rab (atomic distances) to be tested
+      nab_xyz = CEILING(REAL(nab_min, KIND=dp)**(1.0_dp/3.0_dp)-1.0E-06)
+      nab = nab_xyz**3
+
+      ALLOCATE (rab(3, nab))
+      count_ab = 0
+      DO iab = 1, nab_xyz
+         DO jab = 1, nab_xyz
+            DO kab = 1, nab_xyz
+               count_ab = count_ab+1
+               rab(:, count_ab) = [iab*ABS(cell_par(1)), jab*ABS(cell_par(2)), kab*ABS(cell_par(3))]/nab_xyz
+            ENDDO
+         ENDDO
+      ENDDO
+
+      !*** Calculate the SHG integrals
+      nfa = fba%nsgf
+      nfb = fbb%nsgf
+      ALLOCATE (vab_shg(nfa, nfb), dvab_shg(nfa, nfb, 3))
+      ALLOCATE (vab_os(nfa, nfb), dvab_os(nfa, nfb, 3))
+      omega = 2.3_dp
+      dmax_coulomb = 0.0_dp
+      ddmax_coulomb = 0.0_dp
+      dmax_verf = 0.0_dp
+      ddmax_verf = 0.0_dp
+      dmax_verfc = 0.0_dp
+      ddmax_verfc = 0.0_dp
+      dmax_vgauss = 0.0_dp
+      ddmax_vgauss = 0.0_dp
+      dmax_gauss = 0.0_dp
+      ddmax_gauss = 0.0_dp
+      DO irep = 1, nrep
+         DO iab = 1, nab
+            !*** Coulomb: (a|1/r12|b)
+            CALL two_center_shg_integrals(operator_coulomb, vab_shg, dvab_shg, rab(:, iab), &
+                                          fba, fbb, scona_shg, sconb_shg, calculate_forces=.TRUE.)
+            CALL two_center_os_integrals(operator_coulomb, vab_os, dvab_os, rab(:, iab), fba, fbb, &
+                                         calculate_forces=.TRUE.)
+            CALL calculate_deviation(vab_shg, vab_os, dvab_shg, dvab_os, dtemp, ddtemp)
+            dmax_coulomb = MAX(dmax_coulomb, dtemp)
+            ddmax_coulomb = MAX(ddmax_coulomb, ddtemp)
+            !*** verf: (a|erf(omega*r12)/r12|b)
+            CALL two_center_shg_integrals(operator_verf, vab_shg, dvab_shg, rab(:, iab), &
+                                          fba, fbb, scona_shg, sconb_shg, omega, &
+                                          calculate_forces=.TRUE.)
+            CALL two_center_os_integrals(operator_verf, vab_os, dvab_os, rab(:, iab), fba, fbb, &
+                                         omega, calculate_forces=.TRUE.)
+            CALL calculate_deviation(vab_shg, vab_os, dvab_shg, dvab_os, dtemp, ddtemp)
+            dmax_verf = MAX(dmax_verf, dtemp)
+            ddmax_verf = MAX(ddmax_verf, ddtemp)
+            !*** verfc: (a|erfc(omega*r12)/r12|b)
+            CALL two_center_shg_integrals(operator_verfc, vab_shg, dvab_shg, rab(:, iab), &
+                                          fba, fbb, scona_shg, sconb_shg, omega, &
+                                          calculate_forces=.TRUE.)
+            CALL two_center_os_integrals(operator_verfc, vab_os, dvab_os, rab(:, iab), fba, fbb, &
+                                         omega, calculate_forces=.TRUE.)
+            CALL calculate_deviation(vab_shg, vab_os, dvab_shg, dvab_os, dtemp, ddtemp)
+            dmax_verfc = MAX(dmax_verfc, dtemp)
+            ddmax_verfc = MAX(ddmax_verfc, ddtemp)
+            !*** vgauss: (a|exp(omega*r12^2)/r12|b)
+            CALL two_center_shg_integrals(operator_vgauss, vab_shg, dvab_shg, rab(:, iab), &
+                                          fba, fbb, scona_shg, sconb_shg, omega, &
+                                          calculate_forces=.TRUE.)
+            CALL two_center_os_integrals(operator_vgauss, vab_os, dvab_os, rab(:, iab), fba, fbb, &
+                                         omega, calculate_forces=.TRUE.)
+            CALL calculate_deviation(vab_shg, vab_os, dvab_shg, dvab_os, dtemp, ddtemp)
+            dmax_vgauss = MAX(dmax_vgauss, dtemp)
+            ddmax_vgauss = MAX(ddmax_vgauss, ddtemp)
+            !*** gauss: (a|exp(omega*r12^2)|b)
+            CALL two_center_shg_integrals(operator_gauss, vab_shg, dvab_shg, rab(:, iab), &
+                                          fba, fbb, scona_shg, sconb_shg, omega, &
+                                          calculate_forces=.TRUE.)
+            CALL two_center_os_integrals(operator_gauss, vab_os, dvab_os, rab(:, iab), fba, fbb, &
+                                         omega, calculate_forces=.TRUE.)
+            CALL calculate_deviation(vab_shg, vab_os, dvab_shg, dvab_os, dtemp, ddtemp)
+            dmax_gauss = MAX(dmax_gauss, dtemp)
+            ddmax_gauss = MAX(ddmax_gauss, ddtemp)
+         ENDDO
+      ENDDO
+
+      IF (iw > 0) THEN
+         WRITE (iw, FMT="(/,T2,A)") "TEST INFO FOR 2-CENTER SHG and OS INTEGRALS:"
+         WRITE (iw, FMT="(T2,A)") "Maximal deviation between SHG and OS integrals and their derivatives"
+         WRITE (iw, FMT="(T2,A,T53,ES12.5,4X,ES12.5)") "SHG_INTEGRALS | [a|1/r12|b]", &
+            dmax_coulomb, ddmax_coulomb
+         WRITE (iw, FMT="(T2,A,T53,ES12.5,4X,ES12.5)") "SHG_INTEGRALS | [a|erf(omega*r12)/r12|b]", &
+            dmax_verf, ddmax_verf
+         WRITE (iw, FMT="(T2,A,T53,ES12.5,4X,ES12.5)") "SHG_INTEGRALS | [a|erfc(omega*r12)/r12|b]", &
+            dmax_verfc, ddmax_verfc
+         WRITE (iw, FMT="(T2,A,T53,ES12.5,4X,ES12.5)") "SHG_INTEGRALS | [a|exp(-omega*r12^2)/r12|b]", &
+            dmax_vgauss, ddmax_vgauss
+         WRITE (iw, FMT="(T2,A,T53,ES12.5,4X,ES12.5)") "SHG_INTEGRALS | [a|exp(-omega*r12^2)|b]", &
+            dmax_gauss, ddmax_gauss
+
+         IF (acc_check) THEN
+            acc_param = 1.0E-08_dp
+            IF ((dmax_coulomb >= acc_param) .OR. (ddmax_coulomb >= acc_param)) THEN
+               CPABORT("[a|1/r12|b]: Deviation between SHG and OS larger than 1.0E-08")
+            ENDIF
+            IF ((dmax_verf >= acc_param) .OR. (ddmax_verf >= acc_param)) THEN
+               CPABORT("[a|erf(omega*r12)/r12|b]: Deviation between SHG and OS larger than 1.0E-08")
+            ENDIF
+            IF ((dmax_verfc >= acc_param) .OR. (ddmax_verfc >= acc_param)) THEN
+               CPABORT("[a|erfc(omega*r12)/r12|b]: Deviation between SHG and OS larger than 1.0E-08")
+            ENDIF
+            IF ((dmax_vgauss >= acc_param) .OR. (ddmax_vgauss >= acc_param)) THEN
+               CPABORT("[a|exp(-omega*r12^2)/r12|b]: Dev. between SHG and OS larger than 1.0E-08")
+            ENDIF
+            IF ((dmax_gauss >= acc_param) .OR. (ddmax_gauss >= acc_param)) THEN
+               CPABORT("[a|exp(-omega*r12^2)|b]: Deviation between SHG and OS larger than 1.0E-08")
+            ENDIF
+         ENDIF
+      ENDIF
+      DEALLOCATE (scona_shg, sconb_shg, rab, vab_shg, dvab_shg, vab_os, dvab_os)
+      CALL deallocate_gto_basis_set(fba)
+      CALL deallocate_gto_basis_set(fbb)
+
+   END SUBROUTINE shg_integrals_perf_acc_test
+
+! **************************************************************************************************
+!> \brief contraction matrix for the SHG integrals. If these integrals are actually used,
+!>        this routine should be moved, e.g. in basis_set_types.F, where the basis set
+!>        information are initialized
+!> \param basis ...
+!> \param scon_shg contraction matrix
+! **************************************************************************************************
+   SUBROUTINE contraction_matrix_shg(basis, scon_shg)
+
+      TYPE(gto_basis_set_type), POINTER                  :: basis
+      REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: scon_shg
+
+      CHARACTER(len=*), PARAMETER :: routineN = 'contraction_matrix_shg', &
+         routineP = moduleN//':'//routineN
+
+      INTEGER                                            :: ipgf, iset, ishell, l, maxpgf, maxshell, &
+                                                            nset
+      INTEGER, DIMENSION(:), POINTER                     :: npgf, nshell
+      REAL(KIND=dp)                                      :: aif, gcc, norm_w
+      REAL(KIND=dp), DIMENSION(:, :), POINTER            :: norm, zet
+
+      nset = basis%nset
+      npgf => basis%npgf
+      nshell => basis%nshell
+      zet => basis%zet
+
+      maxpgf = SIZE(basis%gcc, 1)
+      maxshell = SIZE(basis%gcc, 2)
+      ALLOCATE (norm(basis%nset, maxshell))
+      ALLOCATE (scon_shg(maxpgf, maxshell, nset))
+      scon_shg = 0.0_dp
+
+      CALL basis_norm_shg(basis, norm)
+
+      DO iset = 1, nset
+         DO ishell = 1, nshell(iset)
+            l = basis%l(ishell, iset)
+            norm_w = norm(iset, ishell)
+            DO ipgf = 1, npgf(iset)
+               aif = 1.0_dp/((2._dp*zet(ipgf, iset))**l)
+               gcc = basis%gcc(ipgf, ishell, iset)
+               scon_shg(ipgf, ishell, iset) = norm(iset, ishell)*gcc*aif
+            END DO
+         END DO
+      END DO
+
+      DEALLOCATE (norm)
+
+   END SUBROUTINE contraction_matrix_shg
+
+!*****************************************************************************
+!> \brief normalization solid harmonic Gaussians (SHG)
+!> \param basis ...
+!> \param norm ...
+! **************************************************************************************************
+   SUBROUTINE basis_norm_shg(basis, norm)
+
+      TYPE(gto_basis_set_type), POINTER                  :: basis
+      REAL(dp), DIMENSION(:, :), POINTER                 :: norm
+
+      CHARACTER(len=*), PARAMETER :: routineN = 'basis_norm_shg', routineP = moduleN//':'//routineN
+
+      INTEGER                                            :: ipgf, iset, ishell, jpgf, l
+      REAL(KIND=dp)                                      :: aai, aaj, cci, ccj, expa, ppl
+
+      norm = 0._dp
+
+      DO iset = 1, basis%nset
+         DO ishell = 1, basis%nshell(iset)
+            l = basis%l(ishell, iset)
+            expa = 0.5_dp*REAL(2*l+3, dp)
+            ppl = fac(2*l+2)*pi**(1.5_dp)/fac(l+1)
+            ppl = ppl/(2._dp**REAL(2*l+1, dp))
+            ppl = ppl/REAL(2*l+1, dp)
+            DO ipgf = 1, basis%npgf(iset)
+               cci = basis%gcc(ipgf, ishell, iset)
+               aai = basis%zet(ipgf, iset)
+               DO jpgf = 1, basis%npgf(iset)
+                  ccj = basis%gcc(jpgf, ishell, iset)
+                  aaj = basis%zet(jpgf, iset)
+                  norm(iset, ishell) = norm(iset, ishell)+cci*ccj*ppl/(aai+aaj)**expa
+               END DO
+            END DO
+            norm(iset, ishell) = 1.0_dp/SQRT(norm(iset, ishell))
+         END DO
+      END DO
+
+   END SUBROUTINE basis_norm_shg
+
+! **************************************************************************************************
+!> \brief Calculation of the deviation between SHG and OS integrals
+!> \param vab_shg integral matrix obtained from the SHG scheme
+!> \param vab_os integral matrix obtained from the OS scheme
+!> \param dvab_shg derivative of the integrals, SHG
+!> \param dvab_os derivative of the integrals, OS
+!> \param dmax maximal deviation of vab matrices
+!> \param ddmax maximal deviation of dvab matrices
+! **************************************************************************************************
+   SUBROUTINE calculate_deviation(vab_shg, vab_os, dvab_shg, dvab_os, dmax, ddmax)
+
+      REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: vab_shg, vab_os
+      REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN)      :: dvab_shg, dvab_os
+      REAL(KIND=dp), INTENT(OUT)                         :: dmax, ddmax
+
+      CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_deviation', &
+         routineP = moduleN//':'//routineN
+
+      INTEGER                                            :: i, j, k
+      REAL(KIND=dp)                                      :: diff
+
+      dmax = 0.0_dp
+      ddmax = 0.0_dp
+
+      ! integrals vab
+      DO j = 1, SIZE(vab_shg, 2)
+         DO i = 1, SIZE(vab_shg, 1)
+            diff = ABS(vab_shg(i, j)-vab_os(i, j))
+            dmax = MAX(dmax, diff)
+         ENDDO
+      ENDDO
+
+      ! derivatives dvab
+      DO k = 1, 3
+         DO j = 1, SIZE(dvab_shg, 2)
+            DO i = 1, SIZE(dvab_shg, 1)
+               diff = ABS(dvab_shg(i, j, k)-dvab_os(i, j, k))
+               ddmax = MAX(ddmax, diff)
+            ENDDO
+         ENDDO
+      ENDDO
+
+   END SUBROUTINE
+END MODULE shg_integrals_test

src/start/input_cp2k.F

@@ -50,6 +50,7 @@ MODULE input_cp2k
    USE pw_grids,                        ONLY: do_pw_grid_blocked_false,&
                                               do_pw_grid_blocked_free,&
                                               do_pw_grid_blocked_true
+   USE shg_integrals_test,              ONLY: create_shg_integrals_test_section
    USE string_utilities,                ONLY: newline,&
                                               s2a
    USE swarm_input,                     ONLY: create_swarm_section
@@ -276,6 +277,10 @@ CONTAINS
       CALL section_add_subsection(section, subsection)
       CALL section_release(subsection)
 
+      CALL create_shg_integrals_test_section(subsection)
+      CALL section_add_subsection(section, subsection)
+      CALL section_release(subsection)
+
    END SUBROUTINE create_test_section
 
 ! **************************************************************************************************

tests/LIBTEST/TEST_FILES

@@ -29,4 +29,5 @@ dbcsr_order_N.inp                                     58    1.0E-14
 test_eri_mme_accuracy.inp                              0
 test_eri_mme_performance.inp                           0
 test_minimax.inp                                       0
+test_shg_integrals_01.inp                              0
 #EOF

tests/LIBTEST/test_shg_integrals_01.inp

@@ -0,0 +1,70 @@
+&GLOBAL
+  PROJECT test_shg_integrals
+  PRINT_LEVEL MEDIUM
+  PROGRAM_NAME TEST
+  RUN_TYPE NONE
+&END GLOBAL
+&TEST
+  &SHG_INTEGRALS_TEST
+    ABC 10.0 10.0 10.0
+    NAB_MIN 8
+    NREP 2
+    CHECK_ACCURACY 
+    &BASIS 
+    10
+    2   0   0   1  1 
+      22.956000679816   1.0   
+    2   0   1   1  1  1  
+      11.437045132575   1.0   1.0 
+    2   0   2   1  1  1  1
+       5.698118029748   1.0   1.0   1.0 
+    2   0   2   1  1  1  1                 
+       2.838893149810   1.0   1.0   1.0
+    2   0   3   1  1  1  1  1               
+       1.414381779030   1.0   1.0   1.0  1.0
+    2   0   3   1  1  1  1  1               
+       0.704667527549   1.0   1.0   1.0  1.0
+    2   0   3   1  1  1  1  1               
+       0.351076584656   1.0   1.0   1.0  1.0
+    2   0   3   1  1  1  1  1               
+       0.174911945670   1.0   1.0   1.0  1.0
+    2   0   3   1  1  1  1  1               
+       0.087143916955   1.0   1.0   1.0  1.0
+    2   0   3   1  1  1  1  1               
+       0.043416487268   1.0   1.0   1.0  1.0
+    &END 
+   &BASIS  
+  15
+   2   0   0   1  1  
+     24.031909411024   1.0
+   2   0   0   1  1  
+     16.167926705922   1.0   
+   2   0   1   1  1  1  
+     10.877281929506   1.0   1.0   
+   2   0   1   1  1  1  
+      7.317899463920   1.0   1.0   
+   2   0   2   1  1  1  1 
+      4.923256831173   1.0   1.0   1.0   
+   2   0   3   1  1  1  1  1  
+      3.312215198528   1.0   1.0   1.0   1.0   
+   2   0   3   1  1  1  1  1  
+      2.228356126354   1.0   1.0   1.0   1.0   
+   2   0   4   1  1  1  1  1  1
+      1.499169204968   1.0   1.0   1.0   1.0   1.0
+   2   0   4   1  1  1  1  1  1
+      1.008594756710   1.0   1.0   1.0   1.0   1.0
+   2   0   4   1  1  1  1  1  1
+      0.678551413605   1.0   1.0   1.0   1.0   1.0
+   2   0   4   1  1  1  1  1  1
+      0.456508441910   1.0   1.0   1.0   1.0   1.0
+   2   0   4   1  1  1  1  1  1
+      0.307124785767   1.0   1.0   1.0   1.0   1.0
+   2   0   4   1  1  1  1  1  1
+      0.206624073890   1.0   1.0   1.0   1.0   1.0
+   2   0   4   1  1  1  1  1  1
+      0.139010297734   1.0   1.0   1.0   1.0   1.0
+   2   0   4   1  1  1  1  1  1 
+    0.093521836600   1.0   1.0   1.0   1.0   1.0
+    &END
+  &END SHG_INTEGRALS_TEST
+&END TEST