gen_double_excit_rs_hub_transcorr Subroutine

    subroutine gen_double_excit_rs_hub_transcorr(nI, ilutI, nJ, ilutJ, ex, tPar, pgen)
        integer, intent(in) :: nI(nel)
        integer(n_int), intent(in) :: ilutI(0:NIfTot)
        integer, intent(out) :: nJ(nel), ex(2, 2)
        integer(n_int), intent(out) :: ilutJ(0:NIfTot)
        logical, intent(out) :: tPar
        real(dp), intent(out) :: pgen
#ifdef DEBUG_
        character(*), parameter :: this_routine = "gen_double_excit_rs_hub_transcorr"
#endif
        integer :: elecs(2), orbs(2), src(2), ind
        real(dp) :: p_elec, cum_arr(nBasis / 2), cum_sum, p_orb_a, p_orb_b, p_orb_switch

        call pick_spin_opp_elecs(nI, elecs, p_elec)

        src = nI(elecs)

        ! pick the first hole at random
        call pick_random_hole(ilutI, orbs(1), p_orb_a)

        if (orbs(1) == 0) then
            nJ(1) = 0
            pgen = 0.0_dp
            return
        end if

        ! create the cum-list for b
        call create_cum_list_rs_hubbard_transcorr_double(nI, ilutI, src, orbs(1), cum_arr, &
                                                         cum_sum)

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

        call pick_from_cum_list(cum_arr, cum_sum, ind, p_orb_b)

        ! pick the right spin
        if (is_beta(orbs(1))) then
            orbs(2) = 2 * ind
        else
            orbs(2) = 2 * ind - 1
        end if

        ! now pick the other way around
        call create_cum_list_rs_hubbard_transcorr_double(nI, ilutI, src, orbs(2), cum_arr, &
                                                         cum_sum, orbs(1), p_orb_switch)

        ! if cum_sum can be 0 here i made something wrong with the cum_sum
        ! check above!
        ASSERT(cum_sum > EPS)

        ! no factor of 2, since we add the a|b b|a already!
        pgen = 1.0_dp * p_elec * p_orb_a * (p_orb_b + p_orb_switch)

        call make_double(nI, nJ, elecs(1), elecs(2), orbs(1), orbs(2), ex, tPar)
        ilutJ = make_ilutJ(ilutI, ex, 2)

    end subroutine gen_double_excit_rs_hub_transcorr