pickOrbs_sym_uniform_ueg_double Subroutine

    subroutine pickOrbs_sym_uniform_ueg_double(ilut, nI, csf_i, excitInfo, pgen)
        ! specific orbital picker for hubbard and UEG type models with
        ! full k-point symmetry
        integer(n_int), intent(in) :: ilut(0:nifguga)
        integer, intent(in) :: nI(nel)
        type(CSF_Info_t), intent(in) :: csf_i
        type(ExcitationInformation_t), intent(out) :: excitInfo
        real(dp), intent(out) :: pgen

        integer :: occ_orbs(2), ind, eleci, elecj, orb, orb_2, orb_arr(nSpatOrbs)
        real(dp) :: cum_sum, cum_arr(nSpatOrbs), pelec, r, cpt1, cpt2
        type(ExcitationInformation_t) :: excit_arr(nBasis)

        ! pick 2 electrons uniformly first:

        ! or since other symmetries are usually not used in the hubbard just:
        ! Pick a pair of electrons (i,j) to generate from.
        ! This uses a triangular mapping to pick them uniformly.
        ind = 1 + int(ElecPairs * genrand_real2_dSFMT())
        elecj = ceiling((1 + sqrt(1 + 8 * real(ind, dp))) / 2)
        eleci = ind - ((elecj - 1) * (elecj - 2)) / 2
        pelec = 1.0_dp / real(ElecPairs, dp)

        ! i pick spatial orbitals out of occupied spin-orbitals ->
        ! so the chance to pick a doubly occupied spatial orbital is twice
        ! as high -> adjust the corresponding probabilities

        ! Obtain the orbitals and their momentum vectors for the given elecs.
        occ_orbs(1) = nI(eleci)
        occ_orbs(2) = nI(elecj)

        call gen_ab_cum_list_ueg(ilut, csf_i, occ_orbs, cum_arr, excit_arr, orb_arr)

        ! then pick a orbital randomly and consider a <> b contribution
        cum_sum = cum_arr(nSpatOrbs)

        if (near_zero(cum_sum)) then
            excitInfo%valid = .false.
            return
        end if

        r = genrand_real2_dSFMT() * cum_sum
        orb = binary_search_first_ge(cum_arr, r)

        ! pick the info
        excitInfo = excit_arr(orb)

        ! and pgens, for that i need second orbital too..
        orb_2 = orb_arr(orb)
        if (orb == 1) then
            cpt1 = cum_arr(1)
        else
            cpt1 = cum_arr(orb) - cum_arr(orb - 1)
        end if

        cpt2 = 0.0_dp
        if (orb /= orb_2) then
            if (orb_2 == 1) then
                cpt2 = cum_arr(1)
            else
                cpt2 = cum_arr(orb_2) - cum_arr(orb_2 - 1)
            end if
        end if
        pgen = pelec * (cpt1 + cpt2) / cum_sum

        if (.not. is_in_pair(occ_orbs(1), occ_orbs(2))) then
            if (csf_i%stepvector(gtID(occ_orbs(1))) == 3) pgen = pgen * 2.0_dp
            if (csf_i%stepvector(gtID(occ_orbs(2))) == 3) pgen = pgen * 2.0_dp
        end if

    end subroutine pickOrbs_sym_uniform_ueg_double