pick_occupied_orbital_single Subroutine

    subroutine pick_occupied_orbital_single(ilut, src, cc_index, part_type, &
                                            pgen, orb, calc_pgen)
        integer, intent(in) :: src, cc_index, part_type
        integer(n_int), intent(in) :: ilut(0:niftot)
        real(dp), intent(out) :: pgen
        integer, intent(out) :: orb
        logical, intent(in), optional :: calc_pgen

        ! routine to pick an orbital from the occupied manifold in the
        ! reference determinant for single excitations
        ! i have to take symmetry into account now..  that complicates
        ! things.. and spin..
        character(*), parameter :: this_routine = "pick_occupied_orbital_single"

        integer :: n_valid, j, occ_orbs(nel), i, ind, norb, label_index

        j = 1
        occ_orbs = -1
        orb = -1

        norb = OrbClassCount(cc_index)
        label_index = SymLabelCounts2(1, cc_index)

        ! damn i did not include symmetries todo
        ! ok do it now with symmetries
        do i = 1, norb
            orb = SymLabelList2(label_index + i - 1)
            if (is_in_ref(orb, part_type) .and. IsNotOcc(ilut, orb)) then

                ASSERT(SpinOrbSymLabel(orb) == SpinOrbSymLabel(src))

                occ_orbs(j) = orb
                j = j + 1
            end if
        end do

        n_valid = j - 1

        if (n_valid == 0) then
            orb = 0
            pgen = 0.0_dp
            return
        end if

        ASSERT(n_valid > 0)
        pgen = 1.0_dp / real(n_valid, dp)

        if (present(calc_pgen)) then
            if (calc_pgen) return
        end if

        ! else pick uniformly from that available list..
        ind = 1 + int(genrand_real2_dSFMT() * n_valid)
        orb = occ_orbs(ind)

        ASSERT(orb > 0)

    end subroutine pick_occupied_orbital_single