GAS_doubles_PCHB_get_pgen Function

    function GAS_doubles_PCHB_get_pgen(this, nI, ilutI, ex, ic, ClassCount2, ClassCountUnocc2) result(pgen)
        class(GAS_PCHB_DoublesSpinorbFullyWeightedExcGenerator_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
        real(dp) :: pgen_particle
        character(*), parameter :: this_routine = 'GAS_doubles_PCHB_get_pgen'

        integer :: IJ, i_sg
        integer :: unoccupied(nBasis - nEl)
        real(dp) :: renorm_first, p_first(2), renorm_second(2), p_second(2)

#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 /scratch/jenkins/jobs/existing_branch_doc/workspace/build_config/gfortran-doc/src/gasci_pch&
                &b_doubles_spinorb_fullyweighted.fpp:275"
            call stop_all (this_routine, "Assert fail: ic == 2")
        end if
    end block
#endif

        IJ = fuseIndex(ex(1, 1), ex(1, 2))
        i_sg = this%indexer%idx_nI(nI)

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

        block
            integer(n_int) :: ilut_unoccupied(0 : nIfD)
            call this%get_unoccupied(ilutI(0 : nIfD), ilut_unoccupied, unoccupied)
        end block
        ! The renormalization for the first hole is the same, regardless of order.
        renorm_first = 1._dp - sum(this%A_sampler%get_prob(IJ, i_sg, nI))
        if (do_direct_calculation(renorm_first)) then
            renorm_first = sum(this%A_sampler%get_prob(IJ, i_sg, unoccupied))
        end if

        associate(A => ex(2, 1), B => ex(2, 2))
            p_first(1) = this%A_sampler%constrained_getProb(IJ, i_sg, unoccupied, renorm_first, A)
            p_first(2) = this%A_sampler%constrained_getProb(IJ, i_sg, unoccupied, renorm_first, B)

            renorm_second(1) = 1._dp - sum(this%B_sampler%get_prob(A, IJ, i_sg, nI))
            if (do_direct_calculation(renorm_second(1))) then
                renorm_second(1) = sum(this%B_sampler%get_prob(A, IJ, i_sg, unoccupied))
            end if
            p_second(1) = this%B_sampler%constrained_getProb(&
                A, IJ, i_sg, unoccupied, renorm_second(1), B)

            renorm_second(2) = 1._dp - sum(this%B_sampler%get_prob(B, IJ, i_sg, nI))
            if (do_direct_calculation(renorm_second(2))) then
                renorm_second(2) = sum(this%B_sampler%get_prob(B, IJ, i_sg, unoccupied))
            end if
            p_second(2) = this%B_sampler%constrained_getProb(&
                B, IJ, i_sg, unoccupied, renorm_second(2), A)
        end associate

        pgen = pgen_particle * sum(p_first * p_second)

    end function GAS_doubles_PCHB_get_pgen