calcFullstopLoweringStochastic Subroutine

    subroutine calcFullstopLoweringStochastic(ilut, csf_i, excitInfo, t, branch_pgen, &
                                              posSwitches, negSwitches, opt_weight)
        integer(n_int), intent(in) :: ilut(0:nifguga)
        type(CSF_Info_t), intent(in) :: csf_i
        type(ExcitationInformation_t), intent(in) :: excitInfo
        integer(n_int), intent(out) :: t(0:nifguga)
        real(dp), intent(out) :: branch_pgen
        real(dp), intent(in) :: posSwitches(nSpatOrbs), negSwitches(nSpatOrbs)
        type(WeightObj_t), intent(in), optional :: opt_weight
        character(*), parameter :: this_routine = "calcFullstopLoweringStochastic"

        real(dp) :: nOpen, tempWeight, bVal, temp_pgen
        HElement_t(dp) :: umat
        type(WeightObj_t) :: weights
        integer :: iOrb, deltaB

        ASSERT(isZero(ilut, excitInfo%fullEnd))
        ASSERT(isProperCSF_ilut(ilut))

        ! could check U matrix element here.. although in the final version
        ! of the orbital picker, it probably should be accessed already
        ! for the cauchy-schwarz criteria
        umat = (get_umat_el(excitInfo%fullEnd, excitInfo%fullEnd, &
                            excitInfo%fullStart, excitInfo%secondStart) + &
                get_umat_el(excitInfo%fullEnd, excitInfo%fullEnd, &
                            excitInfo%secondStart, excitInfo%fullStart)) / 2.0_dp

        ! todo correct combination of sum terms and correct indexcing
        if (near_zero(umat)) then
            branch_pgen = 0.0_dp
            t = 0_n_int
            return
        end if

        ! create weight object here
        ! i think i only need single excitations weights here, since
        ! the semi stop in this case is like an end step...
        if (present(opt_weight)) then
            weights = opt_weight
        else
            weights = init_singleWeight(csf_i, excitInfo%secondStart)
        end if

        ! i only need normal single stochastic start then..
        call createStochasticStart_single(ilut, csf_i, excitInfo, weights, posSwitches, &
                                          negSwitches, t, branch_pgen)

        ! matrix elements cant be 0, but maybe both branch weights were...
        check_abort_excit(branch_pgen, t)

        ! then stochastc single update
        do iOrb = excitInfo%fullStart + 1, excitInfo%secondStart - 1
            call singleStochasticUpdate(ilut, csf_i, iOrb, excitInfo, weights, posSwitches, &
                                        negSwitches, t, temp_pgen)
            branch_pgen = branch_pgen * temp_pgen
            ! matrix elements cant be 0 but branch weight might be..
            check_abort_excit(branch_pgen, t)
        end do

        ! write it for specific lowering semi-start
        iOrb = excitInfo%secondStart
        bVal = csf_i%B_real(iOrb)

        ! hope that the weighs correctly assert that only fitting deltaB
        ! values arrive here, to ensure a deltaB=0 branch in the DE overlap
        deltaB = getDeltaB(t)

        select case (csf_i%stepvector(iOrb))
        case (1)
            ASSERT(getDeltaB(t) == -1)

            ! change 1 -> 0
            clr_orb(t, 2 * iOrb - 1)

            call getDoubleMatrixElement(0, 1, deltaB, gen_type%L, gen_type%L, &
                                        bVal, excitInfo%order, tempWeight)

            call setDeltaB(0, t)

        case (2)
            ASSERT(getDeltaB(t) == 1)

            ! change 2 -> 0
            clr_orb(t, 2 * iOrb)

            call getDoubleMatrixElement(0, 2, deltaB, gen_type%L, gen_type%L, &
                                        bVal, 1.0_dp, tempWeight)

            call setDeltaB(0, t)

        case (3)
            ASSERT(isThree(ilut, iOrb))

            if (deltaB == -1) then

                ! change 3->2
                clr_orb(t, 2 * iOrb - 1)

                call getDoubleMatrixElement(2, 3, deltaB, gen_type%L, gen_type%L, &
                                            bVal, 1.0_dp, tempWeight)

            else
                ! change 3->1
                clr_orb(t, 2 * iOrb)

                call getDoubleMatrixElement(1, 3, deltaB, gen_type%L, gen_type%L, &
                                            bVal, 1.0_dp, tempWeight)

            end if

            call setDeltaB(0, t)

#ifdef DEBUG_
        case (0)
            call stop_all(this_routine, "wrong stepvalue!")
#endif
        end select

        ! only deltaB = 1 branch, so only number of open orbitals important
        ! for matrix element in overlap region

        nOpen = (-1.0_dp)**real(count_open_orbs_ij(csf_i, excitInfo%secondStart + 1, &
                                                   excitInfo%fullEnd - 1), dp)

        iOrb = excitInfo%fullEnd

        ! change 0 -> 3
        set_orb(t, 2 * iOrb)
        set_orb(t, 2 * iOrb - 1)

        call update_matrix_element(t, tempWeight * nOpen * Root2, 1)

        if (tFillingStochRDMOnFly) then
            call encode_stochastic_rdm_info(GugaBits, t, rdm_ind= &
                                            contract_2_rdm_ind(excitInfo%i, excitInfo%j, excitInfo%k, excitInfo%l, &
                                                               excit_lvl=2, excit_typ=excitInfo%typ), x1=0.0_dp, &
                                            x0=extract_matrix_element(t, 1))
        end if

        ! update all matrix element contributions at once
        call encode_matrix_element(t, 0.0_dp, 2)
        call update_matrix_element(t, umat, 1)

    end subroutine calcFullstopLoweringStochastic