subroutine pick_spin_opp_holes(ilutI, orbs, p_orb)
! routine to pick two spin-opposite unoccupied orbitals
integer(n_int), intent(in) :: ilutI(0:NIfTot)
integer, intent(out) :: orbs(2)
real(dp), intent(out) :: p_orb
integer, parameter :: max_trials = 100
integer :: cnt
cnt = 0
orbs = 0
p_orb = 0.0_dp
do while (cnt < max_trials)
orbs(1) = 1 + int(genrand_real2_dsfmt() * nBasis)
if (IsOcc(ilutI, orbs(1))) cycle
do
orbs(2) = 1 + int(genrand_real2_dsfmt() * nBasis)
if (IsOcc(ilutI, orbs(2))) cycle
if (orbs(1) /= orbs(2)) exit
end do
if (.not. same_spin(orbs(1), orbs(2))) exit
end do
orbs = [minval(orbs), maxval(orbs)]
p_orb = 1.0_dp / real((nBasis / 2 - nOccAlpha) * (nBasis / 2 - nOccBeta), dp)
end subroutine pick_spin_opp_holes