GAS_doubles_PCHB_compute_samplers Subroutine

    subroutine GAS_doubles_PCHB_compute_samplers(this, PCHB_particle_selection)
        !! computes and stores values for the alias sampling table
        class(GAS_PCHB_DoublesSpinorbFullyWeightedExcGenerator_t), intent(inout) :: this
        type(PCHB_ParticleSelection_t), intent(in) :: PCHB_particle_selection
        integer :: I, J, IJ, IJ_max, A, B ! Uppercase because they are indexing spin orbitals
        integer(int64) :: memCost
        real(dp), allocatable :: w_A(:), w_B(:), IJ_weights(:, :, :)
        integer, allocatable :: supergroups(:, :)
        integer :: i_sg
        ! possible supergroups
        supergroups = this%indexer%get_supergroups()

        ! number of possible source orbital pairs
        IJ_max = fuseIndex(nBasis, nBasis)

        !> n_supergroup * number_of_fused_indices * 3 * (bytes_per_sampler)
        memCost = size(supergroups, 2, kind=int64) &
                    * (int(IJ_max, int64) * (int(nBasis, int64) * 3_int64 * 8_int64) & ! p(A | IJ)
                        + int(IJ_max, int64) * (int(nBasis, int64)**2 * 3_int64 * 8_int64) & ! P (B | IJ A)
                    )

        write(stdout, *) "Excitation generator requires", real(memCost, dp) / 2.0_dp**30, "GB of memory"
        write(stdout, *) "The number of supergroups is", size(supergroups, 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%A_sampler%shared_alloc([IJ_max, size(supergroups, 2)], nBasis, 'A_sampler_PCHB_spinorb_FullyWeighted')
        call this%B_sampler%shared_alloc([nBasis, IJ_max, size(supergroups, 2)], nBasis, 'A_sampler_PCHB_spinorb_FullyWeighted')

        if (iProcIndex_intra == root) then
            allocate(w_A(nBasis), w_B(nBasis), source=0._dp)
            allocate(IJ_weights(nBasis, nBasis, size(supergroups, 2)), source=0._dp)
        else
            allocate(w_A(0), w_B(0), IJ_weights(0, 0, 0))
        end if

        ! initialize the three samplers
        do i_sg = 1, size(supergroups, 2)
            if (mod(i_sg, 100) == 0) write(stdout, *) 'Still generating the samplers'
            first_particle: do I = 1, nBasis
                second_particle: do J = I + 1, nBasis
                    IJ = fuseIndex(I, J)
                    w_A(:) = 0.0_dp
                    first_hole: do A = 1, nBasis
                        if (any(A == [I, J])) cycle
                        w_B(:) = 0.0_dp
                        second_hole: do B = 1, nBasis
                            if (any(B == [I, J, A])) cycle
                            associate(exc => merge(Excite_2_t(I, A, J, B), Excite_2_t(I, B, J, A), A < B))
                            if (iProcIndex_intra == root) then
                                if (this%GAS_spec%is_allowed(exc, supergroups(:, i_sg)) &
                                        .and. spin_allowed(exc)) then
                                    w_B(B) = get_PCHB_weight(exc)
                                end if
                            end if
                            end associate
                        end do second_hole
                        call this%B_sampler%setup_entry(A, IJ, i_sg, root, w_B(:))
                        if (iProcIndex_intra == root) then
                            w_A(A) = sum(w_B)
                        end if
                    end do first_hole
                    call this%A_sampler%setup_entry(IJ, i_sg, root, w_A(:))

                    if (iProcIndex_intra == root) then
                        IJ_weights(I, J, i_sg) = sum(w_A)
                        IJ_weights(J, I, i_sg) = IJ_weights(I, J, i_sg)
                    end if
                end do second_particle
            end do first_particle
        end do

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

    end subroutine GAS_doubles_PCHB_compute_samplers