pick_random_hole Subroutine

    subroutine pick_random_hole(ilutI, orb, p_orb, spin)
        ! routine to pick an unoccupied spin-orbital (orb) with probability
        ! p_orb
        ! with an additional spin input we can restrict ourself to a
        ! specific parallel spin!
        ! spin = 1 -> beta orbital to be picked
        ! spin = 0 -> alpha orbital to be picked
        integer(n_int), intent(in) :: ilutI(0:NIfTot)
        integer, intent(out) :: orb
        real(dp), intent(out) :: p_orb
        integer, intent(in), optional :: spin
#ifdef DEBUG_
        character(*), parameter :: this_routine = "pick_random_hole"
#endif
        integer, parameter :: max_trials = 100
        integer :: cnt, i

        orb = 0
        cnt = 0
        p_orb = 0.0_dp

        if (present(spin)) then
            ASSERT(spin == 0 .or. spin == 1)

            do while (cnt < max_trials)
                ! create a random spin-orbital of parallel spin
                i = 2 * (1 + int(genrand_real2_dsfmt() * nBasis / 2)) - spin

                if (IsNotOcc(ilutI, i)) then
                    orb = i
                    if (spin == 0) then
                        p_orb = 1.0_dp / real(nBasis / 2 - nOccAlpha, dp)
                    else if (spin == 1) then
                        p_orb = 1.0_dp / real(nBasis / 2 - nOccBeta, dp)
                    end if
                    return
                end if
            end do

        else

            do while (cnt < max_trials)
                cnt = cnt + 1
                i = 1 + int(genrand_real2_dsfmt() * nBasis)

                if (IsNotOcc(ilutI, i)) then
                    orb = i
                    p_orb = 1.0_dp / real(nBasis - nel, dp)
                    return
                end if
            end do
        end if

    end subroutine pick_random_hole