subroutine ints_exchange(basis, schwarz_ints, mu2)
use int2e_rotaxis, only: genr22
use int2e_libint, only: libint2_init_eri, libint2_cleanup_eri
use int2e_libint, only: libint_compute_eri, libint_print_eri
use int2e_libint, only: libint_t, libint2_active
use int2e_rys, only: int2_rys_compute
use types, only: information
use constants, only: NUM_CART_BF
use, intrinsic :: iso_c_binding, only: C_NULL_PTR, C_INT, c_f_pointer
implicit none
type(basis_set), intent(in) :: basis
real(kind=dp), intent(inout) :: schwarz_ints(:,:)
real(kind=dp), optional, intent(in) :: mu2
real(kind=dp), parameter :: &
ic_exchng = 1.0d-15, &
ei1_exchng = 1.0d-17, &
ei2_exchng = 1.0d-17, &
cux_exchng = 50.0
integer :: flips(4)
integer :: shell_ids(4)
integer :: nbf(4)
integer :: lmax
integer :: ish, jsh
integer :: i, j, jmax
integer :: am(4), max_am
integer :: ok
logical :: rotspd, libint, zero_shq, rys
logical :: attenuated
real(kind=dp) :: vmax
real(kind=dp), allocatable, target :: ints(:)
real(kind=dp), pointer :: pints(:,:,:,:)
type(libint_t), allocatable :: erieval(:)
type(int2_rys_data_t) :: gdat
type(int2_cutoffs_t) :: cutoffs
type(int2_pair_storage) :: ppairs
attenuated = present(mu2)
lmax = maxval(basis%am)
if (lmax < 0 .or. lmax > 6) call show_message("Basis set agular momentum exceeds max. supported", WITH_ABORT)
! Set very tight cutoff
call cutoffs%set(ic_exchng, ei1_exchng, ei2_exchng, cux_exchng)
if (libint2_active) then
allocate(erieval(basis%mxcontr**4))
call libint2_init_eri(erieval, int(4, C_INT), C_NULL_PTR)
end if
allocate(ints(NUM_CART_BF(lmax)**4), source=0.0d0)
call gdat%init(lmax, cutoffs, ok)
call ppairs%alloc(basis, cutoffs)
call ppairs%compute(basis, cutoffs)
do ish = 1, basis%nshell
do jsh = 1, ish
shell_ids = [ish, jsh, ish, jsh]
am = basis%am(shell_ids)
max_am = maxval(am)
rotspd = max_am <= 2
libint = .not.rotspd.and.libint2_active.and..not.attenuated
rys = .not.rotspd.and..not.libint
if (rotspd) then
if (attenuated) then
call genr22(basis, ppairs, ints, shell_ids, flips, cutoffs, mu2)
else
call genr22(basis, ppairs, ints, shell_ids, flips, cutoffs)
end if
nbf = am(flips)
nbf = (nbf+1)*(nbf+2)/2
vmax = maxval(abs(ints(1:product(nbf))))
else if (libint) then
nbf = am(flips)
call libint_compute_eri(basis, ppairs, cutoffs, shell_ids, 0, erieval, flips, zero_shq)
call c_f_pointer(erieval(1)%targets(1), pints, shape=nbf([4,3,2,1]))
vmax = maxval(abs(pints))
else if (rys) then
call gdat%set_ids(basis, shell_ids)
if (attenuated) then
call int2_rys_compute(ints, gdat, ppairs, zero_shq, mu2)
else
call int2_rys_compute(ints, gdat, ppairs, zero_shq)
end if
nbf = gdat%nbf
pints(1:nbf(4), 1:nbf(3), 1:nbf(2), 1:nbf(1)) => ints
call normalize_ints(nbf, gdat%am, pints)
vmax = maxval(abs(pints))
end if
schwarz_ints(ish, jsh) = sqrt(vmax)
schwarz_ints(jsh, ish) = sqrt(vmax)
end do
end do
if (libint2_active) then
call libint2_cleanup_eri(erieval)
deallocate(erieval)
end if
call gdat%clean()
end subroutine ints_exchange