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