PropVec_doubles_PCHB_get_pgen Function

    function PropVec_doubles_PCHB_get_pgen(this, nI, ilutI, ex, ic, ClassCount2, ClassCountUnocc2) result(pgen)
        class(PropVec_PCHB_DoublesSpatOrbFastWeightExcGen_t), intent(inout) :: this
        integer, intent(in) :: nI(nel)
        integer(n_int), intent(in) :: ilutI(0:NIfTot)
        integer, intent(in) :: ex(2, maxExcit), ic
        integer, intent(in) :: ClassCount2(ScratchSize), ClassCountUnocc2(ScratchSize)
        real(dp) :: pgen
        character(*), parameter :: this_routine = 'PropVec_doubles_PCHB_get_pgen'

        integer :: ab, ij, nex(2, 2), samplerIndex, i_sg

#ifdef WARNING_WORKAROUND_
        associate(ilutI => ilutI); end associate
        associate(ClassCount2 => ClassCount2); end associate
        associate(ClassCountUnocc2 => ClassCountUnocc2); end associate
#endif
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        use constants, only: stderr
        if (.not. (ic == 2)) then
            write(stderr, *) ""
            write(stderr, *) "Assertion ic == 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:295"
            call stop_all (this_routine, "Assert fail: ic == 2")
        end if
    end block
#endif

        nex = gtID(ex(:, : ic))
        ij = fuse_symm_idx(nex(1, 1), nex(1, 2))
        ab = fuse_symm_idx(nex(2, 1), nex(2, 2))
        i_sg = this%indexer%idx_nI(nI)

        pgen = this%particle_selector%get_pgen(nI, i_sg, ex(1, 1), ex(1, 2))

        ! the probability of picking the two electrons: they are chosen uniformly
        ! check which sampler was used
        if (is_beta(ex(1, 1)) .eqv. is_beta(ex(1, 2))) then
            ! same-spin case
            samplerIndex = SAME_SPIN
        else
            ! excitations without spin-exchange OR to the same spatial orb
            if ((is_beta(ex(1, 1)) .eqv. is_beta(ex(2, 1))) .or. (nex(2, 1) == nex(2, 2))) then
                ! opp spin case without exchange
                samplerIndex = OPP_SPIN_NO_EXCH
                pGen = pGen * (1.0_dp - this%pExch(ij, i_sg))
            else
                ! opp spin case with exchange
                samplerIndex = OPP_SPIN_EXCH
                pGen = pGen * this%pExch(ij, i_sg)
            end if
        end if

        pgen = pgen * this%pchb_samplers%get_prob(ij, samplerIndex, i_sg, ab)
    end function PropVec_doubles_PCHB_get_pgen