get_rand_orb Subroutine

    subroutine get_rand_orb(nI, tgt, ms, nPicked, pgen)
        integer, intent(inout) :: tgt(3)
        integer, intent(in) :: ms, nPicked, nI(nel)
        real(dp), intent(inout) :: pgen

        integer :: pool ! available orbitals
        integer :: iOrb, i
        real(dp) :: r

        ! get the number of possible orbitals
        if (ms > 0) then
            pool = nUnoccAlpha
        else
            pool = nUnoccBeta
        end if
        ! we need to see how many same spin orbs have been picked so far
        do i = 1, nPicked
            if ((ms > 0) .neqv. is_beta(tgt(i))) pool = pool - 1
        end do

        ! pick a random index
        r = genrand_real2_dSFMT()
        iOrb = spinOrb(int(r * pool) + 1)

        ! check if the orb is already targeted
        call skipPicked()
        ! assign the orbital
        tgt(nPicked + 1) = iOrb

        ! adjust the probability
        pgen = pgen / pool

        ! we need to sort tgt (see above)
        tgt(1:(nPicked + 1)) = sort_unique(tgt(1:(nPicked + 1)))

    contains
        pure function spinOrb(orb) result(sorb)
            integer, intent(in) :: orb
            integer :: sorb

            sorb = 2 * orb + (ms - 1) / 2
        end function spinOrb

        subroutine skipPicked()
            integer :: i
            integer :: invalidOrbs(nel + nPicked)

            invalidOrbs(1:nel) = nI
            invalidOrbs((nel + 1):(nel + nPicked)) = tgt(1:nPicked)

            invalidOrbs = sort_unique(invalidOrbs)

            do i = 1, (nPicked + nel)
                ! check if the orb is already targeted
                ! assumes tgt is sorted
                if (invalidOrbs(i) <= iOrb .and. G1(invalidOrbs(i))%ms == ms) then
                    iOrb = iOrb + 2
                end if
            end do
        end subroutine skipPicked
    end subroutine get_rand_orb