pick_orbitals_double_pchb Subroutine

    subroutine pick_orbitals_double_pchb(this, ilut, nI, csf_i, excitInfo, pgen)
        debug_function_name("pick_orbitals_double_pchb")
        class(GugaAliasSampler_t), intent(in) :: this
        integer(n_int), intent(in) :: ilut(0:GugaBits%len_tot)
        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 :: src(2), sym_prod, sum_ml, i, j, a, b, orbs(2), ij, ab
        real(dp) :: pgen_elec, pgen_orbs

        unused_var(ilut)

        ! maybe I will also produce a weighted electron pickin in the
        ! GUGA formalism.. but for now pick them uniformly:
        call pick_elec_pair_uniform_guga(nI, src, sym_prod, sum_ml, pgen_elec)
        ! make a different picker, which does not bias towards doubly
        ! occupied orbitals here too!

        ASSERT( src(1) < src(2) )

        i = gtID(src(1))
        j = gtID(src(2))

        if (i /= j) then
            if (csf_i%stepvector(i) == 3) pgen_elec = 2.0_dp * pgen_elec
            if (csf_i%stepvector(j) == 3) pgen_elec = 2.0_dp * pgen_elec
        end if

        ! use the sampler for this electron pair -> order of src electrons
        ! does not matter
        ij = fuseIndex(i, j)

        ! get a pair of orbitals using the precomputed weights
        call this%alias_sampler%sample(ij,ab,pgen_orbs)

        ! unfortunately, there is a super-rare case when, due to floating point error,
        ! an excitation with pGen=0 is created. Those are invalid, too
        if(near_zero(pgen_orbs)) then
            excitInfo%valid = .false.
            pgen = 0.0_dp
            ! Yes, print. Those events are signficant enough to be always noted in the output
            print *, "WARNING: Generated excitation with probability of 0"
            return
        endif


        ! split the index ab (using a table containing mapping ab -> (a,b))
        orbs = this%tgtOrbs(:,ab)

        a = orbs(1)
        b = orbs(2)

        ! check if the picked orbs are a valid choice - if they are the same, match one
        ! occupied orbital or are zero (maybe because there are no allowed picks for
        ! the given source) abort
        ! these are the 'easy' checks for GUGA.. more checks need to be done
        ! to see if it is actually a valid combination..
        if (any(orbs == 0) &
                .or. csf_i%stepvector(a) == 3 &
                .or. csf_i%stepvector(b) == 3 &
                .or. a == b .and. csf_i%stepvector(b) /= 0) then
            excitInfo%valid = .false.
            return
        end if

        ! setup a getInfo functionality in the sampler!
        call extract_excit_info(this%get_info_entry(ij, ab), excitInfo)

        pGen = pgen_elec * pgen_orbs

    end subroutine pick_orbitals_double_pchb