calc_all_excits_guga_rdm_singles Subroutine

    subroutine calc_all_excits_guga_rdm_singles(ilut, csf_i, i, j, excits, n_excits)
        integer(n_int), intent(in) :: ilut(0:GugaBits%len_tot)
        type(CSF_Info_t), intent(in) :: csf_i
        integer, intent(in) :: i, j
        integer(n_int), intent(out), allocatable :: excits(:, :)
        integer, intent(out) :: n_excits
        character(*), parameter :: this_routine = "calc_all_excits_guga_rdm_singles"

        type(ExcitationInformation_t) :: excitInfo
        type(WeightObj_t) :: weights
        real(dp) :: posSwitches(nSpatOrbs), negSwitches(nSpatOrbs)
        integer :: iEx, iOrb
        integer(n_int), allocatable :: tempExcits(:, :)
        real(dp) :: minusWeight, plusWeight

        ASSERT(i > 0 .and. i <= nSpatOrbs)
        ASSERT(j > 0 .and. j <= nSpatOrbs)

        n_excits = 0

        excitInfo = excitationIdentifier(i, j)

        associate (gen1 => excitInfo%gen1, en => excitInfo%fullEnd, &
                   st => excitInfo%fullStart)

            if (gen1 /= 0) then
                if (csf_i%stepvector(i) == 3 .or. csf_i%stepvector(j) == 0) then
                    allocate(excits(0, 0))
                    return
                end if
            end if

            call calcRemainingSwitches_excitInfo_single(csf_i, excitInfo, posSwitches, negSwitches)

            weights = init_singleWeight(csf_i, en)
            plusWeight = weights%proc%plus(posSwitches(st), csf_i%B_real(st), weights%dat)
            minusWeight = weights%proc%minus(negSwitches(st), csf_i%B_real(st), weights%dat)

            ! check compatibility of chosen indices

            if (csf_i%stepvector(st) == 1 .and. near_zero(plusWeight) &
                    .or. csf_i%stepvector(st) == 2 .and. near_zero(minusWeight) &
                    .or. near_zero(minusWeight + plusWeight)) then
                allocate(excits(0, 0))
                return
            end if

            ! have to give probabilistic weight object as input, to deal
            call createSingleStart(ilut, csf_i, excitInfo, posSwitches, &
                                   negSwitches, weights, tempExcits, n_excits)

            do iOrb = st + 1, en - 1
                call singleUpdate(ilut, csf_i, iOrb, excitInfo, posSwitches, &
                                  negSwitches, weights, tempExcits, n_excits)
            end do

            call singleEnd(ilut, csf_i, excitInfo, tempExcits, &
                           n_excits, excits)

            ! encode the combined RDM-ind in the deltaB position for
            ! communication purposes
            do iEx = 1, n_excits
                ! just for consistency all encode the x0 element in the
                ! ind_rdm_x0 entry
                call encode_stochastic_rdm_info(GugaBits, excits(:, iEx), &
                                                rdm_ind=contract_1_rdm_ind(i, j), &
                                                x0=extract_matrix_element(excits(:, iex), 1), x1=0.0_dp)
            end do

        end associate

    end subroutine calc_all_excits_guga_rdm_singles