calc_pgen_mol_guga_single_orbs Function

    function calc_pgen_mol_guga_single_orbs(ilut, nI, csf_i, excitInfo) result(pgen)
        integer(n_int), intent(in) :: ilut(0:niftot)
        integer, intent(in) :: nI(nel)
        type(CSF_Info_t), intent(in) :: csf_i
        type(ExcitationInformation_t), intent(in) :: excitInfo
        real(dp) :: pgen
        real(dp) :: p_elec, p_orb, cum_sum
        integer :: cc_i, nOrb, ierr, i, a
        real(dp), allocatable :: cum_arr(:)

        call get_orbs_from_excit_info(excitInfo, i, a)

        if (IsDoub(ilut, i)) then
            p_elec = 2.0 / real(nel, dp)
        else if (IsNotOcc(ilut, i)) then
            pgen = 0.0_dp
            return
        else
            p_elec = 1.0 / real(nel, dp)
        end if

        if (IsDoub(ilut, a)) then
            pgen = 0.0_dp
            return
        end if

        if (gtID(i) == gtID(a)) then
            pgen = 0.0_dp
            return
        end if

        cc_i = ClassCountInd(1, SpinOrbSymLabel(a), G1(a)%Ml)

        nOrb = OrbClassCount(cc_i)
        allocate(cum_arr(nOrb), stat=ierr)

        if (IsDoub(ilut, i)) then
            call gen_cum_list_guga_single_3(nI, csf_i, i, cc_i, cum_arr)
        else
            if (is_beta(i)) then
                call gen_cum_list_guga_single_1(nI, csf_i, i, cc_i, cum_arr)
            else
                call gen_cum_list_guga_single_2(nI, csf_i, i, cc_i, cum_arr)
            end if
        end if

        if (near_zero(cum_sum)) then
            pgen = 0.0_dp
            return
        end if

        if (a == 1) then
            p_orb = cum_arr(1) / cum_sum
        else
            p_orb = (cum_arr(a) - cum_arr(a - 1)) / cum_sum
        end if

        pgen = p_orb * p_elec

    end function calc_pgen_mol_guga_single_orbs