init_pcpp_singles_excitgen Subroutine

    subroutine init_pcpp_singles_excitgen()
        ! create the probability distributions for drawing single excitations
        ! The normalization 1/n_{jb} used by Neufeld and Thom seems to be just a constant, we
        ! omit it for now
        implicit none

        call setup_elecs_sampler()
        call setup_holes_sampler()

    contains

        !------------------------------------------------------------------------------------------!

        subroutine setup_elecs_sampler()
            ! the probability distribution for selection of electrons
            ! creates a sampler for electron indices within the reference determinant
            ! these later have to be transferred to the current determinant

            ! even though strictly speaking, these are sums of the hole probabilities,
            ! expressing them in terms of the latter would be an unwanted dependency,
            ! since this relation is merely a matter of choice of the algorithm and should not
            ! be reflected in the design of the implementation
            implicit none
            real(dp) :: w(nel)
            integer :: i, a
            integer :: iEl

            do iEl = 1, nel
                w(iEl) = 0
                i = refDet(iEl)
                do a = 1, nBasis
                    w(iEl) = w(iEl) + acc_doub_matel(i, a)
                end do
            end do
            ! load the probabilites for electron selection into the alias table
            call single_elec_sampler%setup(root, w)

        end subroutine setup_elecs_sampler

        !------------------------------------------------------------------------------------------!

        subroutine setup_holes_sampler()
            ! the probability distributions for selection of holes, given the electron orbital
            implicit none
            real(dp) :: w(nBasis)
            integer :: i, a
            integer :: aerr

            ! there is one table for given source orbital
            allocate(single_hole_sampler(nBasis), stat=aerr)

            ! each table has probabilities for all given virtual orbitals a (some of them might
            ! be 0, this has no additional cost)
            do i = 1, nBasis
                w = 0.0_dp
                do a = 1, nBasis
                    ! we never want to sample the source orbital
                    if (i /= a) &
                        ! store the accumulated matrix elements (= un-normalized probability) with
                        ! the corresponding symmetry (if spins of a/i are different, w is 0)
                        w(a) = acc_doub_matel(i, a)
                end do

                call single_hole_sampler(i)%setup(root, w)
            end do

        end subroutine setup_holes_sampler

        !------------------------------------------------------------------------------------------!

        function acc_doub_matel(src, tgt) result(prob)
            ! Accumulate all single excitation matrix elements connecting
            ! D_(j,src)^(b,tgt) for all (j,b), where D is the reference
            ! Input: src - orbital we want to excite from (any orbital is allowed)
            !        tgt - orbital we want to excite to (any orbital is allowed)
            ! Output: prob - Accumulated matrix elements of singles
            !                If src/tgt have different spin, the result is 0, if they have
            !                different symmetry, it is not necessarily
            !
            ! IMPORTANT: The matrix elements are calculated as
            ! <D_j^b|H|D_(j,src)^(b,tgt)> := h_src^tgt + sum_(k occ in D_j^b) <tgt k|src k> - <tgt k|k src>
            ! That is, the left side is to be understood symbolic, there is no actual excitation
            ! from src to tgt
            ! symmetry has to be preserved
            implicit none
            integer, intent(in) :: src, tgt
            real(dp) :: prob
            integer :: b, j
            integer :: nI(nel)
            type(Excite_1_t) :: ex
            logical :: tPar

            prob = 0

            if (symAllowed(src, tgt)) then
                do b = 1, nBasis
                    ! loop over all non-occupied orbitals
                    if (.not. any(b == refDet(:))) then
                        do j = 1, nel
                            ! get the excited determinant D_j^b used for the matrix element
                            if (symAllowed(refDet(j), b)) then
                                call make_single(refDet(:), nI, j, b, ex%val, tPar)
                                ! this is a symbolic excitation, we do NOT require src to be occupied
                                ! we just use the formula for excitation matrix elements
                                prob = prob + abs(sltcnd_excit(nI, Excite_1_t(src, tgt), tPar, assert_occupation=.false.))
                            end if
                        end do
                    end if
                end do
            end if

        end function acc_doub_matel
        !------------------------------------------------------------------------------------------!
    end subroutine init_pcpp_singles_excitgen