pickRandomOrb_vector Subroutine

    subroutine pickRandomOrb_vector(csf_i, orbRes, pgen, orb, occRes)
        ! this picks a random orb under multiple orbital or occupation
        ! number restrictions
        type(CSF_Info_t), intent(in) :: csf_i
        integer, intent(in) :: orbRes(:)
        real(dp), intent(inout) :: pgen
        integer, intent(out) :: orb
        integer, intent(in), optional :: occRes
        character(*), parameter :: this_routine = "pickRandomOrb_vector"

        integer :: r, nOrbs, ierr, num, i
        logical :: mask(nSpatOrbs)
        integer, allocatable :: resOrbs(:)

        ASSERT(all(orbRes >= 0 .and. orbRes <= nSpatOrbs))

        num = size(orbRes)

        ! do not need to check if orbRes is empty, since i call it explicetly
        ! for multiple orbitals
        mask = .true.

        if (present(occRes)) then
            ASSERT(occRes >= 0 .and. occRes <= 2)
            mask = (csf_i%Occ_int /= occRes)
        end if

        do i = 1, num
            mask = (mask .and. orbitalIndex /= orbRes(i))
        end do

        nOrbs = count(mask)

        if (nOrbs > 0) then
            allocate(resOrbs(nOrbs), stat=ierr)

            resOrbs = pack(orbitalIndex, mask)

            r = 1 + floor(genrand_real2_dSFMT() * real(nOrbs, dp))

            orb = resOrbs(r)

            pgen = pgen / real(nOrbs, dp)
            deallocate(resOrbs)
        else
            orb = 0
            pgen = 0.0_dp
        end if

    end subroutine pickRandomOrb_vector