init_pcpp_doubles_excitgen Subroutine

    subroutine init_pcpp_doubles_excitgen()
        implicit none

        call setup_elec_one_sampler()
        call setup_elec_two_sampler()

        call setup_hole_one_sampler()
        call setup_hole_two_sampler()
    contains

        subroutine setup_elec_one_sampler()
            integer :: i
            integer :: a, b, j
            real(dp) :: w(nel)
            logical :: tPar
            integer :: iEl, jEl

            tPar = .false.
            do iEl = 1, nel
                w(iEl) = 0
                do jEl = 1, nel
                    if (iEl /= jEl) then
                        i = refDet(iEl)
                        j = refDet(jEl)
                        do a = 1, nBasis
                            if (.not. any(a == [i, j])) then
                                do b = 1, nBasis
                                    if (.not. any(b == [a, i, j]) &
                                            .and. calc_spin_raw(i) + calc_spin_raw(j) == calc_spin_raw(a) + calc_spin_raw(b)) then
                                        w(iEl) = w(iEl) + abs(sltcnd_excit(refDet, Excite_2_t(i, a, j, b), tPar, assert_occupation=.false.))
                                    end if
                                end do
                            end if
                        end do
                    end if
                end do
            end do

            call apply_lower_bound(w)
            call double_elec_one_sampler%setup(root, w)
        end subroutine setup_elec_one_sampler

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

        subroutine setup_elec_two_sampler()
            implicit none
            real(dp) :: w(nel)
            type(Excite_2_t) :: ex
            logical :: tPar
            integer :: aerr
            integer :: i, j, a, b
            integer :: jEl

            allocate(double_elec_two_sampler(nBasis), stat=aerr)
            tPar = .false.
            do i = 1, nBasis
                w = 0.0_dp
                do jEl = 1, nel
                    j = refDet(jEl)
                    if (i /= j) then
                        do a = 1, nBasis
                            if (.not. any(a == [i, j])) then
                                do b = 1, nBasis
                                    if (.not. any(b == [a, i, j]) &
                                            .and. calc_spin_raw(i) + calc_spin_raw(j) == calc_spin_raw(a) + calc_spin_raw(b)) then
                                            w(jEl) = w(jEl) + abs(sltcnd_excit(refDet, Excite_2_t(i, a, j, b), tPar, assert_occupation=.false.))
                                    end if
                                end do
                            end if
                        end do
                    end if
                end do

                ! to prevent bias, a lower bound for the probabilities is set
                call apply_lower_bound(w)

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

        end subroutine setup_elec_two_sampler
        !------------------------------------------------------------------------------------------!

        subroutine setup_hole_one_sampler()
            ! generate precomputed probabilities for picking a hole given a selected electron
            ! this is for picking the first hole where no symmetry restrictions apply
            implicit none
            real(dp) :: w(nBasis, 0:spinMax)
            integer :: i, a, iSpin
            integer :: aerr

            allocate(double_hole_one_sampler(nBasis, 0:spinMax), stat=aerr)
            do i = 1, nBasis
                w = 0.0_dp
                do a = 1, nBasis
                    ! we will be requesting orbitals with a defined spin, store it along
                    if (a /= i) &
                        w(a, getSpinIndex(a)) = pp_weight_function(i, a)
                end do

                do iSpin = 0, spinMax
                    call double_hole_one_sampler(i, iSpin)%setup(root, w(:, iSpin))
                end do
            end do
        end subroutine setup_hole_one_sampler

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

        subroutine setup_hole_two_sampler()
            ! generate precomputed probabilities for picking hole number 2 given a selected electron
            ! this is for picking the second hole where symmetry restrictions apply
            implicit none
            real(dp) :: w(nBasis, 0:symmax - 1, 0:spinMax)
            integer :: j, b, iSym, iSpin
            integer :: aerr

            ! there is one table for each symmetry and each starting orbital
            allocate(double_hole_two_sampler(nBasis, 0:symmax - 1, 0:spinMax), stat=aerr)
            do j = 1, nBasis
                w = 0.0_dp
                do b = 1, nBasis
                    ! only same-spin and symmetry-allowed excitations from j -> b
                    if (b /= j) &
                        w(b, G1(b)%Sym%s, getSpinIndex(b)) = pp_weight_function(j, b)
                end do

                do iSpin = 0, spinMax
                    do iSym = 0, symmax - 1
                        call double_hole_two_sampler(j, iSym, iSpin)%setup(root, w(:, iSym, iSpin))
                    end do
                end do
            end do

        end subroutine setup_hole_two_sampler

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

    end subroutine init_pcpp_doubles_excitgen