subroutine gen_double_ueg(nI, ilutI, nJ, ilutJ, tPar, ex, pgen)
integer, intent(in) :: nI(nel)
integer(n_int), intent(in) :: ilutI(0:niftot)
integer, intent(out) :: nJ(nel), ex(2, maxExcit)
integer(n_int), intent(out) :: ilutJ(0:niftot)
logical, intent(out) :: tPar
real(dp), intent(out) :: pgen
integer :: elecs(2), eleci, elecj, ispn, orba, orbb, orbi, orbj
real(dp) :: pelec, cum_sum, elem, porb, r, cum_arr(nBasis)
! Pick a pair of electrons (i,j) to generate from.
! This uses a triangular mapping to pick them uniformly.
call pick_uniform_elecs(elecs, pelec)
eleci = elecs(1)
elecj = elecs(2)
! Obtain the orbitals and their momentum vectors for the given elecs.
orbi = nI(eleci)
orbj = nI(elecj)
ex(1, 1) = orbi
ex(1, 2) = orbj
! Loop through all available orbitals, A. If it is unoccupied, then
! find the orbital, B, that will complete the excitation given i,j.
! If this is also unoccupied, then contribute to the cumulative list
! for making selections
if (TContact) then
call create_ab_list_ua(nI, ilutI, [orbi, orbj], cum_arr, cum_sum)
else
call create_ab_list_ueg(ilutI, [orbi, orbj], cum_arr, cum_sum)
end if
! If there are no available excitations, then we need to reject this
! excitation
if (cum_sum < EPS) then
nJ(1) = 0
pgen = 0.0_dp
return
end if
! Pick a pair of orbitals A,B according to the cumulative probabilites
! already generated.
r = genrand_real2_dSFMT() * cum_sum
orba = binary_search_first_ge(cum_arr, r)
! write a cleaner implementation
orbb = get_orb_from_kpoints(orbi, orbj, orba)
! Calculate the orbital selection probability
elem = cum_arr(orba)
if (orba > 1) elem = elem - cum_arr(orba - 1)
porb = elem / cum_sum
! Construct and return the determinant
call make_double(nI, nJ, eleci, elecj, orba, orbb, ex, tPar)
ilutJ = ilutI
clr_orb(ilutJ, orbi)
clr_orb(ilutJ, orbj)
set_orb(ilutj, orba)
set_orb(ilutj, orbb)
pgen = pelec * porb
end subroutine gen_double_ueg