PropVec_doubles_PCHB_spinorb_gen_exc Subroutine

    subroutine PropVec_doubles_PCHB_spinorb_gen_exc(&
                this, nI, ilutI, nJ, ilutJ, exFlag, ic, &
                ex, tParity, pGen, hel, store, part_type)
        !! given the initial determinant (both as nI and ilut), create a random
        !! doubles excitation using the Hamiltonian matrix elements as weights
        class(PropVec_PCHB_DoublesSpinOrbFastWeightExcGen_t), intent(inout) :: this
            !! the exctitation generator
        integer, intent(in) :: nI(nel), exFlag
            !! determinant to excite from
            !! unused in this generator
        integer(n_int), intent(in) :: ilutI(0:NIfTot)
            !! determint from which to excite, ilut format
        integer, intent(out) :: nJ(nel), ic, ex(2, maxExcit)
            !! the excited determinant upon return
            !! excitation order (for doubles generator, always == 2)
            !! excitation matrix nI -> nJ
        integer(n_int), intent(out) :: ilutJ(0:NifTot)
            !! excited determinant, ilut format
        real(dp), intent(out) :: pGen
            !! probability of generating the excitation
        logical, intent(out) :: tParity
            !! the parity of the excitation nI -> nJ
        HElement_t(dp), intent(out) :: hel
            !! matrix element Hijab
        type(excit_gen_store_type), intent(inout), target :: store
            !!
        integer, intent(in), optional :: part_type
            !! unused in this generator
        character(*), parameter :: this_routine = 'PropVec_doubles_PCHB_spinorb_gen_exc'

        integer :: i_sg ! prop_vec index
        integer :: src(2) ! particles (I, J)
        integer :: elecs(2) ! particle indices in nI
        integer :: tgt(2)

        integer :: ij
        logical :: invalid

        integer :: ab
        real(dp) :: pGenHoles


#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._dp) ! macro
#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

        ! pick two random electrons
        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_spinorb_fastweighted.fpp:188"
            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.
        ij = fuse_symm_idx(src(1), src(2))

        ! get a pair of orbitals using the precomputed weights
        call this%AB_sampler%sample(ij, 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%AB_sampler%get_prob(ij, i_sg, ab)))) then
            write(stderr, *) ""
            write(stderr, *) "Assertion ab /= 0 .implies. (pGenHoles .isclose. this%AB_sampler%get_prob(ij, 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_spinorb_fastweighted.fpp:195"
            call stop_all (this_routine, "Assert fail: ab /= 0 .implies. (pGenHoles .isclose. this%AB_sampler%get_prob(ij, i_sg,&
                & ab))")
        end if
    end block
#endif

        if (ab == 0) then
            invalid = .true.
            tgt = 0
        else
            ! split ab -> a,b
            tgt = this%tgtOrbs(:, ab)
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        use constants, only: stderr
        if (.not. (all(tgt /= 0) .implies. tgt(1) /= tgt(2))) then
            write(stderr, *) ""
            write(stderr, *) "Assertion all(tgt /= 0) .implies. tgt(1) /= tgt(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_spinorb_fastweighted.fpp:203"
            call stop_all (this_routine, "Assert fail: all(tgt /= 0) .implies. tgt(1) /= tgt(2)")
        end if
    end block
#endif
            invalid = any(tgt == 0) .or. any(tgt(1) == nI) .or. any(tgt(2) == nI)
        end if

        ! as in the spatially-resolved case, there is a very rare case where (due
        ! to floating point error) an excitation with pgen=0 is created. Invalidate.
        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 returned, no excitations -> nulldet
            call invalidate()
        else
            ! construct the determinant given the selected orbitals and electrons
            call make_double(nI, nJ, elecs(1), elecs(2), tgt(1), tgt(2), ex, tParity)

            ilutJ = exciteIlut(ilutI, src, tgt)

            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_spinorb_fastweighted.fpp:228"
            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) = tgt
        end subroutine invalidate

    end subroutine PropVec_doubles_PCHB_spinorb_gen_exc