calc_pgen_rs_hubbard Function

    function calc_pgen_rs_hubbard(ilutI, ex, ic) result(pgen)
        ! i also need a pgen recalculator.. specifically for the HPHF
        ! implementation and i need to take the transcorrelated keyword
        ! into account here!
        integer, intent(in) :: 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"
#endif
        integer :: src, tgt
        real(dp) :: p_elec, p_orb, cum_sum
        real(dp), allocatable :: cum_arr(:)

        ! only single excitations in the real-space hubbard
        if (ic /= 1) then
            pgen = 0.0_dp
            return
        end if

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

        ! can i assert the same spin of the 2 involved orbitals?
        ! just return 0 if both have different spin?
        ASSERT(same_spin(src, tgt))

        ! and assert that we actually take a valid excitation:
        ASSERT(any(tgt == lat%get_spinorb_neighbors(src)))
        ASSERT(IsOcc(ilutI, src))
        ASSERT(IsNotOcc(ilutI, tgt))

        p_elec = 1.0_dp / real(nel, dp)

        call create_cum_list_rs_hubbard(ilutI, src, lat%get_spinorb_neighbors(src), &
                                        cum_arr, cum_sum, tgt, p_orb)
        if (cum_sum < EPS) then
            pgen = 0.0_dp
            return
        end if

        pgen = p_elec * p_orb

    end function calc_pgen_rs_hubbard