PropVec_doubles_PCHB_compute_samplers Subroutine

    subroutine PropVec_doubles_PCHB_compute_samplers(this, nBI, PCHB_particle_selection)
        !! computes and stores values for the alias sampling table
        class(PropVec_PCHB_DoublesSpatOrbFastWeightExcGen_t), intent(inout) :: this
        integer, intent(in) :: nBI
        type(PCHB_ParticleSelection_t), intent(in) :: PCHB_particle_selection
        integer :: i, j, ij, ijMax
        integer :: a, b, ab, abMax
        integer :: ex(2, 2)
        integer(int64) :: memCost
        real(dp), allocatable :: w(:), pNoExch(:), IJ_weights(:, :, :)
        integer, allocatable :: prop_vecs(:, :)
        integer :: i_sg, i_exch
        ! possible prop_vecs
        prop_vecs = this%indexer%get_prop_vecs()

        ! number of possible source orbital pairs
        ijMax = fuse_symm_idx(nBI, nBI)
        abMax = ijMax
        ! allocate the bias for picking an exchange excitation
        allocate(this%pExch(ijMax, size(prop_vecs, 2)), source=0.0_dp)
        ! temporary storage for the unnormalized prob of not picking an exchange excitation

        !> n_prop_vec * number_of_fused_indices * 3 * (bytes_per_sampler)
        memCost = size(prop_vecs, 2, kind=int64) &
                    * int(ijMax, int64) &
                    * 3_int64 &
                    * (int(abMax, int64) * 3_int64 * 8_int64)

        write(stdout, *) "Excitation generator requires", real(memCost, dp) / 2.0_dp**30, "GB of memory"
        write(stdout, *) "The number of prop_vecs is", size(prop_vecs, 2)
        write(stdout, *) "Generating samplers for PCHB excitation generator"
        write(stdout, *) "Depending on the number of prop_vecs this can take up to 10min."
        call this%pchb_samplers%shared_alloc([ijMax, 3, size(prop_vecs, 2)], abMax, 'PCHB_RHF')

        ! One could allocate only on the intra-node-root here, if memory
        ! at initialization ever becomes an issue.
        ! Look at `gasci_pchb_doubles_spin_fulllyweighted.fpp` for inspiration.
        allocate(w(abMax))
        allocate(IJ_weights(nBI * 2, nBI * 2, size(prop_vecs, 2)), source=0._dp)
        over_prop_vec: do i_sg = 1, size(prop_vecs, 2)
            if (mod(i_sg, 100) == 0) write(stdout, *) 'Still generating the samplers'
            pNoExch = 1.0_dp - this%pExch(:, i_sg)
            over_spin_type: do i_exch = 1, 3
                particle_1: do i = 1, nBI
                    ex(1, 1) = to_spin_orb(i, is_alpha=.true.)
                    particle_2: do j = i, nBi
                        if (i_exch == SAME_SPIN .and. i == j) cycle
                        ij = fuse_symm_idx(i, j)
                        w(:) = 0.0_dp
                        ex(1, 2) = to_spin_orb(j, i_exch == SAME_SPIN)
                        hole_1: do a = 1, nBI
                            ex(2, 1) = to_spin_orb(a, any(i_exch == [SAME_SPIN, OPP_SPIN_NO_EXCH]))
                            if (any(ex(2, 1) == ex(1, :))) cycle
                            hole_2: do b = a, nBi
                                if (i_exch == OPP_SPIN_EXCH .and. a == b) cycle
                                ab = fuse_symm_idx(a, b)
                                ex(2, 2) = to_spin_orb(b, any(i_exch == [SAME_SPIN, OPP_SPIN_EXCH]))
                                if (any(ex(2, 2) == ex(1, :)) .or. ex(2, 2) == ex(2, 1)) cycle

                                associate(exc => canonicalize(Excite_2_t(ex)))
                                if (this%indexer%is_allowed(exc, prop_vecs(:, i_sg))) then
                                    w(ab) = get_PCHB_weight(exc)
                                end if
                                end associate
                            end do hole_2
                        end do hole_1

                        call this%pchb_samplers%setup_entry(ij, i_exch, i_sg, root, w)
                        if (i_exch == OPP_SPIN_EXCH) this%pExch(ij, i_sg) = sum(w)
                        if (i_exch == OPP_SPIN_NO_EXCH) pNoExch(ij) = sum(w)

                        associate(I => ex(1, 1), J => ex(1, 2))
                            IJ_weights(I, J, i_sg) = IJ_weights(I, J, i_sg) + sum(w)
                            IJ_weights(J, I, i_sg) = IJ_weights(I, J, i_sg)
                        end associate
                        if (i /= j) then
                            ! sum over alpha and beta of the same orbital
                            if (i_exch == SAME_SPIN) then
                                associate(I => ex(1, 1) - 1, J => ex(1, 2) - 1)
                                    IJ_weights(I, J, i_sg) = IJ_weights(I, J, i_sg) + sum(w)
                                    IJ_weights(J, I, i_sg) = IJ_weights(I, J, i_sg)
                                end associate
                            else
                                associate(I => ex(1, 1) - 1, J => ex(1, 2) + 1)
                                    IJ_weights(I, J, i_sg) = IJ_weights(I, J, i_sg) + sum(w)
                                    IJ_weights(J, I, i_sg) = IJ_weights(I, J, i_sg)
                                end associate
                            end if
                        end if

                    end do particle_2
                end do particle_1
            end do over_spin_type
            ! normalize the exchange bias (where normalizable)
            where (near_zero(this%pExch(:, i_sg) + pNoExch))
                this%pExch(:, i_sg) = 0._dp
            else where
                this%pExch(:, i_sg) = this%pExch(:, i_sg) / (this%pExch(:, i_sg) + pNoExch)
            end where
        end do over_prop_vec


        call allocate_and_init(PCHB_particle_selection, this%indexer, &
            IJ_weights, root, this%use_lookup, this%particle_selector)

    contains
        elemental function to_spin_orb(orb, is_alpha) result(sorb)
            ! map spatial orbital to the spin orbital matching the current samplerIndex
            ! Input: orb - spatial orbital to be mapped
            ! Output: sorb - corresponding spin orbital
            integer, intent(in) :: orb
            logical, intent(in) :: is_alpha
            integer :: sorb

            sorb = merge(2 * orb, 2 * orb - 1, is_alpha)
        end function to_spin_orb
    end subroutine PropVec_doubles_PCHB_compute_samplers