#include "macros.h"
module guga_prop_vec_pchb_singles_main
    use util_mod, only: EnumBase_t, stop_all, operator(.implies.), &
        operator(.div.), operator(.isclose.), near_zero
    use SystemData, only: nEl, nSpatOrbs, nBasis
    use MPI_wrapper, only: root
    use constants, only: n_int, dp, bits_n_int, stdout
    use dSFMT_interface, only: genrand_real2_dSFMT
    use property_vector_index, only: AlsoGUGA_PropertyIndexer_t, lookup_property_indexer
    use excitation_types, only: Excite_1_t
    use fortran_strings, only: to_upper
    use bit_rep_data, only: nIfD, GugaBits
    use bit_reps, only: decode_bit_det, UnoccupiedGetter_t
    use aliasSampling, only: AliasSampler_1D_t, AliasSampler_2D_t, do_direct_calculation, AliasSampler_3D_t
    use UMatCache, only: gtID, get_umat_el

    use guga_base_class, only: SinglesGUGABase_t
    use guga_bitrepops, only: CSF_Info_t
    use guga_data, only: ExcitationInformation_t, gen_type
    use guga_pchb_singles_weights_mod, only: get_weight_t, PropVec_PCHB_SinglesWeighting_t, &
        PropVec_PCHB_SinglesWeighting_vals_t, set_get_weight_pointer
    use guga_excitations, only: assign_excitinfo_values_single, &
        checkcompatibility_single
    better_implicit_none


    private
    public :: allocate_and_init, PropVec_PCHB_SinglesOptions_t, PropVec_PCHB_SinglesOptions_vals_t, &
        scale_down_exch_singles

    type, extends(EnumBase_t) :: PropVec_PCHB_SinglesAlgorithm_t
        character(9) :: str
    contains
        procedure :: to_str => to_str_algorithm
    end type

    real(dp) :: scale_down_exch_singles = 1._dp

    type :: PropVec_PCHB_SinglesAlgorithm_vals_t
        type(PropVec_PCHB_SinglesAlgorithm_t) :: &
            UNIF_UNIF = PropVec_PCHB_SinglesAlgorithm_t(1, 'UNIF:UNIF'), &
            FULL_FULL = PropVec_PCHB_SinglesAlgorithm_t(2, 'FULL:FULL'), &
            UNIF_FULL = PropVec_PCHB_SinglesAlgorithm_t(3, 'UNIF:FULL'), &
            UNIF_XNEW = PropVec_PCHB_SinglesAlgorithm_t(4, 'UNIF:XNEW')
    contains
        procedure, nopass :: from_str => alg_from_str
    end type

    type :: PropVec_PCHB_SinglesOptions_t
        type(PropVec_PCHB_SinglesAlgorithm_t) :: algorithm
        type(PropVec_PCHB_SinglesWeighting_t) :: weighting
    contains
        procedure :: is_valid
        procedure :: to_str => to_str_singles
    end type

    type :: PropVec_PCHB_SinglesOptions_vals_t
        type(PropVec_PCHB_SinglesAlgorithm_vals_t) :: &
            algorithm = PropVec_PCHB_SinglesAlgorithm_vals_t()
        type(PropVec_PCHB_SinglesWeighting_vals_t) :: &
            weighting = PropVec_PCHB_SinglesWeighting_vals_t()
    end type

    type(PropVec_PCHB_SinglesOptions_vals_t), parameter :: &
        option_vals = PropVec_PCHB_SinglesOptions_vals_t()


    !> The precomputed PropVec uniform excitation generator
    type, extends(SinglesGUGABase_t) :: PropVec_UniformExcGenerator_t
        private
        class(AlsoGUGA_PropertyIndexer_t), allocatable :: indexer
        logical, public :: use_lookup = .false.
            !! Use a lookup for the prop_vec index in global_det_data.
        logical, public :: create_lookup = .false.
            !! Create **and** manage! the prop_vec index lookup in global_det_data.

        integer(n_int), allocatable :: allowed_holes(:, :, :)
            !! Bitmask for the allowed holes for a single excitation
        integer :: L_spat_bits = -1
            !! The number of  integer(n_int) numbers to
            !! store information about spatial orbitals as bitmask
    contains
        private
        procedure, public :: init => init_PropVec_UniformExcGenerator_t
        procedure, public :: finalize => PropVec_singles_uniform_finalize

        procedure :: get_valid_orbs

        procedure, public :: calc_pgen => calc_pgen_uniform_singles
        procedure, public :: pickOrbitals => pickOrbitals_single
    end type

    !> The precomputed PropVec uniform excitation generator
    type, extends(SinglesGUGABase_t), abstract :: PropVec_Weighted_t
        private
        class(AlsoGUGA_PropertyIndexer_t), allocatable :: indexer
        logical, public :: use_lookup = .false.
            !! Use a lookup for the prop_vec index in global_det_data.
        logical, public :: create_lookup = .false.
            !! Create **and** manage! the prop_vec index lookup in global_det_data.


        type(AliasSampler_1D_t) :: I_sampler
            !! The shape is (n_prop_vec) -> 2 * nSpatOrbs
            !!
            !! @note
            !! To make doubly occupied orbitals more probable we have doubly the
            !! amount of spatial orbitals and treat it as if there were spin-orbitals.
            !! @endnote
            !!
        type(AliasSampler_2D_t) :: A_sampler
            !! The shape is (nSpatOrbs, n_prop_vec) -> 2 * nSpatOrbs
            !!
            !! @note
            !! To make empty orbitals more probable than singly occpied ones,
            !! we have doubly the
            !! amount of spatial orbitals and treat it as if there were spin-orbitals.
            !! @endnote
            !!
        procedure(get_weight_t), pointer, nopass :: get_weight => null()

        type(UnoccupiedGetter_t) :: unoccupied

    contains
        private
        procedure, public :: init => init_PropVec_Weighted_t
        procedure, public :: finalize => finalize_PropVec_Weighted_t
    end type

    type, extends(PropVec_Weighted_t) :: SinglesPropVec_FullFull_t
    contains
        procedure, public :: calc_pgen => calc_pgen_FullFull
        procedure, public :: pickOrbitals => pickOrbitals_FullFull
    end type

    type, extends(PropVec_Weighted_t) :: SinglesPropVec_UnifFull_t
    contains
        procedure, public :: calc_pgen => calc_pgen_UnifFull
        procedure, public :: pickOrbitals => pickOrbitals_UnifFull
    end type

    !> The precomputed PropVec uniform excitation generator
    type, extends(SinglesGUGABase_t) :: UnifXnew_t
        private
        class(AlsoGUGA_PropertyIndexer_t), allocatable :: indexer
        logical, public :: use_lookup = .false.
        logical, public :: create_lookup = .false.


        type(AliasSampler_3D_t) :: A_sampler
            !! The shape is (nSpatOrbs: i, stepvector: i, n_prop_vec: i_sg) -> 2 * nSpatOrbs
        procedure(get_weight_t), pointer, nopass :: get_weight => null()

        type(UnoccupiedGetter_t) :: unoccupied

    contains
        private
        procedure, public :: init => init_UnifXnew_t
        procedure, public :: finalize => finalize_UnifXnew_t

        procedure, public :: pickOrbitals => pickOrbitals_UnifXnew_t
        procedure, public :: calc_pgen => calc_pgen_UnifXnew_t
    end type

contains

    elemental logical function is_valid(this)
        class(PropVec_PCHB_SinglesOptions_t), intent(in) :: this
        is_valid = (this%algorithm == option_vals%algorithm%UNIF_UNIF &
                            .implies. this%weighting == option_vals%weighting%UNIFORM)
    end function


    subroutine allocate_and_init(indexer, options, use_lookup, create_lookup, generator)
        class(AlsoGUGA_PropertyIndexer_t), intent(in) :: indexer
        type(PropVec_PCHB_SinglesOptions_t), intent(in) :: options
        logical, intent(in) :: use_lookup, create_lookup
        class(SinglesGUGABase_t), allocatable, intent(inout) :: generator
        routine_name("allocate_and_init")

#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        use constants, only: stderr
        if (.not. (options%is_valid())) then
            write(stderr, *) ""
            write(stderr, *) "Assertion options%is_valid()"
            write(stderr, *) "failed in /home/cicdadmin/fkfest-gh-runners/altest-runner-7/_work/NECI_autoforward/NECI_autoforward/s&
                &rc/GUGA/guga_propvec_pchb_singles_main.fpp:182"
            call stop_all (this_routine, "Assert fail: options%is_valid()")
        end if
    end block
#endif

        if (allocated(generator)) then
            call generator%finalize()
            deallocate(generator)
        end if

        if (options%algorithm == option_vals%algorithm%UNIF_UNIF) then
            allocate(PropVec_UniformExcGenerator_t :: generator)
            select type(generator)
            type is(PropVec_UniformExcGenerator_t)
                call generator%init(indexer, use_lookup, create_lookup)
            end select
        else if (options%algorithm == option_vals%algorithm%FULL_FULL) then
            allocate(SinglesPropVec_FullFull_t :: generator)
            select type(generator)
            type is(SinglesPropVec_FullFull_t)
                call generator%init(options%weighting, indexer, use_lookup, create_lookup)
            end select
        else if (options%algorithm == option_vals%algorithm%UNIF_FULL) then
            allocate(SinglesPropVec_UnifFull_t :: generator)
            select type(generator)
            type is(SinglesPropVec_UnifFull_t)
                call generator%init(options%weighting, indexer, use_lookup, create_lookup)
            end select
        else if (options%algorithm == option_vals%algorithm%UNIF_XNEW) then
            allocate(UnifXnew_t :: generator)
            select type(generator)
            type is(UnifXnew_t)
                call generator%init(options%weighting, indexer, use_lookup, create_lookup, scale_down_exch_singles)
            end select
        else
            call stop_all(this_routine, "unsupported singles option.")
        end if
    end subroutine



    subroutine init_PropVec_UniformExcGenerator_t(this, indexer, use_lookup, create_lookup)
        class(PropVec_UniformExcGenerator_t), intent(inout) :: this
        class(AlsoGUGA_PropertyIndexer_t), intent(in) :: indexer
        logical, intent(in) :: use_lookup, create_lookup
        integer, allocatable :: prop_vecs(:, :)
        routine_name('init_PropVec_UniformExcGenerator_t')


        integer :: i_sg, src, tgt

        allocate(this%indexer, source=indexer)
        this%create_lookup = create_lookup
        if (create_lookup) then
            if (allocated(lookup_property_indexer)) then
                call stop_all(this_routine, 'Someone else is already managing the prop_vec lookup.')
            else
                write(stdout, *) 'PropVec singles is creating and managing the prop_vec lookup'
                allocate(lookup_property_indexer, source=this%indexer)
            end if
        end if
        this%use_lookup = use_lookup
        if (use_lookup) write(stdout, *) 'PropVec singles is using the prop_vec lookup'
        ! possible prop_vecs
        prop_vecs = this%indexer%get_prop_vecs()

        this%L_spat_bits = nSpatOrbs .div. bits_n_int

        allocate(this%allowed_holes(0:this%L_spat_bits, nSpatOrbs, size(prop_vecs, 2)), &
                 source=0_n_int)


        do i_sg = 1, size(prop_vecs, 2)
            do src = 1, nSpatOrbs
                do tgt = 1, nSpatOrbs
                block
                    type(Excite_1_t) :: exc
                    exc = Excite_1_t(2 * src, 2 * tgt)
                    if (src /= tgt &
                            .and. this%indexer%is_allowed(exc, prop_vecs(:, i_sg)) &
                            .and. symmetry_allowed(exc)) then
                        call my_set_orb(this%allowed_holes(:, src, i_sg), tgt)
                    end if
                end block
                end do
            end do
        end do

    contains

        ! Ugly Fortran syntax rules forces us to not use the macro.
        ! Even associate would not help here
        ! https://stackoverflow.com/questions/65734764/non-one-indexed-array-from-associate
        pure subroutine my_set_orb(ilut, orb)
            integer(n_int), intent(inout) :: ilut(0:this%L_spat_bits)
            integer, intent(in) :: orb
            set_orb(ilut, orb)
        end subroutine

    end subroutine

    logical pure function symmetry_allowed(exc)
        use SymExcitDataMod, only: SpinOrbSymLabel
        type(Excite_1_t), intent(in) :: exc
        symmetry_allowed = SpinOrbSymLabel(exc%val(1, 1)) == SpinOrbSymLabel(exc%val(2, 1))
    end function

    subroutine PropVec_singles_uniform_finalize(this)
        class(PropVec_UniformExcGenerator_t), intent(inout) :: this
        deallocate(this%allowed_holes, this%indexer)
        if (this%create_lookup) deallocate(lookup_property_indexer)
    end subroutine


    subroutine pickOrbitals_single(this, nI, ilut, csf_i, excitInfo, pgen)
        debug_function_name("pickOrbitals_single")
        class(PropVec_UniformExcGenerator_t), intent(in) :: this
        integer, intent(in) :: nI(nel)
        integer(n_int), intent(in) :: ilut(0:GugaBits%len_tot)
        type(CSF_Info_t), intent(in) :: csf_i
        type(ExcitationInformation_t), intent(out) :: excitInfo
        real(dp), intent(out) :: pgen

        integer :: src_orb, tgt_orb, i
        integer, allocatable :: occ_spat_orbs(:)
            !! The occupied spatial orbitals
        integer, allocatable :: valid_orbs(:)
        unused_var(ilut); unused_var(nI)

        occ_spat_orbs = pack([(i, i=1, size(csf_i%Occ_int))], csf_i%Occ_int /= 0)
        src_orb = occ_spat_orbs(1 + int(genrand_real2_dSFMT() * size(occ_spat_orbs)))

        valid_orbs = this%get_valid_orbs(csf_i, src_orb)

        if (size(valid_orbs) == 0) then
            pgen = 0.0_dp
            tgt_orb = 0
            excitInfo%valid = .false.
        else
            pgen = 1._dp / real(size(occ_spat_orbs) * size(valid_orbs), dp)
            tgt_orb = valid_orbs(1 + floor(genrand_real2_dSFMT() * size(valid_orbs)))


            if (tgt_orb < src_orb) then
                excitInfo = assign_excitinfo_values_single( &
                            gen_type%R, tgt_orb, src_orb, tgt_orb, src_orb)
            else
                excitInfo = assign_excitinfo_values_single( &
                            gen_type%L, tgt_orb, src_orb, src_orb, tgt_orb)
            end if
            excitInfo%valid = checkCompatibility_single(csf_i, excitInfo)
            if (excitInfo%valid) then
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        use constants, only: stderr
        if (.not. (pgen .isclose. this%calc_pgen(nI, ilut, csf_i, excitInfo))) then
            write(stderr, *) ""
            write(stderr, *) "Assertion pgen .isclose. this%calc_pgen(nI, ilut, csf_i, excitInfo)"
            write(stderr, *) "failed in /home/cicdadmin/fkfest-gh-runners/altest-runner-7/_work/NECI_autoforward/NECI_autoforward/s&
                &rc/GUGA/guga_propvec_pchb_singles_main.fpp:331"
            call stop_all (this_routine, "Assert fail: pgen .isclose. this%calc_pgen(nI, ilut, csf_i, excitInfo)")
        end if
    end block
#endif
            else
                pgen = 0._dp
            end if
        end if
    end subroutine pickOrbitals_single


    function get_valid_orbs(this, csf_i, src_orb) result(valid_orbs)
        class(PropVec_UniformExcGenerator_t), intent(in) :: this
        type(CSF_Info_t), intent(in) :: csf_i
        integer, intent(in) :: src_orb
        integer, allocatable :: valid_orbs(:)

        integer(n_int) :: &
            ilut_doubly_occupied(0:this%L_spat_bits), ilut_valid(0:this%L_spat_bits)
        integer :: i_sg

        if (this%use_lookup) then
            i_sg = this%indexer%lookup_prop_vec_idx(csf_I)
        else
            i_sg = this%indexer%idx_csf(csf_i)
        end if

        block
            integer, allocatable :: doubly_occupied(:)
            integer :: i
            doubly_occupied = pack([(i, i=1, size(csf_i%stepvector))], csf_i%stepvector == 3)
            ilut_doubly_occupied = 0_n_int
            do i = 1, size(doubly_occupied)
                set_orb(ilut_doubly_occupied, doubly_occupied(i))
            end do
        endblock

        ilut_valid = iand(this%allowed_holes(:, src_orb, i_sg), not(ilut_doubly_occupied))
        allocate(valid_orbs(sum(popcnt(ilut_valid))))
        if (size(valid_orbs) /= 0) then
            call decode_bit_det(valid_orbs, ilut_valid)
        end if
    end function get_valid_orbs


    function calc_pgen_uniform_singles(this, nI, ilutI, csf_i, excitInfo) result(pgen)
        debug_function_name("calc_pgen_uniform_singles")
        class(PropVec_UniformExcGenerator_t), intent(in) :: this
        integer, intent(in) :: nI(nel)
        integer(n_int), intent(in) :: ilutI(0:GugaBits%len_tot)
        type(CSF_Info_t), intent(in) :: csf_i
        type(ExcitationInformation_t), intent(in) :: excitInfo
        real(dp) :: pgen

        integer :: nOcc, nUnocc
        integer, allocatable :: valid_orbs(:)

#ifdef WARNING_WORKAROUND_
        associate(nI => nI); end associate
        associate(ilutI => ilutI); end associate
#endif

        ASSERT(1 <= excitInfo%i .and. excitInfo%i <= nSpatOrbs)
        ASSERT(1 <= excitInfo%j .and. excitInfo%j <= nSpatOrbs)

        if (excitInfo%i == excitInfo%j &
               .or. csf_i%stepvector(excitInfo%i) == 3 &
               .or. csf_i%stepvector(excitInfo%j) == 0) then
            pgen = 0.0_dp
        else
            nOcc = count(csf_i%Occ_int /= 0)
            valid_orbs = this%get_valid_orbs(csf_i, excitInfo%j)
            nUnocc = size(valid_orbs)
            pgen = 1.0_dp / real(nOcc * nUnocc, dp)
        end if
    end function calc_pgen_uniform_singles


    elemental function alg_from_str(str) result(res)
        character(*), intent(in) :: str
        type(PropVec_PCHB_SinglesAlgorithm_t) :: res
        routine_name("alg_from_str")

        select case(to_upper(str))
        case(option_vals%algorithm%UNIF_UNIF%str)
            res = option_vals%algorithm%UNIF_UNIF
        case(option_vals%algorithm%FULL_FULL%str)
            res = option_vals%algorithm%FULL_FULL
        case(option_vals%algorithm%UNIF_FULL%str)
            res = option_vals%algorithm%UNIF_FULL
        case(option_vals%algorithm%UNIF_XNEW%str)
            res = option_vals%algorithm%UNIF_XNEW
        case default
            call stop_all(this_routine, 'Invalid input: '//str)
        end select
    end function

    subroutine init_PropVec_Weighted_t(this, weighting, indexer, use_lookup, create_lookup)
        class(PropVec_Weighted_t), intent(inout) :: this
        type(PropVec_PCHB_SinglesWeighting_t), intent(in) :: weighting
        class(AlsoGUGA_PropertyIndexer_t), intent(in) :: indexer
        logical, intent(in) :: use_lookup, create_lookup
        integer, allocatable :: prop_vecs(:, :)
        routine_name("init_PropVec_Weighted_t")

        integer :: i_sg, src, tgt, nBi
        real(dp), allocatable :: I_weights(:), A_weights(:)

        allocate(this%indexer, source=indexer)
        nBi = this%indexer%n_spin_orbs()
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        use constants, only: stderr
        if (.not. (nBi == 2 * nSpatOrbs)) then
            write(stderr, *) ""
            write(stderr, *) "Assertion nBi == 2 * nSpatOrbs"
            write(stderr, *) "failed in /home/cicdadmin/fkfest-gh-runners/altest-runner-7/_work/NECI_autoforward/NECI_autoforward/s&
                &rc/GUGA/guga_propvec_pchb_singles_main.fpp:435"
            call stop_all (this_routine, "Assert fail: nBi == 2 * nSpatOrbs")
        end if
    end block
#endif
        this%create_lookup = create_lookup
        if (create_lookup) then
            if (allocated(lookup_property_indexer)) then
                call stop_all(this_routine, 'Someone else is already managing the prop_vec lookup.')
            else
                write(stdout, *) 'PropVec singles is creating and managing the prop_vec lookup'
                allocate(lookup_property_indexer, source=this%indexer)
            end if
        end if
        this%use_lookup = use_lookup
        if (use_lookup) write(stdout, *) 'PropVec singles is using the prop_vec lookup'
        ! possible prop_vecs
        prop_vecs = this%indexer%get_prop_vecs()

        this%unoccupied = UnoccupiedGetter_t(nBasis)


        call set_get_weight_pointer(weighting, this%get_weight)

        call this%I_sampler%shared_alloc(size(prop_vecs, 2), nBi, 'PCHB')
        call this%A_sampler%shared_alloc([nSpatOrbs, size(prop_vecs, 2)], nBi, 'PCHB')

        allocate(I_weights(nBi), A_weights(nBi), source=0._dp)
        do i_sg = 1, size(prop_vecs, 2)
            I_weights(:) = 0._dp
            do src = 1, nSpatOrbs
                A_weights(:) = 0._dp
                do tgt = 1, nSpatOrbs
                block
                    type(Excite_1_t) :: exc
                    exc = Excite_1_t(2 * src, 2 * tgt)
                    if (src /= tgt &
                            .and. this%indexer%is_allowed(exc, prop_vecs(:, i_sg)) &
                            .and. symmetry_allowed(exc)) then
                        A_weights(2 * tgt) = abs(this%get_weight(src, tgt))
                        A_weights(2 * tgt - 1) = A_weights(2 * tgt)
                    end if
                end block
                end do
                I_weights(2 * src) = sum(A_weights)
                I_weights(2 * src - 1) = I_weights(2 * src)
                call this%A_sampler%setup_entry(src, i_sg, root, A_weights)
            end do
            call this%I_sampler%setup_entry(i_sg, root, I_weights)
        end do

    end subroutine

    subroutine finalize_PropVec_Weighted_t(this)
        class(PropVec_Weighted_t), intent(inout) :: this
        if (allocated(this%indexer)) then
            ! Yes I assume that either all or none are allocated
            call this%I_sampler%finalize()
            call this%A_sampler%finalize()
            nullify(this%get_weight)
            deallocate(this%indexer)
            if (this%create_lookup) deallocate(lookup_property_indexer)
        end if
    end subroutine


    subroutine pickOrbitals_FullFull(this, nI, ilut, csf_i, excitInfo, pgen)
        class(SinglesPropVec_FullFull_t), intent(in) :: this
        integer, intent(in) :: nI(nel)
        integer(n_int), intent(in) :: ilut(0:GugaBits%len_tot)
        type(CSF_Info_t), intent(in) :: csf_i
        type(ExcitationInformation_t), intent(out) :: excitInfo
        real(dp), intent(out) :: pgen
        routine_name("pickOrbitals_PropVec_Weighted_t")

        integer :: i_sg, src, tgt
        real(dp) :: p_src, p_tgt

        if (this%use_lookup) then
            i_sg = this%indexer%lookup_prop_vec_idx(csf_I)
        else
            i_sg = this%indexer%idx_nI(nI)
        end if

        select_particle: block
            real(dp) :: renorm_src
            integer :: spin_src, elec
            renorm_src = sum(this%I_sampler%get_prob(i_sg, nI))
            call this%I_sampler%constrained_sample(i_sg, nI, ilut(0 : nIfD), renorm_src, elec, spin_src, p_src)
            if (spin_src == 0) then
                call make_invalid()
                return
            end if
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        use constants, only: stderr
        if (.not. (nI(elec) == spin_src)) then
            write(stderr, *) ""
            write(stderr, *) "Assertion nI(elec) == spin_src"
            write(stderr, *) "failed in /home/cicdadmin/fkfest-gh-runners/altest-runner-7/_work/NECI_autoforward/NECI_autoforward/s&
                &rc/GUGA/guga_propvec_pchb_singles_main.fpp:524"
            call stop_all (this_routine, "Assert fail: nI(elec) == spin_src")
        end if
    end block
#endif
            src = gtId(spin_src)
            p_src = csf_I%occ_int(src) * p_src
        end block select_particle

        select_hole: block
            real(dp) :: renorm_tgt
            integer :: spin_tgt, dummy, unoccupied(nBasis - nEl)
            integer(n_int) :: ilut_unoccupied(0 : nIfD)
            call this%unoccupied%get(ilut(0 : nIfD), ilut_unoccupied, unoccupied)
            renorm_tgt = 1._dp - sum(this%A_sampler%get_prob(src, i_sg, nI))
            if (do_direct_calculation(renorm_tgt)) then
                renorm_tgt = sum(this%A_sampler%get_prob(src, i_sg, unoccupied))
            end if

            call this%A_sampler%constrained_sample(&
                 src, i_sg, unoccupied, ilut_unoccupied, renorm_tgt, dummy, spin_tgt, p_tgt)
            if (spin_tgt == 0) then
                call make_invalid()
                return
            end if
            tgt = gtId(spin_tgt)
            p_tgt = (2 - csf_I%occ_int(tgt)) * p_tgt
        end block select_hole
        pgen = p_src * p_tgt

        if (tgt < src) then
            excitInfo = assign_excitinfo_values_single( &
                        gen_type%R, tgt, src, tgt, src)
        else
            excitInfo = assign_excitinfo_values_single( &
                        gen_type%L, tgt, src, src, tgt)
        end if
        excitInfo%valid = checkCompatibility_single(csf_i, excitInfo)
        if (excitInfo%valid) then
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        use constants, only: stderr
        if (.not. (pgen .isclose. this%calc_pgen(nI, ilut, csf_i, excitInfo))) then
            write(stderr, *) ""
            write(stderr, *) "Assertion pgen .isclose. this%calc_pgen(nI, ilut, csf_i, excitInfo)"
            write(stderr, *) "failed in /home/cicdadmin/fkfest-gh-runners/altest-runner-7/_work/NECI_autoforward/NECI_autoforward/s&
                &rc/GUGA/guga_propvec_pchb_singles_main.fpp:559"
            call stop_all (this_routine, "Assert fail: pgen .isclose. this%calc_pgen(nI, ilut, csf_i, excitInfo)")
        end if
    end block
#endif
        else
            pgen = 0._dp
        end if

        contains

            subroutine make_invalid()
                pgen = 0.0_dp
                tgt = 0
                excitInfo%valid = .false.
            end subroutine
    end subroutine


    function calc_pgen_FullFull(this, nI, ilutI, csf_i, excitInfo) result(pgen)
        class(SinglesPropVec_FullFull_t), intent(in) :: this
        integer, intent(in) :: nI(nel)
        integer(n_int), intent(in) :: ilutI(0:GugaBits%len_tot)
        type(CSF_Info_t), intent(in) :: csf_i
        type(ExcitationInformation_t), intent(in) :: excitInfo
        real(dp) :: pgen
        routine_name("calc_pgen_FullFull")

        real(dp) :: p_src, p_tgt
        integer :: i_sg

#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        use constants, only: stderr
        if (.not. (1 <= excitInfo%i .and. excitInfo%i <= nSpatOrbs)) then
            write(stderr, *) ""
            write(stderr, *) "Assertion 1 <= excitInfo%i .and. excitInfo%i <= nSpatOrbs"
            write(stderr, *) "failed in /home/cicdadmin/fkfest-gh-runners/altest-runner-7/_work/NECI_autoforward/NECI_autoforward/s&
                &rc/GUGA/guga_propvec_pchb_singles_main.fpp:586"
            call stop_all (this_routine, "Assert fail: 1 <= excitInfo%i .and. excitInfo%i <= nSpatOrbs")
        end if
    end block
#endif
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        use constants, only: stderr
        if (.not. (1 <= excitInfo%j .and. excitInfo%j <= nSpatOrbs)) then
            write(stderr, *) ""
            write(stderr, *) "Assertion 1 <= excitInfo%j .and. excitInfo%j <= nSpatOrbs"
            write(stderr, *) "failed in /home/cicdadmin/fkfest-gh-runners/altest-runner-7/_work/NECI_autoforward/NECI_autoforward/s&
                &rc/GUGA/guga_propvec_pchb_singles_main.fpp:587"
            call stop_all (this_routine, "Assert fail: 1 <= excitInfo%j .and. excitInfo%j <= nSpatOrbs")
        end if
    end block
#endif

        if (excitInfo%i == excitInfo%j &
               .or. csf_i%stepvector(excitInfo%i) == 3 &
               .or. csf_i%stepvector(excitInfo%j) == 0) then
            pgen = 0.0_dp
            return
        end if

        if (this%use_lookup) then
            i_sg = this%indexer%lookup_prop_vec_idx(csf_I)
        else
            i_sg = this%indexer%idx_csf(csf_i)
        end if

        select_particle: block
            real(dp) :: renorm_src
            integer :: spin_src
            spin_src = 2 * excitInfo%j - 1 + ((csf_I%stepvector(excitInfo%j) + 1) .div. 3)
            renorm_src = sum(this%I_sampler%get_prob(i_sg, nI))
            p_src = csf_I%occ_int(excitInfo%j) &
                    * this%I_sampler%constrained_getProb(i_sg, nI, renorm_src, spin_src)
        end block select_particle

        select_hole: block
            real(dp) :: renorm_tgt
            integer :: spin_tgt, unoccupied(nBasis - nEl)
            integer(n_int) :: ilut_unoccupied(0 : nIfD)
            spin_tgt = 2 * excitInfo%i - ((csf_I%stepvector(excitInfo%i) + 1) .div. 3)

            call this%unoccupied%get(ilutI(0 : nIfD), ilut_unoccupied, unoccupied)
            renorm_tgt = 1._dp - sum(this%A_sampler%get_prob(excitInfo%j, i_sg, nI))
            if (do_direct_calculation(renorm_tgt)) then
                renorm_tgt = sum(this%A_sampler%get_prob(excitInfo%j, i_sg, unoccupied))
            end if


            p_tgt = (2 - csf_I%occ_int(excitInfo%i)) &
                    * this%A_sampler%constrained_getProb(excitInfo%j, i_sg, unoccupied, renorm_tgt, spin_tgt)
        end block select_hole
        pgen = p_src * p_tgt
    end function

    subroutine pickOrbitals_UnifFull(this, nI, ilut, csf_i, excitInfo, pgen)
        class(SinglesPropVec_UnifFull_t), intent(in) :: this
        integer, intent(in) :: nI(nel)
        integer(n_int), intent(in) :: ilut(0:GugaBits%len_tot)
        type(CSF_Info_t), intent(in) :: csf_i
        type(ExcitationInformation_t), intent(out) :: excitInfo
        real(dp), intent(out) :: pgen
        routine_name("pickOrbitals_PropVec_Weighted_t")

        integer :: i_sg, src, tgt
        real(dp) :: p_src, p_tgt

        if (this%use_lookup) then
            i_sg = this%indexer%lookup_prop_vec_idx(csf_I)
        else
            i_sg = this%indexer%idx_nI(nI)
        end if

        select_particle: block
            src = gtID(nI(int(genrand_real2_dSFMT() * nEl) + 1))
            p_src = csf_I%occ_real(src) / real(nEl, dp)
        end block select_particle

        select_hole: block
            real(dp) :: renorm_tgt
            integer :: spin_tgt, dummy, unoccupied(nBasis - nEl)
            integer(n_int) :: ilut_unoccupied(0 : nIfD)
            call this%unoccupied%get(ilut(0 : nIfD), ilut_unoccupied, unoccupied)
            renorm_tgt = 1._dp - sum(this%A_sampler%get_prob(src, i_sg, nI))
            if (do_direct_calculation(renorm_tgt)) then
                renorm_tgt = sum(this%A_sampler%get_prob(src, i_sg, unoccupied))
            end if

            call this%A_sampler%constrained_sample(&
                 src, i_sg, unoccupied, ilut_unoccupied, renorm_tgt, dummy, spin_tgt, p_tgt)
            if (spin_tgt == 0) then
                call make_invalid()
                return
            end if
            tgt = gtId(spin_tgt)
            p_tgt = (2 - csf_I%occ_int(tgt)) * p_tgt
        end block select_hole
        pgen = p_src * p_tgt

        if (tgt < src) then
            excitInfo = assign_excitinfo_values_single( &
                        gen_type%R, tgt, src, tgt, src)
        else
            excitInfo = assign_excitinfo_values_single( &
                        gen_type%L, tgt, src, src, tgt)
        end if
        excitInfo%valid = checkCompatibility_single(csf_i, excitInfo)
        if (excitInfo%valid) then
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        use constants, only: stderr
        if (.not. (pgen .isclose. this%calc_pgen(nI, ilut, csf_i, excitInfo))) then
            write(stderr, *) ""
            write(stderr, *) "Assertion pgen .isclose. this%calc_pgen(nI, ilut, csf_i, excitInfo)"
            write(stderr, *) "failed in /home/cicdadmin/fkfest-gh-runners/altest-runner-7/_work/NECI_autoforward/NECI_autoforward/s&
                &rc/GUGA/guga_propvec_pchb_singles_main.fpp:683"
            call stop_all (this_routine, "Assert fail: pgen .isclose. this%calc_pgen(nI, ilut, csf_i, excitInfo)")
        end if
    end block
#endif
        else
            pgen = 0._dp
        end if

        contains

            subroutine make_invalid()
                pgen = 0.0_dp
                tgt = 0
                excitInfo%valid = .false.
            end subroutine
    end subroutine


    function calc_pgen_UnifFull(this, nI, ilutI, csf_i, excitInfo) result(pgen)
        class(SinglesPropVec_UnifFull_t), intent(in) :: this
        integer, intent(in) :: nI(nel)
        integer(n_int), intent(in) :: ilutI(0:GugaBits%len_tot)
        type(CSF_Info_t), intent(in) :: csf_i
        type(ExcitationInformation_t), intent(in) :: excitInfo
        real(dp) :: pgen
        routine_name("calc_pgen_UnifFull")

        real(dp) :: p_src, p_tgt
        integer :: i_sg

#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        use constants, only: stderr
        if (.not. (1 <= excitInfo%i .and. excitInfo%i <= nSpatOrbs)) then
            write(stderr, *) ""
            write(stderr, *) "Assertion 1 <= excitInfo%i .and. excitInfo%i <= nSpatOrbs"
            write(stderr, *) "failed in /home/cicdadmin/fkfest-gh-runners/altest-runner-7/_work/NECI_autoforward/NECI_autoforward/s&
                &rc/GUGA/guga_propvec_pchb_singles_main.fpp:710"
            call stop_all (this_routine, "Assert fail: 1 <= excitInfo%i .and. excitInfo%i <= nSpatOrbs")
        end if
    end block
#endif
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        use constants, only: stderr
        if (.not. (1 <= excitInfo%j .and. excitInfo%j <= nSpatOrbs)) then
            write(stderr, *) ""
            write(stderr, *) "Assertion 1 <= excitInfo%j .and. excitInfo%j <= nSpatOrbs"
            write(stderr, *) "failed in /home/cicdadmin/fkfest-gh-runners/altest-runner-7/_work/NECI_autoforward/NECI_autoforward/s&
                &rc/GUGA/guga_propvec_pchb_singles_main.fpp:711"
            call stop_all (this_routine, "Assert fail: 1 <= excitInfo%j .and. excitInfo%j <= nSpatOrbs")
        end if
    end block
#endif

        if (excitInfo%i == excitInfo%j &
               .or. csf_i%stepvector(excitInfo%i) == 3 &
               .or. csf_i%stepvector(excitInfo%j) == 0) then
            pgen = 0.0_dp
            return
        end if

        if (this%use_lookup) then
            i_sg = this%indexer%lookup_prop_vec_idx(csf_I)
        else
            i_sg = this%indexer%idx_csf(csf_i)
        end if

        select_particle: block
            p_src = csf_I%occ_real(excitInfo%j) / real(nEl, dp)
        end block select_particle

        select_hole: block
            real(dp) :: renorm_tgt
            integer :: spin_tgt, unoccupied(nBasis - nEl)
            integer(n_int) :: ilut_unoccupied(0 : nIfD)
            spin_tgt = 2 * excitInfo%i - ((csf_I%stepvector(excitInfo%i) + 1) .div. 3)

            call this%unoccupied%get(ilutI(0 : nIfD), ilut_unoccupied, unoccupied)
            renorm_tgt = 1._dp - sum(this%A_sampler%get_prob(excitInfo%j, i_sg, nI))
            if (do_direct_calculation(renorm_tgt)) then
                renorm_tgt = sum(this%A_sampler%get_prob(excitInfo%j, i_sg, unoccupied))
            end if


            p_tgt = (2 - csf_I%occ_int(excitInfo%i)) &
                    * this%A_sampler%constrained_getProb(excitInfo%j, i_sg, unoccupied, renorm_tgt, spin_tgt)
        end block select_hole
        pgen = p_src * p_tgt

    end function



    subroutine init_UnifXnew_t(this, weighting, indexer, use_lookup, create_lookup, scale_down)
        class(UnifXnew_t), intent(inout) :: this
        type(PropVec_PCHB_SinglesWeighting_t), intent(in) :: weighting
        class(AlsoGUGA_PropertyIndexer_t), intent(in) :: indexer
        logical, intent(in) :: use_lookup, create_lookup
        real(dp), intent(in) :: scale_down
            !! Scale down the "exchange" single excitations
        integer, allocatable :: prop_vecs(:, :)
        routine_name("init_UnifXnew_t")

        integer :: i_sg, src, tgt, nBi
        real(dp), allocatable :: spat_w_A(:), A_weights(:)

        allocate(this%indexer, source=indexer)
        nBi = this%indexer%n_spin_orbs()
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        use constants, only: stderr
        if (.not. (nBi == 2 * nSpatOrbs)) then
            write(stderr, *) ""
            write(stderr, *) "Assertion nBi == 2 * nSpatOrbs"
            write(stderr, *) "failed in /home/cicdadmin/fkfest-gh-runners/altest-runner-7/_work/NECI_autoforward/NECI_autoforward/s&
                &rc/GUGA/guga_propvec_pchb_singles_main.fpp:767"
            call stop_all (this_routine, "Assert fail: nBi == 2 * nSpatOrbs")
        end if
    end block
#endif
        this%create_lookup = create_lookup
        if (create_lookup) then
            if (allocated(lookup_property_indexer)) then
                call stop_all(this_routine, 'Someone else is already managing the prop_vec lookup.')
            else
                lookup_property_indexer = this%indexer
            end if
        end if
        this%use_lookup = use_lookup

        prop_vecs = this%indexer%get_prop_vecs()

        this%unoccupied = UnoccupiedGetter_t(nBasis)

        call set_get_weight_pointer(weighting, this%get_weight)

        call this%A_sampler%shared_alloc([nSpatOrbs, 3, size(prop_vecs, 2)], nBi, 'PCHB')

        allocate(spat_w_A(nSpatOrbs), A_weights(nBi), source=0._dp)
        do i_sg = 1, size(prop_vecs, 2)
            do src = 1, nSpatOrbs
                A_weights(:) = 0._dp
                spat_w_A(:) = 0._dp
                do tgt = 1, nSpatOrbs
                block
                    type(Excite_1_t) :: exc
                    exc = Excite_1_t(2 * src, 2 * tgt)
                    if (src /= tgt &
                            .and. this%indexer%is_allowed(exc, prop_vecs(:, i_sg)) &
                            .and. symmetry_allowed(exc)) then
                        spat_w_A(tgt) = abs(this%get_weight(src, tgt))
                    end if
                end block
                end do

                start_from_3: block
                    A_weights(1::2) = spat_w_A(:)
                    A_weights(2::2) = spat_w_A(:)
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        use constants, only: stderr
        if (.not. (near_zero(A_weights(2 * src - 1)) .and. near_zero(A_weights(2 * src)))) then
            write(stderr, *) ""
            write(stderr, *) "Assertion near_zero(A_weights(2 * src - 1)) .and. near_zero(A_weights(2 * src))"
            write(stderr, *) "failed in /home/cicdadmin/fkfest-gh-runners/altest-runner-7/_work/NECI_autoforward/NECI_autoforward/s&
                &rc/GUGA/guga_propvec_pchb_singles_main.fpp:806"
            call stop_all (this_routine, "Assert fail: near_zero(A_weights(2 * src - 1)) .and. near_zero(A_weights(2 * src))")
        end if
    end block
#endif
                    call this%A_sampler%setup_entry(src, 3, i_sg, root, A_weights)
                end block start_from_3

                start_from_1: block
                    A_weights(1::2) = spat_w_A(:)
                    A_weights(2::2) = spat_w_A(:) * scale_down
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        use constants, only: stderr
        if (.not. (near_zero(A_weights(2 * src - 1)) .and. near_zero(A_weights(2 * src)))) then
            write(stderr, *) ""
            write(stderr, *) "Assertion near_zero(A_weights(2 * src - 1)) .and. near_zero(A_weights(2 * src))"
            write(stderr, *) "failed in /home/cicdadmin/fkfest-gh-runners/altest-runner-7/_work/NECI_autoforward/NECI_autoforward/s&
                &rc/GUGA/guga_propvec_pchb_singles_main.fpp:813"
            call stop_all (this_routine, "Assert fail: near_zero(A_weights(2 * src - 1)) .and. near_zero(A_weights(2 * src))")
        end if
    end block
#endif
                    if (src > 1) A_weights(2 * (src - 1)) = 0._dp
                    if (src < nSpatOrbs) A_weights(2 * (src + 1)) = 0._dp
                    call this%A_sampler%setup_entry(src, 1, i_sg, root, A_weights)
                end block start_from_1

                start_from_2: block
                    A_weights(1::2) = spat_w_A(:) * scale_down
                    A_weights(2::2) = spat_w_A(:)
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        use constants, only: stderr
        if (.not. (near_zero(A_weights(2 * src - 1)) .and. near_zero(A_weights(2 * src)))) then
            write(stderr, *) ""
            write(stderr, *) "Assertion near_zero(A_weights(2 * src - 1)) .and. near_zero(A_weights(2 * src))"
            write(stderr, *) "failed in /home/cicdadmin/fkfest-gh-runners/altest-runner-7/_work/NECI_autoforward/NECI_autoforward/s&
                &rc/GUGA/guga_propvec_pchb_singles_main.fpp:822"
            call stop_all (this_routine, "Assert fail: near_zero(A_weights(2 * src - 1)) .and. near_zero(A_weights(2 * src))")
        end if
    end block
#endif
                    if (src > 1) A_weights(2 * (src - 1) - 1) = 0._dp
                    if (src < nSpatOrbs) A_weights(2 * (src + 1) - 1) = 0._dp
                    call this%A_sampler%setup_entry(src, 2, i_sg, root, A_weights)
                end block start_from_2
            end do
        end do
    end subroutine

    subroutine finalize_UnifXnew_t(this)
        class(UnifXnew_t), intent(inout) :: this
        call this%A_sampler%finalize()
        nullify(this%get_weight)
        deallocate(this%indexer)
        if (this%create_lookup) deallocate(lookup_property_indexer)
    end subroutine

    subroutine pickOrbitals_UnifXnew_t(this, nI, ilut, csf_i, excitInfo, pgen)
        class(UnifXnew_t), intent(in) :: this
        integer, intent(in) :: nI(nel)
        integer(n_int), intent(in) :: ilut(0:GugaBits%len_tot)
        type(CSF_Info_t), intent(in) :: csf_i
        type(ExcitationInformation_t), intent(out) :: excitInfo
        real(dp), intent(out) :: pgen
        routine_name("pickOrbitals_PropVec_Weighted_t")

        integer :: i_sg, src, tgt
        real(dp) :: p_src, p_tgt
        real(dp) :: renorm_tgt
        integer :: spin_tgt, dummy, unoccupied(nBasis - nEl)
        integer(n_int) :: ilut_unoccupied(0 : nIfD)

        if (this%use_lookup) then
            i_sg = this%indexer%lookup_prop_vec_idx(csf_I)
        else
            i_sg = this%indexer%idx_nI(nI)
        end if

        select_particle: block
            src = gtID(nI(int(genrand_real2_dSFMT() * nEl) + 1))
            p_src = csf_I%occ_real(src) / real(nEl, dp)
        end block select_particle

        select_hole: block

            call this%unoccupied%get(ilut(0 : nIfD), ilut_unoccupied, unoccupied)

            renorm_tgt = 1._dp - sum(this%A_sampler%get_prob(src, csf_i%stepvector(src), i_sg, nI))
            if (do_direct_calculation(renorm_tgt)) then
                renorm_tgt = sum(this%A_sampler%get_prob(src, csf_i%stepvector(src), i_sg, unoccupied))
            end if

            call this%A_sampler%constrained_sample(&
                 src, csf_i%stepvector(src), i_sg, unoccupied, ilut_unoccupied, &
                 renorm_tgt, dummy, spin_tgt, p_tgt)
            if (spin_tgt == 0) then
                call make_invalid()
                return
            end if
            tgt = gtId(spin_tgt)
            if (csf_I%occ_int(tgt) == 0) then
                if (mod(spin_tgt, 2) == 1) then
                    p_tgt = p_tgt &
                            + this%A_sampler%constrained_getProb(&
                                    src, csf_i%stepvector(src), i_sg, &
                                    unoccupied, renorm_tgt, 2 * tgt)
                else
                    p_tgt = p_tgt &
                            + this%A_sampler%constrained_getProb(&
                                    src, csf_i%stepvector(src), i_sg, &
                                    unoccupied, renorm_tgt, 2 * tgt - 1)
                end if
            end if
        end block select_hole
        pgen = p_src * p_tgt

        if (tgt < src) then
            excitInfo = assign_excitinfo_values_single( &
                        gen_type%R, tgt, src, tgt, src)
        else
            excitInfo = assign_excitinfo_values_single( &
                        gen_type%L, tgt, src, src, tgt)
        end if
        excitInfo%valid = checkCompatibility_single(csf_i, excitInfo)
        if (excitInfo%valid) then
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        use constants, only: stderr
        if (.not. (pgen .isclose. this%calc_pgen(nI, ilut, csf_i, excitInfo))) then
            write(stderr, *) ""
            write(stderr, *) "Assertion pgen .isclose. this%calc_pgen(nI, ilut, csf_i, excitInfo)"
            write(stderr, *) "failed in /home/cicdadmin/fkfest-gh-runners/altest-runner-7/_work/NECI_autoforward/NECI_autoforward/s&
                &rc/GUGA/guga_propvec_pchb_singles_main.fpp:907"
            call stop_all (this_routine, "Assert fail: pgen .isclose. this%calc_pgen(nI, ilut, csf_i, excitInfo)")
        end if
    end block
#endif
        else
            pgen = 0._dp
        end if

        contains

            subroutine make_invalid()
                pgen = 0.0_dp
                tgt = 0
                excitInfo%valid = .false.
            end subroutine
    end subroutine


    function calc_pgen_UnifXnew_t(this, nI, ilutI, csf_i, excitInfo) result(pgen)
        class(UnifXnew_t), intent(in) :: this
        integer, intent(in) :: nI(nel)
        integer(n_int), intent(in) :: ilutI(0:GugaBits%len_tot)
        type(CSF_Info_t), intent(in) :: csf_i
        type(ExcitationInformation_t), intent(in) :: excitInfo
        real(dp) :: pgen
        routine_name("calc_pgen_UnifXnew_t")

        real(dp) :: p_src, p_tgt
        integer :: i_sg

        associate(tgt => excitInfo%i, src => excitInfo%j)

#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        use constants, only: stderr
        if (.not. (1 <= src .and. src <= nSpatOrbs)) then
            write(stderr, *) ""
            write(stderr, *) "Assertion 1 <= src .and. src <= nSpatOrbs"
            write(stderr, *) "failed in /home/cicdadmin/fkfest-gh-runners/altest-runner-7/_work/NECI_autoforward/NECI_autoforward/s&
                &rc/GUGA/guga_propvec_pchb_singles_main.fpp:936"
            call stop_all (this_routine, "Assert fail: 1 <= src .and. src <= nSpatOrbs")
        end if
    end block
#endif
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        use constants, only: stderr
        if (.not. (1 <= tgt .and. tgt <= nSpatOrbs)) then
            write(stderr, *) ""
            write(stderr, *) "Assertion 1 <= tgt .and. tgt <= nSpatOrbs"
            write(stderr, *) "failed in /home/cicdadmin/fkfest-gh-runners/altest-runner-7/_work/NECI_autoforward/NECI_autoforward/s&
                &rc/GUGA/guga_propvec_pchb_singles_main.fpp:937"
            call stop_all (this_routine, "Assert fail: 1 <= tgt .and. tgt <= nSpatOrbs")
        end if
    end block
#endif


        if (src == tgt &
               .or. csf_i%occ_int(tgt) == 2 &
               .or. csf_i%occ_int(src) == 0) then
            pgen = 0.0_dp
            return
        end if

        if (this%use_lookup) then
            i_sg = this%indexer%lookup_prop_vec_idx(csf_I)
        else
            i_sg = this%indexer%idx_csf(csf_i)
        end if

        select_particle: block
            p_src = csf_I%occ_real(src) / real(nEl, dp)
        end block select_particle

        select_hole: block
            real(dp) :: renorm_tgt
            integer :: unoccupied(nBasis - nEl)
            integer(n_int) :: ilut_unoccupied(0 : nIfD)

            call this%unoccupied%get(ilutI(0 : nIfD), ilut_unoccupied, unoccupied)

            renorm_tgt = 1._dp - sum(this%A_sampler%get_prob(src, csf_i%stepvector(src), i_sg, nI))
            if (do_direct_calculation(renorm_tgt)) then
                renorm_tgt = sum(this%A_sampler%get_prob(src, csf_i%stepvector(src), i_sg, unoccupied))
            end if

            if (csf_i%occ_int(tgt) == 0) then
                p_tgt = this%A_sampler%constrained_getProb(src, csf_i%stepvector(src), i_sg, unoccupied, renorm_tgt, 2 * tgt - 1) &
                        + this%A_sampler%constrained_getProb(src, csf_i%stepvector(src), i_sg, unoccupied, renorm_tgt, 2 * tgt)
            else if (csf_i%stepvector(tgt) == 1) then
                p_tgt = this%A_sampler%constrained_getProb(src, csf_i%stepvector(src), i_sg, unoccupied, renorm_tgt, 2 * tgt)
            else if (csf_i%stepvector(tgt) == 2) then
                p_tgt = this%A_sampler%constrained_getProb(src, csf_i%stepvector(src), i_sg, unoccupied, renorm_tgt, 2 * tgt - 1)
            else
                call stop_all(this_routine, "should not be here.")
            end if

        end block select_hole
        pgen = p_src * p_tgt
        end associate
    end function

    pure function to_str_algorithm(algorithm) result(res)
        class(PropVec_PCHB_SinglesAlgorithm_t), intent(in) :: algorithm
        character(9) :: res
        routine_name("to_str_algorithm")
        if (algorithm == option_vals%algorithm%UNIF_UNIF) then
            res = option_vals%algorithm%UNIF_UNIF%str
        else if (algorithm == option_vals%algorithm%UNIF_FULL) then
            res = option_vals%algorithm%UNIF_FULL%str
        else if (algorithm == option_vals%algorithm%FULL_FULL) then
            res = option_vals%algorithm%FULL_FULL%str
        else if (algorithm == option_vals%algorithm%UNIF_XNEW) then
            res = option_vals%algorithm%UNIF_XNEW%str
        else
            call stop_all(this_routine, "Should not be here.")
        end if
    end function

    pure function to_str_singles(singles_options) result(res)
        class(PropVec_PCHB_SinglesOptions_t), intent(in) :: singles_options
        character(:), allocatable :: res
        res = singles_options%algorithm%to_str() // ' ' // singles_options%weighting%to_str()
    end function

end module
