int2_rys_compute Subroutine

  subroutine int2_rys_compute(ints, gdat, ppairs, zero_shq, mu2)

    use int2_pairs, only: int2_pair_storage
    implicit none

    type(int2_rys_data_t) :: gdat
    type(int2_pair_storage), intent(in) :: ppairs
    real(kind=dp), intent(inout) :: ints(*)
    real(kind=dp), intent(in), optional :: mu2
    logical :: zero_shq

    integer :: ijg, klg, maxgg, mmax, ng
    integer :: nmax
    real(kind=dp) :: aa, ab, aandb1, bb, da, db, test
    real(kind=dp) :: pfac, rho
    real(kind=dp) :: p(3), q(3)
    logical :: first
    integer :: id1, id2, ppid_p, ppid_q, npp_p, npp_q
    real(kind=dp) :: mu2_1

    ! Range-separation parameter for Erfc-attenuated integrals
    mu2_1 = 0
    if (present(mu2)) mu2_1 = 1.0d0/mu2

!   Prepare shell block
    call set_shells(gdat)


    id1 = maxval(gdat%id(1:2))
    id2 = minval(gdat%id(1:2))
    npp_p = ppairs%ppid(1,id1*(id1-1)/2+id2)
    ppid_p = ppairs%ppid(2,id1*(id1-1)/2+id2)

    id1 = maxval(gdat%id(3:4))
    id2 = minval(gdat%id(3:4))
    npp_q = ppairs%ppid(1,id1*(id1-1)/2+id2)
    ppid_q = ppairs%ppid(2,id1*(id1-1)/2+id2)

    zero_shq = npp_p*npp_q == 0
    if (zero_shq) return

    nmax = gdat%am(1)+gdat%am(2)+1
    mmax = gdat%am(3)+gdat%am(4)+1

    maxgg = MAXCONTR/gdat%nroots

!   Pair of k,l primitives

    first = .true.
    zero_shq = .true.
    ng = 0

    do klg = 1, npp_q
      db = ppairs%k(ppid_q-1+klg)*ppairs%ginv(ppid_q-1+klg)
      bb = ppairs%g(ppid_q-1+klg)
      q = ppairs%P(:,ppid_q-1+klg)

!     Pair of i,j primitives
      do ijg = 1, npp_p
        da = ppairs%k(ppid_p-1+ijg)*ppairs%ginv(ppid_p-1+ijg)
        aa = ppairs%g(ppid_p-1+ijg)
        p = ppairs%P(:,ppid_p-1+ijg)

        ! 2nd term is used for Erfc-attenuated integrals
        ! It is zero for regular integrals
        ab = (aa+bb) + aa*bb*mu2_1

        pfac = da*db
        test = pfac*pfac

        if (test<gdat%quartet_cutoff*ab) cycle

        ng = ng+1

        aandb1 = 1.0_dp/ab
        rho = aa*bb*aandb1

        gdat%abv(1,ng) = ppairs%ginv(ppid_p-1+ijg)
        gdat%abv(2,ng) = ppairs%ginv(ppid_q-1+klg)
        gdat%abv(3,ng) = rho
        gdat%abv(4,ng) = pfac*sqrt(aandb1)
        gdat%abv(5,ng) = aandb1
        gdat%abv(6,ng) = rho*sum((p-q)**2)

        gdat%PQ(:,ng) = p-q

        if (nmax>1) gdat%pb(:,ng) = ppairs%PB(:,ppid_p-1+ijg)
        if (mmax>1) gdat%qd(:,ng) = ppairs%PB(:,ppid_q-1+klg)

        gdat%dij(:,ng) = ppairs%PA(:,ppid_p-1+ijg)-ppairs%PB(:,ppid_p-1+ijg)
        gdat%dkl(:,ng) = ppairs%PA(:,ppid_q-1+klg)-ppairs%PB(:,ppid_q-1+klg)

        if (ng==maxgg) then
          if (ng/=0) then
            if (first) call clear_ints(gdat, ints)
            first = .false.
            call compute(gdat, ng, nmax, mmax, ints)
            ng = 0
            zero_shq = .false.
          end if
        end if

      end do
    end do
    if (ng/=0) then
      if (first) call clear_ints(gdat, ints)
      first = .false.
      call compute(gdat, ng, nmax, mmax, ints)
      ng = 0
      zero_shq = .false.
    end if


  end subroutine int2_rys_compute