calc_pgen_rs_hubbard_transcorr Function

    function calc_pgen_rs_hubbard_transcorr(nI, ilutI, ex, ic) result(pgen)
        integer, intent(in) :: nI(nel), ex(2, 2), ic
        integer(n_int), intent(in) :: ilutI(0:NIfTot)
        real(dp) :: pgen
#ifdef DEBUG_
        character(*), parameter :: this_routine = "calc_pgen_rs_hubbard_transcorr"
#endif
        integer :: src(2), tgt(2)
        real(dp) :: p_elec, p_orb, cum_arr(nBasis / 2), cum_sum, p_hole_1, &
                    p_orb_a, p_orb_b

        src = ex(1, :)
        tgt = ex(2, :)

        if (ic == 1) then

            ASSERT(same_spin(src(1), tgt(1)))

            p_elec = 1.0_dp / real(nel, dp)

            call create_cum_list_rs_hubbard_transcorr_single(nI, ilutI, src(1), &
                                                             cum_arr, cum_sum, tgt(1), p_orb)

            if (cum_sum < EPS) then
                pgen = 0.0_dp
                return
            end if

            pgen = p_elec * p_orb * (1.0_dp - pDoubles)

        else if (ic == 2) then

            if (same_spin(src(1), src(2)) .or. same_spin(tgt(1), tgt(2))) then
                pgen = 0.0_dp
                return
            end if

            p_elec = 1.0_dp / real(nOccAlpha * nOccBeta, dp)

            ! we need two holes..
            ! pick the first at random..
            p_hole_1 = 1.0_dp / real(nBasis - nel, dp)

            call create_cum_list_rs_hubbard_transcorr_double(nI, ilutI, src, tgt(1), &
                                                             cum_arr, cum_sum, tgt(2), p_orb_a)

            if (cum_sum < EPS) then
                pgen = 0.0_dp
                return
            end if

            ! and the other way around
            call create_cum_list_rs_hubbard_transcorr_double(nI, ilutI, src, tgt(2), &
                                                             cum_arr, cum_sum, tgt(1), p_orb_b)

            if (cum_sum < EPS) then
                pgen = 0.0_dp
                return
            end if

            pgen = 1.0_dp * p_elec * p_hole_1 * (p_orb_a + p_orb_b) * pDoubles

        else
            ! every other ic is 0
            pgen = 0.0_dp

        end if

    end function calc_pgen_rs_hubbard_transcorr