pick_a_orb Function

    function pick_a_orb(ilut, src, ispn, cpt, cum_sum, cum_arr) result(orb)

        integer(n_int), intent(in) :: ilut(0:NifTot)
        integer, intent(in) :: src(2), iSpn
        real(dp), intent(out) :: cpt, cum_sum
        real(dp), intent(inout) :: cum_arr(nbasis)
        integer :: orb
        character(*), parameter :: t_r = 'pick_a_orb'
        character(*), parameter :: this_routine = t_r

        real(dp) :: r, cum_tst, cpt_tst
        integer :: start_ind, srcid(2)

        logical :: beta, parallel, valid

        ! Just in case. eeep.
        ! Note that we scale spin/spatial orbitals here with factors of 2
        ! which need more care if using spin orbitals
        if (tUHF .or. tStoreSpinOrbs) &
            call stop_all(this_routine, "ASSUMES RHF orbitals")

        ! Generate the cumulative list for making the selection, and bail if
        ! there is no selection avaialable
        call gen_a_orb_cum_list(ilut, src, ispn, cum_arr)
        cum_sum = cum_arr(nbasis)
        if (cum_sum < EPS) then
            orb = 0
            return
        end if

        ! Pick the orbital, and extract the relevant list components for
        ! probability generation purposes
        r = genrand_real2_dSFMT() * cum_sum
        orb = binary_search_first_ge(cum_arr, r)
        if (orb == 1) then
            cpt = cum_arr(1)
        else
            cpt = cum_arr(orb) - cum_arr(orb - 1)
        end if

#ifdef DEBUG_
        call pgen_select_a_orb(ilut, src, orb, iSpn, cpt_tst, cum_tst, &
                               cum_arr, .false.)
        if (abs(cpt_tst - cpt) > 1e-6 .or. abs(cum_tst - cum_sum) > 1e-6) then
            call stop_all(t_r, 'Calculated probability does not match')
        end if
#endif

    end function pick_a_orb