calc_pgen_k_space_hubbard Function

    function calc_pgen_k_space_hubbard(nI, ilutI, ex, ic) result(pgen)
        integer(n_int), intent(in) :: ilutI(0:niftot)
        integer, intent(in) :: nI(nel), ex(2, 2), ic
        real(dp) :: pgen
#ifdef DEBUG_
        real(dp) :: test
#endif
        real(dp) :: p_elec, p_orb, cum_arr(nbasis), cum_sum
        integer :: orb_list(nbasis, 2), src(2)

        if (ic /= 2) then
            pgen = 0.0_dp
            return
        end if

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

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

        src = get_src(ex)

        call create_ab_list_hubbard(nI, ilutI, src, orb_list, cum_arr, cum_sum, &
                                    ex(2, 1), p_orb)

#ifdef DEBUG_
        call create_ab_list_hubbard(nI, ilutI, src, orb_list, cum_arr, cum_sum, &
                                    ex(2, 2), test)

        if (abs(test - p_orb) > 1.e-8) then
            print *, "pgen assumption wrong:!"
            print *, "p_orb: ", p_orb
            print *, "test: ", test
            print *, "ex(2,:): ", ex(2, :)
        end if

#endif

        ! i do not need to recalc, the p(b|ij) since it is the same..
        ! but i need a factor of 2 somewhere.. figure that out!
        pgen = p_elec * p_orb * 2.0_dp

    end function calc_pgen_k_space_hubbard