PropVec_doubles_PCHB_gen_exc Subroutine

    subroutine PropVec_doubles_PCHB_gen_exc(&
                    this, nI, ilutI, nJ, ilutJ, exFlag, ic, &
                    ex, tParity, pGen, hel, store, part_type)
        class(PropVec_PCHB_DoublesSpatOrbFastWeightExcGen_t), intent(inout) :: this
        integer, intent(in) :: nI(nel), exFlag
        integer(n_int), intent(in) :: ilutI(0:NIfTot)
        integer, intent(out) :: nJ(nel), ic, ex(2, maxExcit)
        integer(n_int), intent(out) :: ilutJ(0:NifTot)
        real(dp), intent(out) :: pGen
        logical, intent(out) :: tParity
        HElement_t(dp), intent(out) :: hel
        type(excit_gen_store_type), intent(inout), target :: store
        integer, intent(in), optional :: part_type
        character(*), parameter :: this_routine = 'PropVec_doubles_PCHB_gen_exc'

        integer :: elecs(2), src(2), ij
        integer :: orbs(2), ab
        real(dp) :: pGenHoles
        logical :: invalid
        integer :: spin(2), samplerIndex
        integer :: i_sg

#ifdef WARNING_WORKAROUND_
        associate(exFlag => exFlag); end associate
        associate(part_type => part_type); end associate
#endif
        ic = 2
#ifdef WARNING_WORKAROUND_
        hel = h_cast(0.0_dp)
#endif

        if (this%use_lookup) then
            i_sg = this%indexer%lookup_prop_vec_idx(store%idx_curr_dets, nI)
        else
            i_sg = this%indexer%idx_nI(nI)
        end if


        ! first, pick two random elecs
        call this%particle_selector%draw(nI, ilutI, i_sg, elecs, src, pGen)
        if (src(1) == 0) then
            call invalidate()
            return
        end if


#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        use constants, only: stderr
        if (.not. (pGen .isclose. this%particle_selector%get_pgen(nI, i_sg, src(1), src(2)))) then
            write(stderr, *) ""
            write(stderr, *) "Assertion pGen .isclose. this%particle_selector%get_pgen(nI, i_sg, src(1), src(2))"
            write(stderr, *) "failed in /home/cicdadmin/fkfest-gh-runners/altest-runner-7/_work/NECI_autoforward/NECI_autoforward/s&
                &rc/property_vector_constraint/property_vector_pchb_doubles_spatorb_fastweighted.fpp:195"
            call stop_all (this_routine, "Assert fail: pGen .isclose. this%particle_selector%get_pgen(nI, i_sg, src(1), src(2))")
        end if
    end block
#endif

        invalid = .false.
        ! use the sampler for this electron pair -> order of src electrons does not matter
        ij = fuse_symm_idx(gtID(src(1)), gtID(src(2)))
        ! the spin of the electrons: 0 - alpha, 1 - beta
        spin = getSpinIndex(src)
        ! determine type of spin-excitation: same-spin, opp spin w exchange, opp spin w/o exchange
        if (spin(1) == spin(2)) then
            ! same spin
            samplerIndex = SAME_SPIN
        else
            ! else, pick exchange with...some ij-spin bias
            if (genrand_real2_dSFMT() < this%pExch(ij, i_sg)) then
                samplerIndex = OPP_SPIN_EXCH
                ! adjust pgen
                pGen = pGen * this%pExch(ij, i_sg)
                ! the spins of the target are the opposite of the source spins
                call swap(spin(1), spin(2))
            else
                samplerIndex = OPP_SPIN_NO_EXCH
                ! adjust pgen
                pGen = pGen * (1.0_dp - this%pExch(ij, i_sg))
            end if
        end if
        ! get a pair of orbitals using the precomputed weights
        call this%pchb_samplers%sample(ij, samplerIndex, i_sg, ab, pGenHoles)
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        use constants, only: stderr
        if (.not. (ab /= 0 .implies. (pGenHoles .isclose.  this%pchb_samplers%get_prob(ij, samplerIndex, i_sg, ab)))) then
            write(stderr, *) ""
            write(stderr, *) "Assertion ab /= 0 .implies. (pGenHoles .isclose.  this%pchb_samplers%get_prob(ij, samplerIndex,&
                & i_sg, ab))"
            write(stderr, *) "failed in /home/cicdadmin/fkfest-gh-runners/altest-runner-7/_work/NECI_autoforward/NECI_autoforward/s&
                &rc/property_vector_constraint/property_vector_pchb_doubles_spatorb_fastweighted.fpp:222"
            call stop_all (this_routine, "Assert fail: ab /= 0 .implies. (pGenHoles .isclose.  this%pchb_samplers%get_prob(ij,&
                & samplerIndex, i_sg, ab))")
        end if
    end block
#endif


        if (ab == 0) then
            invalid = .true.
            orbs = 0
        else
            ! split the index ab (using a table containing mapping ab -> (a,b))
            orbs = this%tgtOrbs(:, ab)
            ! convert orbs to spin-orbs with the same spin
            orbs = 2 * orbs - spin
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        use constants, only: stderr
        if (.not. (all(orbs /= 0) .implies. orbs(1) /= orbs(2))) then
            write(stderr, *) ""
            write(stderr, *) "Assertion all(orbs /= 0) .implies. orbs(1) /= orbs(2)"
            write(stderr, *) "failed in /home/cicdadmin/fkfest-gh-runners/altest-runner-7/_work/NECI_autoforward/NECI_autoforward/s&
                &rc/property_vector_constraint/property_vector_pchb_doubles_spatorb_fastweighted.fpp:233"
            call stop_all (this_routine, "Assert fail: all(orbs /= 0) .implies. orbs(1) /= orbs(2)")
        end if
    end block
#endif
            ! note: nI is an array, not a scalar, so we need two `any` checks below
            invalid = any(orbs == 0) .or. any(orbs(1) == nI) .or. any(orbs(2) == nI)
        end if

        ! 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 (.not. invalid .and. near_zero(pGenHoles)) then
            invalid = .true.
            ! Yes, print. Those events are signficant enough to be always noted in the output
            write(stdout, *) "WARNING: Generated excitation with probability of 0"
        end if

        if (invalid) then
            ! if 0 is returned, there are no excitations for the chosen elecs
            ! -> return nulldet
            call invalidate()
        else
            ! else, construct the det from the chosen orbs/elecs

            call make_double(nI, nJ, elecs(1), elecs(2), orbs(1), orbs(2), ex, tParity)

            ilutJ = exciteIlut(ilutI, src, orbs)

            pGen = pGen * pGenHoles

            block
                integer :: ClassCount2(ScratchSize), ClassCountUnocc2(ScratchSize)
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        use constants, only: stderr
        if (.not. (pgen .isclose. this%get_pgen(nI, ilutI, ex, ic, ClassCount2, ClassCountUnocc2))) then
            write(stderr, *) ""
            write(stderr, *) "Assertion pgen .isclose. this%get_pgen(nI, ilutI, ex, ic, ClassCount2, ClassCountUnocc2)"
            write(stderr, *) "failed in /home/cicdadmin/fkfest-gh-runners/altest-runner-7/_work/NECI_autoforward/NECI_autoforward/s&
                &rc/property_vector_constraint/property_vector_pchb_doubles_spatorb_fastweighted.fpp:261"
            call stop_all (this_routine, "Assert fail: pgen .isclose. this%get_pgen(nI, ilutI, ex, ic, ClassCount2,&
                & ClassCountUnocc2)")
        end if
    end block
#endif
            end block

        end if

    contains

        subroutine invalidate()
            nJ = 0
            ilutJ = 0_n_int
            ex(1, 1 : 2) = src
            ex(2, 1 : 2) = orbs
        end subroutine invalidate
    end subroutine PropVec_doubles_PCHB_gen_exc