pick_a_orb_guga_mol Subroutine

    subroutine pick_a_orb_guga_mol(csf_i, occ_orbs, contrib, cum_sum, cum_arr, orb_a)
        ! general routine, which picks orbital a for a  double excitation in
        ! the guga formalism, with symmetry restrictions and weighted
        ! with the FCIDUMP integrals. This is for MOLECULAR calculations,
        ! since in Hubbard and UEG type calculation with existing k-point
        ! restrictions, there is a more efficient and direct way to pick
        ! weighted with the actual matrix elemetn
        type(CSF_Info_t), intent(in) :: csf_i
        integer, intent(in) :: occ_orbs(2)
        real(dp), intent(out) :: contrib, cum_sum, cum_arr(nSpatOrbs)
        integer, intent(out) :: orb_a

        real(dp) :: r
        ! there are no additional GUGA restrictions on the orbital a yet, only
        ! that it has to be a non-occupied "spin"-orbital
        ! so do it in the same way as already implemented in NECI in
        ! symrandexcit5!

        ! generate the cummulative pgen list:
        if (tGen_guga_weighted) then
            call gen_a_orb_cum_list_guga_mol(csf_i, occ_orbs, cum_arr)
        else
            cum_arr = csf_i%cum_list
        end if

        cum_sum = cum_arr(nSpatOrbs)
        ! check if no excitation is possible
        if (near_zero(cum_sum)) then
            orb_a = 0
            return
        end if

        ! then pick the orbitals according to the list
        r = genrand_real2_dSFMT() * cum_sum
        orb_a = binary_search_first_ge(cum_arr, r)

        if (orb_a == 1) then
            contrib = cum_arr(1)
        else
            contrib = cum_arr(orb_a) - cum_arr(orb_a - 1)
        end if

        ! maybe similar to simon to a DEBUG check if the pgens are correct
        ! TODO

    end subroutine pick_a_orb_guga_mol