PropVec_doubles_PCHB_spinorb_compute_samplers Subroutine

    subroutine PropVec_doubles_PCHB_spinorb_compute_samplers(this, nBI, PCHB_particle_selection)
        !! computes and stores values for the alias (spin-independent) sampling table
        class(PropVec_PCHB_DoublesSpinOrbFastWeightExcGen_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(int64) :: memCost
            !! n_prop_vec * num_fused_indices * bytes_per_sampler
        real(dp), allocatable :: w(:), IJ_weights(:, :, :)
        integer, allocatable :: prop_vecs(:, :)
        integer :: i_sg
        ! possible prop_vecs
        prop_vecs = this%indexer%get_prop_vecs()

        ! number of possible spin orbital pairs
        ijMax = fuse_symm_idx(nBI, nBI)
        abMax = ijMax

        memCost = size(prop_vecs, 2, kind=int64) &
                  * int(ijMax, int64) &
                  * int(abMax, int64) * 8_int64

        write(stdout, *) "Excitation generator requires", real(memCost, dp) / 2.0_dp**30, "GB of memory"
        write(stdout, *) "The number of supergroups is", size(prop_vecs, 2)
        write(stdout, *) "Generating samplers for PCHB excitation generator"
        write(stdout, *) "Depending on the number of supergroups this can take up to 10min."

        call this%AB_sampler%shared_alloc([ijMax, size(prop_vecs, 2)], abMax, 'AB_PCHB_spinorb_sampler')

        ! 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), IJ_weights(nBI, nBI, size(prop_vecs, 2)), source=0._dp)

        prop_vec: do i_sg = 1, size(prop_vecs, 2)
            if (mod(i_sg, 100) == 0) write(stdout, *) 'Still generating the samplers'
            first_particle: do I = 1, nBI
                ! A,B are ordered, so we can assume I < J
                second_particle: do J = I + 1, nBI
                    w(:) = 0._dp
                    ! for each (I,J), get all matrix elements <IJ|H|AB> and use them as
                    ! weights to prepare the sampler
                    first_hole: do A = 1, nBI
                        if (any(A == [I, J])) cycle
                        second_hole: do B = 1, nBI
                            if (any(B == [I, J, A])) cycle
                            AB = fuse_symm_idx(A, B)
                            associate(exc => merge(Excite_2_t(I, A, J, B), Excite_2_t(I, B, J, A), A < B))
                                if (this%indexer%is_allowed(exc, prop_vecs(:, i_sg)) .and. spin_allowed(exc)) then
                                    w(AB) = get_PCHB_weight(exc)
                                end if
                            end associate
                        end do second_hole
                    end do first_hole
                    IJ = fuse_symm_idx(I, J)
                    call this%AB_sampler%setup_entry(IJ, i_sg, root, w)

                    IJ_weights(I, J, i_sg) = sum(w)
                    IJ_weights(J, I, i_sg) = IJ_weights(I, J, i_sg)
                end do second_particle
            end do first_particle
        end do prop_vec

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

    end subroutine PropVec_doubles_PCHB_spinorb_compute_samplers