pickRandomOrb_scalar Subroutine

    subroutine pickRandomOrb_scalar(csf_i, orbRes, pgen, orb, occRes)
        ! routine to pick a random orbital under certain orbital and/or
        ! occupation number restrictions and gives the probability to pick
        ! this orbital
        ! these orbitals for now are chosen uniformly and not with any
        ! cauchy schwarz like criteria involving the one- and two-particle
        ! integrals
        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_scalar"

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

        ASSERT(orbRes >= 0 .and. orbRes <= nSpatOrbs)
        ! could make this a global, permanent variable... should even!
!         orbInd = [ (r, r = 1, nBasis/2) ]

        ! create a list of the available orbitals under the provided
        ! restrictions:
        if (orbRes == 0) then
            ! in this case a occRes should be provided or else one could
            ! just create a random number between 1, nBasis/2 without this
            ! function
            ASSERT(present(occRes))
            ASSERT(occRes >= 0 .and. occRes <= 2)

            mask = (csf_i%Occ_int /= occRes)

        else
            ! check i occRes is present
            if (present(occRes)) then
                ASSERT(occRes >= 0 .and. occRes <= 2)
                mask = ((csf_i%Occ_int /= occRes) .and. (orbitalIndex /= orbRes))

            else
                ! only orbital restriction
                mask = (orbitalIndex /= orbRes)

            end if
        end if

        nOrbs = count(mask)
        ! here i could use nOrbs to indicate that no excitation is possible
        ! no -> just check outside if resorbs is size 0
        ! but i get a error further down if i access resOrbs(1)

        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)

            ! TODO: modify the pgen probabilities!
            ! here it should be 1/nOrbs, since its uniform isnt it?
            pgen = pgen / real(nOrbs, dp)

            deallocate(resOrbs)
        else
            ! indicate that no orbital is fitting with r = 0
            orb = 0

            ! do i have to set pgen to zero? probably
            pgen = 0.0_dp
        end if

    end subroutine pickRandomOrb_scalar