grd2_rys_compute Subroutine

  subroutine grd2_rys_compute(gdat, ppairs, dab, dabmax, mu2)

    use int2_pairs, only: int2_pair_storage, int2_cutoffs_t
    implicit none

    type(grd2_int_data_t) :: gdat
    type(int2_pair_storage), intent(in) :: ppairs
    real(kind=dp), intent(in) :: dabmax
    real(kind=dp), intent(in) :: dab(*)
    real(kind=dp), intent(in), optional :: mu2

    integer :: ijg, klg, maxgg, mmax, ng
    integer :: nimax, njmax, nkmax, nlmax, nmax
    real(kind=dp) :: aa, ab, aandb1, bb, da, db, test
    real(kind=dp) :: pfac, rho
    real(kind=dp) :: p(3), q(3)
    logical :: last
    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)
    if (npp_p*npp_q == 0) return

    nimax = gdat%am(1) + gdat%der(1) + 1
    njmax = gdat%am(2) + gdat%der(2) + 1
    nkmax = gdat%am(3) + gdat%der(3) + 1
    nlmax = gdat%am(4) + gdat%der(4) + 1

    nmax = gdat%am(1)+gdat%am(2)+1 + min(gdat%der(1)+gdat%der(2),gdat%nder)
    mmax = gdat%am(3)+gdat%am(4)+1 + min(gdat%der(3)+gdat%der(4),gdat%nder)

    maxgg = MAXCONTR/gdat%nroots

!   Pair of k,l primitives

    gdat%fd = 0
    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%dtol*ab) cycle
        if (test*dabmax*dabmax<gdat%dabcut*ab) cycle

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

        ng = ng+1

        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%ai(ng) = 2*ppairs%alpha_a(ppid_p-1+ijg)
        gdat%aj(ng) = 2*ppairs%alpha_b(ppid_p-1+ijg)
        gdat%ak(ng) = 2*ppairs%alpha_a(ppid_q-1+klg)
        gdat%al(ng) = 2*ppairs%alpha_b(ppid_q-1+klg)

        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)

        last = klg==npp_q .and. ijg==npp_p

        if (ng==maxgg .or. last) then
          if (ng==0) return

          call compute_grd_ints(gdat, dab, &
                  ng, nmax, mmax, nimax, njmax, nkmax, nlmax)

          ng = 0
        end if

      end do
    end do

!   Process derivative integrals
    call apply_translation_invariance(gdat)

  end subroutine grd2_rys_compute