guga_excitations.F90 Source File


Contents

Source Code


Source Code

#include "macros.h"
! GUGA excitations module:
! contains as much excitation related functionality as possible
module guga_excitations
    ! modules
    use CalcData, only: t_trunc_guga_pgen, t_trunc_guga_matel, &
                        trunc_guga_pgen, trunc_guga_matel, t_trunc_guga_pgen_noninits, &
                        t_matele_cutoff, matele_cutoff

    use SystemData, only: nEl, nBasis, ElecPairs, G1, nSpatOrbs, &
                          tGen_guga_weighted, nBasisMax, treal, &
                          t_approx_exchange, t_approx_exchange_noninits, &
                          is_init_guga, t_heisenberg_model, t_tJ_model, t_mixed_hubbard, &
                          t_guga_pchb, thub

    use constants, only: dp, n_int, bits_n_int, Root2, &
                         EPS, bni_, bn2_, stdout, int_rdm

    use bit_reps, only: decode_bit_det

    use bit_rep_data, only: GugaBits, niftot, nifguga, nifd

    use DetBitOps, only: count_open_orbs, return_ms

    use guga_data, only: ExcitationInformation_t, getSingleMatrixElement, &
                         getDoubleMatrixElement, getMixedFullStop, &
                         WeightData_t, orbitalIndex, &
                         tag_excitations, tag_tmp_excits, &
                         excit_type, gen_type, excit_names

    use guga_bitRepOps, only: isProperCSF_ilut, calcB_vector_ilut, getDeltaB, &
                              setDeltaB, count_open_orbs_ij, &
                              encode_matrix_element, update_matrix_element, &
                              extract_matrix_element, write_det_guga, &
                              write_guga_list, add_guga_lists, count_alpha_orbs_ij, &
                              count_beta_orbs_ij, findFirstSwitch, findLastSwitch, &
                              calcStepvector, find_switches, &
                              calcOcc_vector_int, &
                              extract_h_element, encode_stochastic_rdm_info, &
                              get_preceeding_opposites, &
                              CSF_Info_t, csf_ref

    use OneEInts, only: GetTMatEl

    use procedure_pointers, only: get_umat_el

    use dSFMT_interface, only: genrand_real2_dSFMT

    use FciMCData, only: ilutRef, tFillingStochRDMOnFly

    use util_mod, only: binary_search_first_ge, abs_l1, operator(.isclose.), &
                        operator(.div.), near_zero, stop_all, neci_flush

    use GenRandSymExcitNUMod, only: RandExcitSymLabelProd

    use SymExcitDataMod, only: SpinOrbSymLabel, OrbClassCount, SymLabelCounts2, &
                               KPointToBasisFn, sym_label_list_spat

    use sym_general_mod, only: ClassCountInd

    use excit_gens_int_weighted, only: get_paired_cc_ind

    use umatcache, only: gtID

    use guga_procedure_pointers, only: &
                calc_orbital_pgen_contr, &
                calc_orbital_pgen_contrib_start, calc_orbital_pgen_contrib_end

    use Umatcache, only: UMat2D

    use timing_neci, only: timer, set_timer, halt_timer, get_total_time

    use guga_types, only: WeightObj_t

    use MemoryManager, only: LogMemAlloc, LogMemDealloc

    use sym_mod, only: MomPbcSym

    use guga_bitRepOps, only: contract_2_rdm_ind

    use lattice_models_utils, only: create_all_open_shell_dets

    use guga_matrixElements, only:  calc_mixed_start_contr_sym, calc_mixed_end_contr_sym, &
        calc_mixed_contr_sym, get_forced_zero_double, getminus_double, getminus_semistart, &
        getplus_double, getplus_semistart, init_doubleweight, &
        init_semistartweight, calc_guga_matrix_element, calc_mixed_contr_integral, &
        calcremainingswitches_excitinfo_double, calcremainingswitches_excitinfo_single, &
        calcstartprob, calcstayingprob, endfx, endgx, &
        init_fullstartweight, init_singleweight, calc_mixed_start_contr_integral, &
        calc_mixed_end_contr_integral, calc_mixed_contr_pgen, calc_mixed_start_contr_pgen, &
        calc_mixed_end_contr_pgen

    better_implicit_none

    private
    public :: global_excitinfo, print_excitinfo, &
              assign_excitinfo_values_double, assign_excitinfo_values_single, &
              actHamiltonian, calcdoubleexcitation_withweight, &
              calcnonoverlapdouble, calcsingleoverlaplowering, calcsingleoverlapraising, &
              calcsingleoverlapmixed, calcdoublelowering, calcdoubleraising, &
              calcdoublel2r, calcdoubler2l, calcfullstoplowering, calcfullstopraising, &
              calcfullstopl2r, calcfullstopr2l, calcfullstartlowering, &
              calcfullstartraising, calcfullstartl2r, calcfullstartr2l, &
              calcfullstartfullstopalike, calcfullstartfullstopmixed, &
              calcremainingswitches_excitinfo_double, checkcompatibility, &
              createsinglestart, singleupdate, singleend, init_singleweight, &
              calcremainingswitches_excitinfo_single, excitationIdentifier, &
              pickorbs_sym_uniform_ueg_single, pickorbs_sym_uniform_ueg_double, &
              pickorbs_sym_uniform_mol_single, pickorbs_sym_uniform_mol_double, &
              calc_orbital_pgen_contr_ueg, calc_orbital_pgen_contr_mol, &
              calc_mixed_x2x_ueg, &
              calc_off_diag_guga_ref_direct, &
              pickorbs_real_hubbard_single, pickorbs_real_hubbard_double, &
              excitationIdentifier_single, excitationIdentifier_double, &
              init_doubleWeight, init_semiStartWeight, init_fullStartWeight, &
              calcremainingswitches_single, calcallexcitations_single, &
              calcallexcitations_double, &
              createstochasticstart_single, &
              pickrandomorb_scalar, pickrandomorb_forced, pickrandomorb_vector, &
              pickrandomorb_restricted, singlestochasticupdate, &
              singlestochasticend, &
              calcfullstartfullstopmixedstochastic, mixedfullstartstochastic, &
              doubleupdatestochastic, calcfullstartraisingstochastic, &
              calcfullstartloweringstochastic, calcfullstopraisingstochastic, &
              calcfullstoploweringstochastic, calcsingleoverlapmixedstochastic, &
              mixedfullstopstochastic, calcloweringsemistartstochastic, &
              calcraisingsemistartstochastic, calcloweringsemistopstochastic, &
              calcraisingsemistopstochastic, calcfullstartl2r_stochastic, &
              calcfullstartr2l_stochastic, calcfullstopl2r_stochastic, &
              calcfullstopr2l_stochastic, calcdoubleloweringstochastic, &
              calcdoubleraisingstochastic, calcdoublel2r2l_stochastic, &
              calcdoubler2l2r_stochastic, calcdoublel2r_stochastic, &
              calcdoubler2l_stochastic, calcallexcitations, &
              pick_elec_pair_uniform_guga, get_guga_integral_contrib, &
              calc_pgen_mol_guga_single, get_excit_level_from_excitInfo, &
              get_guga_integral_contrib_spat, calc_orbital_pgen_contrib_start_def, &
              calc_orbital_pgen_contrib_end_def, create_hamiltonian_guga, &
              csf_to_sds_ilut, csf_vector_to_sds, checkCompatibility_single

    ! use a global excitationInformation type variable to store information
    ! about the last generated excitation to analyze matrix elements and
    ! pgens
    type(ExcitationInformation_t) :: global_excitInfo

    abstract interface
        function calc_pgen_general(csf_i, i) result(pgen)
            import :: dp, CSF_Info_t
            type(CSF_Info_t), intent(in) :: csf_i
            integer, intent(in) :: i
            real(dp) :: pgen
        end function calc_pgen_general
    end interface

    interface excitationIdentifier
        module procedure excitationIdentifier_single
        module procedure excitationIdentifier_double
    end interface excitationIdentifier

    interface calcAllExcitations
        module procedure calcAllExcitations_single
        module procedure calcAllExcitations_double
        module procedure calcAllExcitations_excitInfo_single
    end interface calcAllExcitations

contains

    subroutine csf_vector_to_sds(csfs, csf_coeffs, sds, sd_coeffs, ms)
        integer(n_int), intent(in) :: csfs(:,:)
        real(dp), intent(in) :: csf_coeffs(:)
        real(dp), intent(in), optional :: ms
        integer(n_int), intent(out), allocatable :: sds(:,:)
        real(dp), intent(out), allocatable :: sd_coeffs(:)

        real(dp) :: ms_
        integer :: n_sds, spin, n_tot, i
        integer(n_int), allocatable :: temp_all(:,:), temp_sds(:,:)
        real(dp), allocatable :: temp_coeffs(:)


        spin = abs(return_ms(csfs(:,1)))
        def_default(ms_, ms, spin/2.)

        n_sds = 2 ** nSpatorbs
        allocate(temp_all(0:GugaBits%len_tot,n_sds), source = 0_n_int)

        n_tot = 0

        do i = 1, size(csfs,2)
            call csf_to_sds_ilut(csfs(:,i), temp_sds, temp_coeffs, ms_, csf_coeffs(i))
            call add_guga_lists(n_tot, size(temp_sds,2), temp_all, temp_sds)
        end do

        allocate(sds(0:GugaBits%len_tot, n_tot), source = temp_all(:,1:n_tot))
        allocate(sd_coeffs(n_tot), source = 0.0_dp)

        do i = 1, n_tot
            sd_coeffs(i) = extract_matrix_element(sds(:,i), 1)
        end do

    end subroutine csf_vector_to_sds

    subroutine csf_to_sds_ilut(csf, sds, weights, ms, coeff)
        integer(n_int), intent(in) :: csf(0:GugaBits%len_tot)
        real(dp), intent(in), optional :: ms
        real(dp), intent(in), optional :: coeff
        integer(n_int), intent(out), allocatable :: sds(:,:)
        real(dp), intent(out), allocatable :: weights(:)
        character(*), parameter :: this_routine = "csf_to_sds_ilut"

        real(dp) :: ms_
        integer :: spin, n_alpha, n_beta
        integer :: i, j, step(nSpatorbs), delta_k
        integer :: nI(nel)
        real(dp) :: bVec(nSpatorbs), aVec(nSpatorbs), lambda_k, x, coeff_
        integer(n_int), allocatable :: all_sds(:,:)
        real(dp), allocatable :: all_weights(:)

        def_default(coeff_, coeff, 1.0_dp)


        ! only works for heisenberg model for now..
        if (any(calcOcc_vector_int(csf) /= 1)) then
            call stop_all(this_routine, "only implemented for heisenberg for now")
        end if
        if (near_zero(coeff_)) then
            allocate(sds(0:GugaBits%len_tot, 0))
            allocate(weights(0))
            return
        end if
        spin = abs(return_ms(csf))

        def_default(ms_, ms, spin/2.)

        ! construct SDs by attaching (N/2 + ms) up spins and (N/2 - ms) down spins
        n_alpha = nint(nSpatorbs / 2. + ms_)
        n_beta = nint(nSpatorbs / 2. - ms_)

        all_sds = create_all_open_shell_dets(nSpatorbs, n_beta, n_alpha)

        allocate(all_weights(size(all_sds,2)), source = 0.0_dp)

        step = calcStepvector(csf(0:GugaBits%len_orb))
        bVec = calcB_vector_ilut(csf(0:GugaBits%len_orb))
        aVec = real(([(i, i = 1,nSpatorbs)] - bVec),dp) / 2.0

        do i = 1, size(all_sds,2)
            x = 1.0_dp

            call decode_bit_det(nI, all_sds(:,i))

            do j = 1, nSpatorbs
                if (.not. near_zero(x)) then
                    delta_k = mod(get_spin(nI(j)),2)
                    lambda_k = get_preceeding_opposites(nI, j)

                    select case (step(j))

                    case (1)
                        x = x * sqrt((aVec(j) + bVec(j) - lambda_k)/bVec(j))

                    case (2)
                        x = x * (-1.0_dp)**(bVec(j)+delta_k) * &
                            sqrt((lambda_k - aVec(j) + 1.0_dp)/(bVec(j)+2.0_dp))

                    case (3)

                        x = x * (-1.0_dp)**bVec(j)

                    end select

                end if
            end do

            all_weights(i) =  x

        end do


        j = 1
        do i = 1, size(all_sds,2)
            if (.not. near_zero(all_weights(i))) then
                all_sds(:,j) = all_sds(:,i)
                all_weights(j) = all_weights(i)
                j = j + 1
            end if
        end do

        ! i might have to normalize
        all_weights(1:j-1) = coeff_ * all_weights(1:j-1) / sqrt(sum(all_weights(1:j-1)**2))

        allocate(sds(0:GugaBits%len_tot, j - 1), source = 0_n_int)
        allocate(weights(j-1), source = all_weights(1:j-1))

        do i = 1, j - 1
            sds(0:0, i) = all_sds(:, i)
            call encode_matrix_element(sds(:,i), all_weights(i), 1)
        end do

    end subroutine csf_to_sds_ilut


    function calc_off_diag_guga_ref_direct(ilut, csf_i, run, exlevel) result(hel)
        integer(n_int), intent(in) :: ilut(0:niftot)
        type(CSF_Info_t), intent(in) :: csf_i
        integer, intent(in), optional :: run
        integer, intent(out), optional :: exlevel
        HElement_t(dp) :: hel

        integer :: run_
        type(ExcitationInformation_t) :: excitInfo

        def_default(run_, run, 1)
        call calc_guga_matrix_element(ilut, csf_i, ilutRef(0:niftot, run_), csf_ref(run_), excitInfo, hel, .true.)

        if (present(exlevel)) then
            if (excitInfo%valid) then
                exlevel = merge(1, 2, excitInfo%typ == excit_type%single)
            else
                ! non-valid > 3 excit
                exlevel = nel
            end if
        end if
    end function calc_off_diag_guga_ref_direct

    function calc_guga_mat_wrapper(ilutI, csf_i, ilutJ, csf_j) result(mat_ele)
        integer(n_int), intent(in) :: ilutI(0:niftot), ilutJ(0:niftot)
        type(CSF_Info_t), intent(in) :: csf_i, csf_j
        HElement_t(dp) :: mat_ele

        type(ExcitationInformation_t) :: excitInfo

        call calc_guga_matrix_element(ilutI, csf_i, ilutJ, csf_j, excitInfo, mat_ele, t_hamil = .true.)

    end function calc_guga_mat_wrapper

    function create_hamiltonian_guga(ilut_list) result(hamil)
        integer(n_int), intent(in) :: ilut_list(:,:)
        HElement_t(dp) :: hamil(size(ilut_list,2), size(ilut_list,2))

        type(CSF_Info_t) :: csf_i, csf_j
        integer :: i, j

        do i = 1, size(ilut_list,2)
            csf_i = CSF_Info_t(ilut_list(:, i))
            do j = 1, size(ilut_list,2)
                csf_j = CSF_Info_t(ilut_list(:, j))
                hamil(i,j) = calc_guga_mat_wrapper(ilut_list(:, j), csf_j, ilut_list(:,i), csf_i)
            end do
        end do

    end function create_hamiltonian_guga


    ! write up all the specific stochastic excitation routines

    subroutine calcFullStartFullStopMixedStochastic(ilut, csf_i, excitInfo, t, 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(inout) :: excitInfo
        integer(n_int), intent(out) :: t(0:nifguga)
        real(dp), intent(out) :: pgen
        real(dp), intent(in) :: posSwitches(nSpatOrbs), negSwitches(nSpatOrbs)
        type(WeightObj_t), intent(in), optional :: opt_weight
        character(*), parameter :: this_routine = "calcFullStartFullStopMixedStochastic"

        type(WeightObj_t) :: weights
        real(dp) ::  branch_pgen, temp_pgen, above_cpt, below_cpt, rdm_mat, p_orig
        HElement_t(dp) :: integral
        integer :: iOrb, i, j, k, l, typ

        ASSERT(.not. isZero(ilut, excitInfo%fullStart))
        ASSERT(.not. isThree(ilut, excitInfo%fullStart))
        ASSERT(.not. isZero(ilut, excitInfo%fullEnd))
        ASSERT(.not. isThree(ilut, excitInfo%fullEnd))

        if (present(opt_weight)) then
            weights = opt_weight
        else
            if (t_approx_exchange .or. (t_approx_exchange_noninits .and. &
                    (.not. is_init_guga))) then
                weights = init_forced_end_exchange_weight(csf_i, excitInfo%fullEnd)
            else
                weights = init_doubleWeight(csf_i, excitInfo%fullEnd)
            end if
        end if

        if (t_approx_exchange .or. (t_approx_exchange_noninits .and. &
                (.not. is_init_guga))) then
            call forced_mixed_start(ilut, csf_i, excitInfo, t, branch_pgen)
        else
            call mixedFullStartStochastic(ilut, csf_i, excitInfo, weights, posSwitches, &
                                          negSwitches, t, branch_pgen)
        end if

        ! set pgen to 0 in case of early exit
        pgen = 0.0_dp

        ! in this excitation there has to be a switch somewhere in the double
        ! overlap region, so only x1 matrix element counts essentially
        ! and this can be 0 at the full-start already -> so check that here
        ! and in case abort
        ! should i do that in the mixedFullStartStochastic routine

        ! do that x1 matrix element in the routine and only check probWeight here
        check_abort_excit(branch_pgen, t)

        temp_pgen = 1.0_dp

        do iOrb = excitInfo%fullStart + 1, excitInfo%fullEnd - 1
            call doubleUpdateStochastic(ilut, csf_i, iOrb, excitInfo, weights, negSwitches, &
                                        posSwitches, t, branch_pgen)

            ! zero x1 - elements can also happen in the double update
            if (near_zero(extract_matrix_element(t, 2)) .or. near_zero(branch_pgen)) then
                t = 0_n_int
                return
            end if
        end do

        call mixedFullStopStochastic(ilut, csf_i, excitInfo, t)

        ! check if there was a change in the stepvector in the double
        ! overlap region
        if (.not. near_zero(extract_matrix_element(t, 1))) then
            t = 0_n_int
            return
        end if

        ! if we do RDMs also store the x0 and x1 coupling coeffs
        ! and I need to do it before the routines below since excitInfo
        ! gets changed there
        if (tFillingStochRDMOnFly) then
            ! i need to unbias against the total pgen later on in the
            ! RDM sampling otherwise the rdm-bias factor is not correct!
            ! encode the necessary information in the rdm-matele!
            i = excitInfo%i
            j = excitInfo%j
            k = excitInfo%k
            l = excitInfo%l
            typ = excitInfo%typ
            rdm_mat = extract_matrix_element(t, 2)
            call calc_orbital_pgen_contr(csf_i, [2 * i, 2 * j], above_cpt, &
                                         below_cpt)
            p_orig = (above_cpt + below_cpt) * branch_pgen
            if (.not. (t_heisenberg_model .or. t_tJ_model)) then
                p_orig = p_orig / real(ElecPairs, dp)
            end if
        end if

        global_excitInfo = excitInfo

        if (t_approx_exchange .or. (t_approx_exchange_noninits .and. (.not. is_init_guga))) then

            if (getDeltaB(t) == 0) then
                t = 0_n_int
                return
            end if

            call calc_mixed_contr_integral(ilut, csf_i, t, excitInfo%fullStart, &
                excitInfo%fullEnd, integral)

            pgen = branch_pgen

            ! just to be save that a switch always happens at the end
            ! print that out for now
        else
            call calc_mixed_contr_integral(ilut, csf_i, t, excitInfo%fullStart, &
                excitInfo%fullEnd, integral)
            if (.not. near_zero(integral)) then
                call calc_mixed_contr_pgen(ilut, csf_i, t, excitInfo, pgen)
            end if

!             call calc_mixed_contr_sym(ilut, csf_i, t, excitInfo, pgen, integral)
        end if

        if (near_zero(integral)) then
            t = 0_n_int
            pgen = 0.0_dp
            return
        end if

        if (tFillingStochRDMOnFly) then
            if (.not. near_zero(p_orig)) then
                call encode_stochastic_rdm_info(GugaBits, t, rdm_ind= &
                           contract_2_rdm_ind(i, j, k, l, excit_lvl = 2, &
                           excit_typ=typ), x0 = 0.0_dp, x1 = rdm_mat * pgen / p_orig)
            end if
        end if

        call encode_matrix_element(t, 0.0_dp, 2)
        call encode_matrix_element(t, integral, 1)

    end subroutine calcFullStartFullStopMixedStochastic

    subroutine calc_orbital_pgen_contr_mol(csf_i, occ_orbs, cpt_a, cpt_b)
        ! calculates the cumulatice probability list for different
        ! full-start -> full-stop mixed excitations, used in the recalculation of
        ! contrbuting pgens from different picked orbitals
        type(CSF_Info_t), intent(in) :: csf_i
        integer, intent(in) :: occ_orbs(2)
        real(dp), intent(out) :: cpt_a, cpt_b

        integer :: i, j, orb
        real(dp) :: cum_sum, cpt_ab, cpt_ba, ba_sum, ab_sum

        ! given 2 already picked electrons, this routine creates a list of
        ! p(a)*p(b|a) probabilities to pick the already determined holes
        ! is this so easy??
        ! what do i know? i know there is a possible switch between i and j
        ! or otherwise those indices would not have been taken, since an
        ! excitation was already created
        ! i know electrons i and j are from different orbitals and singly
        ! occupied
        ! have to loop once over all orbitals to create list for p(x)
        ! most def, but that should be easy since there are not so many
        ! restrictions, but to determine p(b|a) i have to take into account
        ! the symmetries and the guga restrictions..
        ! but i need the whole information to get cum_sum correct, for the
        ! given electrons i and j..
        ! maybe split up the loop to different sectors, where i know more info

        ! UPDATE: calculate it more directly! also for UEG models!
        ! just output the nececcary p(a)*p(b|a) and vv. and not a whole
        ! list!

        ! chanve that routine now, to use the general pre-generated cum-list
        ! for the current CSF
        i = gtID(minval(occ_orbs))
        j = gtID(maxval(occ_orbs))

        cum_sum = 0.0_dp

        if (tGen_guga_weighted) then
            do orb = 1, i - 1
                ! calc. the p(a)
                if (csf_i%stepvector(orb) /= 3) then
                    cum_sum = cum_sum + get_guga_integral_contrib(occ_orbs, orb, -1)

                end if

            end do
        end if

        ! deal with orb (i) in specific way:
        cpt_a = get_guga_integral_contrib(occ_orbs, i, -1)

        cum_sum = cum_sum + cpt_a

        ! also get p(b|a)
        ! did i get that the wrong way around??
        call pgen_select_orb_guga_mol(csf_i, occ_orbs, i, j, cpt_ba, ba_sum, i, .true.)

        if (tGen_guga_weighted) then
            do orb = i + 1, j - 1
                if (csf_i%stepvector(orb) /= 3) then
                    cum_sum = cum_sum + get_guga_integral_contrib(occ_orbs, orb, -1)
                end if
            end do
        end if

        ! deal with j also speciallly
        cpt_b = get_guga_integral_contrib(occ_orbs, j, -1)

        cum_sum = cum_sum + cpt_b

        ! and get p(a|b)
        call pgen_select_orb_guga_mol(csf_i, occ_orbs, j, i, cpt_ab, ab_sum, -j, .true.)

        ! and deal with rest:

        if (tGen_guga_weighted) then
            do orb = j + 1, nSpatOrbs
                if (csf_i%stepvector(orb) /= 3) then
                    cum_sum = cum_sum + get_guga_integral_contrib(occ_orbs, orb, -1)
                end if
            end do
        end if

        if (.not. tGen_guga_weighted) then
            cum_sum = csf_i%cum_list(nSpatOrbs)
        end if

        if (near_zero(cum_sum) .or. near_zero(ab_sum) .or. near_zero(ba_sum)) then
            cpt_a = 0.0_dp
            cpt_b = 0.0_dp
        else
            ! and get hopefully correct final values:
            cpt_a = cpt_a / cum_sum * cpt_ba / ba_sum
            cpt_b = cpt_b / cum_sum * cpt_ab / ab_sum
        end if

    end subroutine calc_orbital_pgen_contr_mol

    subroutine calc_orbital_pgen_contr_ueg(csf_i, occ_orbs, above_cpt, below_cpt)
        type(CSF_Info_t), intent(in) :: csf_i
        integer, intent(in) :: occ_orbs(2)
        real(dp), intent(out) :: above_cpt, below_cpt

        real(dp) :: cum_sum, cum_arr(nSpatOrbs)
        type(ExcitationInformation_t) :: tmp_excitInfo(nSpatOrbs)
        integer :: tmp_orbArr(nSpatOrbs)
        integer :: i, j

        call gen_ab_cum_list_1_1(csf_i, occ_orbs, cum_arr, tmp_excitInfo, tmp_orbArr)
        cum_sum = cum_arr(nSpatOrbs)

        i = gtID(occ_orbs(1))
        j = gtID(occ_orbs(2))

        ! get orbital probability at beginning: question is can they be
        ! 0, and can they be independently 0, and if yes how is the
        ! influence...
        above_cpt = (cum_arr(j) - cum_arr(j - 1)) / cum_sum

        if (i == 1) then
            below_cpt = cum_arr(1) / cum_sum
        else
            below_cpt = (cum_arr(i) - cum_arr(i - 1)) / cum_sum
        end if

    end subroutine calc_orbital_pgen_contr_ueg

    subroutine calcDoubleR2L_stochastic(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(inout) :: 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 = "calcDoubleR2L_stochastic"

        integer :: iOrb, start2, ende1, ende2, start1, switch
        type(WeightObj_t) :: weights
        real(dp) :: temp_pgen
        HElement_t(dp) :: integral

        ASSERT(.not. isThree(ilut, excitInfo%fullStart))
        ASSERT(.not. isZero(ilut, excitInfo%secondStart))
        ASSERT(.not. isZero(ilut, excitInfo%firstEnd))
        ASSERT(.not. isThree(ilut, excitInfo%fullEnd))

        start1 = excitInfo%fullStart
        start2 = excitInfo%secondStart
        ende1 = excitInfo%firstEnd
        ende2 = excitInfo%fullEnd

        ! : create correct weights:
        if (present(opt_weight)) then
            weights = opt_weight
        else
            weights = init_fullDoubleWeight(csf_i, start2, ende1, ende2, negSwitches(start2), &
                                            negSwitches(ende1), posSwitches(start2), posSwitches(ende1), &
                                            csf_i%B_real(start2), csf_i%B_real(ende1))
        end if

        call createStochasticStart_single(ilut, csf_i, excitInfo, weights, posSwitches, &
                                          negSwitches, t, branch_pgen)

        ! check validity
        check_abort_excit(branch_pgen, t)

        do iOrb = start1 + 1, start2 - 1
            call singleStochasticUpdate(ilut, csf_i, iOrb, excitInfo, weights, posSwitches, &
                                        negSwitches, t, temp_pgen)

            branch_pgen = branch_pgen * temp_pgen
            ! check validity
            check_abort_excit(branch_pgen, t)

        end do

        ! change weights... maybe need both single and double type weights
        ! then do lowering semi start
        weights = weights%ptr

        call calcLoweringSemiStartStochastic(ilut, csf_i, excitInfo, weights, negSwitches, &
                                             posSwitches, t, branch_pgen)

        ! check validity
        check_abort_excit(branch_pgen, t)

        do iOrb = start2 + 1, ende1 - 1
            call doubleUpdateStochastic(ilut, csf_i, iOrb, excitInfo, weights, negSwitches, &
                                        posSwitches, t, branch_pgen)
            ! check validity
            check_abort_excit(branch_pgen, t)

        end do

        ! then update weights and and to lowering semi-stop
        weights = weights%ptr

        call calcRaisingSemiStopStochastic(ilut, csf_i, excitInfo, weights, negSwitches, &
                                           posSwitches, t, branch_pgen)

        ! check validity
        check_abort_excit(branch_pgen, t)

        excitInfo%currentGen = excitInfo%lastGen

        do iOrb = ende1 + 1, ende2 - 1
            call singleStochasticUpdate(ilut, csf_i, iOrb, excitInfo, weights, posSwitches, &
                                        negSwitches, t, temp_pgen)

            branch_pgen = branch_pgen * temp_pgen
            ! check validity
            check_abort_excit(branch_pgen, t)

        end do

        ! and finally to end step
        call singleStochasticEnd(csf_i, excitInfo, t)

        ! if we do RDMs also store the x0 and x1 coupling coeffs
        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), &
                                    x0=extract_matrix_element(t, 1), &
                                    x1=extract_matrix_element(t, 2))
        end if

        ! for the additional contributing integrals:
        ! i have to consider that there might be a non-overlap excitation,
        ! which leads to the same excitation if there is no change in the
        ! stepvector in the overlap region!
        ! where i have to parts of the matrix elements, but which can be
        ! expressed in terms of the already calulated one i think..
        ! atleast in the case if R2L and L2R, as both the bottom and top
        ! parts are the same and only the semi-start ans semi-stop have to
        ! be modified..
        ! and also the x1 element has to be dismissed i guess..
        ! so i should change how i deal with the x1 elements and maybe
        ! keep them out here in these functions..

        ! determine if a switch happended:
        switch = findFirstSwitch(ilut, t, start2 + 1, ende1)
        ! i have to rethink the matrix element contribution by the
        ! equivalent non-overlap excitations, since its possible that
        ! a +-2 excitations stays on the +-2 branch in the double overlap
        ! region all the time and then also can be produced by a non-overlap
        ! but then the matrix element is different then the calculation used
        ! right now..
        ! no!!! sinnce the excitations leading to such a csf wouldn be valid
        ! single excitations-.. since the deltaB values at the end would not
        ! match

        if (switch > 0) then
            ! a switch happened and only mixed overlap contributes
            ! after determining how to deal with different x0 and x1 parts
            ! wait: if a switch happened i know that the x0-element is 0
            ! and only the x1-element of the overlap region counts!
            integral = extract_matrix_element(t, 2) * (get_umat_el(start1, ende2, ende1, start2) + &
                                                       get_umat_el(ende2, start1, start2, ende1)) / 2.0_dp

            if (near_zero(integral)) then
                branch_pgen = 0.0_dp
                t = 0
            else
                call encode_matrix_element(t, 0.0_dp, 2)
                call encode_matrix_element(t, integral, 1)
            end if

        else
            ! no switch happened: so i have to think about the
            ! non-overlap contribution
            ! in the R2L and L2R case it is really easy if i keep the
            ! x0 and x1 elements seperate up until here.
            ! since the types of generators are the same and only the x0
            ! element(which does not get changed in the overlap region)
            ! are nececarry for the non-overlap excitation , it only gets
            ! a -t^2 = -1/2 factor...
            ! so to get the non-overlap version i need to multiply by -2.0!
            integral = (-extract_matrix_element(t, 1) * (get_umat_el(start1, ende2, start2, ende1) + &
                                                     get_umat_el(ende2, start1, ende1, start2)) * 2.0_dp + (extract_matrix_element(t, 1) + &
                                                              extract_matrix_element(t, 2)) * (get_umat_el(start1, ende2, ende1, start2) + &
                                                                                        get_umat_el(ende2, start1, start2, ende1))) / 2.0_dp

            if (near_zero(integral)) then
                branch_pgen = 0.0_dp
                t = 0_n_int
            else
                call encode_matrix_element(t, 0.0_dp, 2)
                call encode_matrix_element(t, integral, 1)
            end if
        end if
    end subroutine calcDoubleR2L_stochastic

    subroutine calcDoubleL2R_stochastic(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(inout) :: 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 = "calcDoubleL2R_stochastic"

        integer :: iOrb, start2, ende1, ende2, start1, switch
        type(WeightObj_t) :: weights
        real(dp) :: temp_pgen
        HElement_t(dp) :: integral

        ASSERT(.not. isZero(ilut, excitInfo%fullStart))
        ASSERT(.not. isThree(ilut, excitInfo%secondStart))
        ASSERT(.not. isThree(ilut, excitInfo%firstEnd))
        ASSERT(.not. isZero(ilut, excitInfo%fullEnd))

        start1 = excitInfo%fullStart
        start2 = excitInfo%secondStart
        ende1 = excitInfo%firstEnd
        ende2 = excitInfo%fullEnd

        ! : create correct weights:
        if (present(opt_weight)) then
            weights = opt_weight
        else
            weights = init_fullDoubleWeight(csf_i, start2, ende1, ende2, negSwitches(start2), &
                                            negSwitches(ende1), posSwitches(start2), posSwitches(ende1), &
                                            csf_i%B_real(start2), csf_i%B_real(ende1))
        end if

        call createStochasticStart_single(ilut, csf_i, excitInfo, weights, posSwitches, &
                                          negSwitches, t, branch_pgen)

        ! check validity
        check_abort_excit(branch_pgen, t)

        do iOrb = start1 + 1, start2 - 1
            call singleStochasticUpdate(ilut, csf_i, iOrb, excitInfo, weights, posSwitches, &
                                        negSwitches, t, temp_pgen)
            branch_pgen = branch_pgen * temp_pgen

            ! check validity
            check_abort_excit(branch_pgen, t)

        end do

        ! change weights... maybe need both single and double type weights
        ! then do lowering semi start
        weights = weights%ptr

        call calcRaisingSemiStartStochastic(ilut, csf_i, excitInfo, weights, negSwitches, &
                                            posSwitches, t, branch_pgen)

        ! check validity
        check_abort_excit(branch_pgen, t)

        do iOrb = start2 + 1, ende1 - 1
            call doubleUpdateStochastic(ilut, csf_i, iOrb, excitInfo, weights, negSwitches, &
                                        posSwitches, t, branch_pgen)
            ! check validity
            check_abort_excit(branch_pgen, t)

        end do

        ! then update weights and and to lowering semi-stop
        weights = weights%ptr

        call calcLoweringSemiStopStochastic(ilut, csf_i, excitInfo, weights, negSwitches, &
                                            posSwitches, t, branch_pgen)

        ! check validity
        check_abort_excit(branch_pgen, t)

        excitInfo%currentGen = excitInfo%lastGen

        do iOrb = ende1 + 1, ende2 - 1
            call singleStochasticUpdate(ilut, csf_i, iOrb, excitInfo, weights, posSwitches, &
                                        negSwitches, t, temp_pgen)
            branch_pgen = branch_pgen * temp_pgen
            ! check validity
            check_abort_excit(branch_pgen, t)

        end do

        ! and finally to end step
        call singleStochasticEnd(csf_i, excitInfo, t)

        ! if we do RDMs also store the x0 and x1 coupling coeffs
        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), &
                                            x0=extract_matrix_element(t, 1), &
                                            x1=extract_matrix_element(t, 2))
        end if

        ! for the additional contributing integrals:
        ! i have to consider that there might be a non-overlap excitation,
        ! which leads to the same excitation if there is no change in the
        ! stepvector in the overlap region!
        ! where i have to parts of the matrix elements, but which can be
        ! expressed in terms of the already calulated one i think..
        ! atleast in the case if R2L and L2R, as both the bottom and top
        ! parts are the same and only the semi-start ans semi-stop have to
        ! be modified..
        ! and also the x1 element has to be dismissed i guess..
        ! so i should change how i deal with the x1 elements and maybe
        ! keep them out here in these functions..
        switch = findFirstSwitch(ilut, t, start2 + 1, ende1)
        if (switch > 0) then
            integral = extract_matrix_element(t, 2) * (get_umat_el(ende1, start2, start1, ende2) + &
                                                       get_umat_el(start2, ende1, ende2, start1)) / 2.0_dp

            if (near_zero(integral)) then
                branch_pgen = 0.0_dp
                t = 0
            else
                call encode_matrix_element(t, 0.0_dp, 2)
                call encode_matrix_element(t, integral, 1)
            end if
        else
            integral = (-extract_matrix_element(t, 1) * (get_umat_el(start2, ende1, start1, ende2) + &
                                                     get_umat_el(ende1, start2, ende2, start1)) * 2.0_dp + (extract_matrix_element(t, 1) + &
                                                              extract_matrix_element(t, 2)) * (get_umat_el(ende1, start2, start1, ende2) + &
                                                                                        get_umat_el(start2, ende1, ende2, start1))) / 2.0_dp

            if (near_zero(integral)) then
                branch_pgen = 0.0_dp
                t = 0_n_int
            else
                call encode_matrix_element(t, 0.0_dp, 2)
                call encode_matrix_element(t, integral, 1)
            end if
        end if
    end subroutine calcDoubleL2R_stochastic

    subroutine calcDoubleL2R2L_stochastic(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 = "calcDoubleL2R2L_stochastic"

        integer :: iOrb, start2, ende1, ende2, start1, switch
        type(WeightObj_t) :: weights
        real(dp) :: temp_pgen
        HElement_t(dp) :: integral

        ! have to create this additional routine to more efficiently
        ! incorporate the integral contributions, since it is vastly different
        ! when involving mixed generators, but thats still todo!
        ASSERT(.not. isZero(ilut, excitInfo%fullStart))
        ASSERT(.not. isThree(ilut, excitInfo%secondStart))
        ASSERT(.not. isZero(ilut, excitInfo%firstEnd))
        ASSERT(.not. isThree(ilut, excitInfo%fullEnd))

        start1 = excitInfo%fullStart
        start2 = excitInfo%secondStart
        ende1 = excitInfo%firstEnd
        ende2 = excitInfo%fullEnd

        if (present(opt_weight)) then
            weights = opt_weight
        else
            ! : create correct weights:
            weights = init_fullDoubleWeight(csf_i, start2, ende1, ende2, negSwitches(start2), &
                                            negSwitches(ende1), posSwitches(start2), posSwitches(ende1), &
                                            csf_i%B_real(start2), csf_i%B_real(ende1))
        end if

        call createStochasticStart_single(ilut, csf_i, excitInfo, weights, posSwitches, &
                                          negSwitches, t, branch_pgen)

        ! check validity (defined in macros.h)
        check_abort_excit(branch_pgen, t)

        do iOrb = start1 + 1, start2 - 1
            call singleStochasticUpdate(ilut, csf_i, iOrb, excitInfo, weights, posSwitches, &
                                        negSwitches, t, temp_pgen)
            ! check validity
            branch_pgen = branch_pgen * temp_pgen

            check_abort_excit(branch_pgen, t)

        end do

        ! change weights... maybe need both single and double type weights
        ! then do lowering semi start
        weights = weights%ptr

        call calcRaisingSemiStartStochastic(ilut, csf_i, excitInfo, weights, negSwitches, &
                                            posSwitches, t, branch_pgen)

        ! check validity
        check_abort_excit(branch_pgen, t)

        do iOrb = start2 + 1, ende1 - 1
            call doubleUpdateStochastic(ilut, csf_i, iOrb, excitInfo, weights, negSwitches, &
                                        posSwitches, t, branch_pgen)
            ! check validity
            check_abort_excit(branch_pgen, t)

        end do

        ! then update weights and and to lowering semi-stop
        weights = weights%ptr

        call calcRaisingSemiStopStochastic(ilut, csf_i, excitInfo, weights, negSwitches, &
                                           posSwitches, t, branch_pgen)

        ! check validity
        check_abort_excit(branch_pgen, t)

        do iOrb = ende1 + 1, ende2 - 1
            call singleStochasticUpdate(ilut, csf_i, iOrb, excitInfo, weights, posSwitches, &
                                        negSwitches, t, temp_pgen)
            ! check validity
            branch_pgen = branch_pgen * temp_pgen

            check_abort_excit(branch_pgen, t)

        end do

        ! and finally to end step
        call singleStochasticEnd(csf_i, excitInfo, t)

        ! if we do RDMs also store the x0 and x1 coupling coeffs
        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), &
                                            x0=extract_matrix_element(t, 1), &
                                            x1=extract_matrix_element(t, 2))
        end if

        ! todo: think about the additional integral contributions and the
        ! relative sign of different order influences...
        ! and not sure yet if in this case i can use this function generally
        ! for R, RR, RR, R
        ! and L -> RL -> RL -> L
        ! since the integral/order influences are different maybe... todo!
        ! see above too!
        ! but yeah matrix elements definitly depends on excitation here.
        ! if there is no change in the overlap region, there is also the
        ! non-overlap contribution! otherwise not. maybe set some flag to
        ! indicate if such a stepvector change happened or not.

        ! update: on the matrix elements..
        ! i have to consider the non-overlap excitation, which can lead to
        ! the same excitation if there is no change in the stepvector
        ! in the overlap region
        ! the problem with the L2R2L and R2L2R funcitons is that the
        ! generator type changes for the non-overlap excitation so the
        ! matrix elements have to be changed more than in the R2L and L2R case

        ! if a switch happend -> also no additional contribution
        switch = findFirstSwitch(ilut, t, start2 + 1, ende1)

        if (switch > 0) then
            integral = extract_matrix_element(t, 2) * (get_umat_el(ende2, start2, start1, ende1) + &
                                                       get_umat_el(start2, ende2, ende1, start1)) / 2.0_dp

            if (near_zero(integral)) then
                branch_pgen = 0.0_dp
                t = 0_n_int
            else
                call encode_matrix_element(t, 0.0_dp, 2)
                call encode_matrix_element(t, integral, 1)
            end if

        else
            ! the no switch happened i have to get the additional contributions
            ! by the non-overlap version, but now the generator type at the
            ! end is different as the on-the fly calculated ...
            ! so the x0 matrix element changes by more (or even recalculate?)
            ! the bottom contribution, stays the same -sqrt(2) to cancel the
            ! -t
            ! the contribution at the semi-stop stays the same +t
            ! so those to get to -2.0_dp
            ! and bullshit that generator changes!! is also the same
            ! so it is exactly the same as in the R2L and L2R cases!! phew

            ! have also check if the non-overlap matrix elements is 0...
            ! this can happen unfortunately
            integral = (-extract_matrix_element(t, 1) * (get_umat_el(start2, ende2, start1, ende1) + &
                                                     get_umat_el(ende2, start2, ende1, start1)) * 2.0_dp + (extract_matrix_element(t, 1) + &
                                                              extract_matrix_element(t, 2)) * (get_umat_el(ende2, start2, start1, ende1) + &
                                                                                        get_umat_el(start2, ende2, ende1, start1))) / 2.0_dp

            if (near_zero(integral)) then
                branch_pgen = 0.0_dp
                t = 0_n_int
            else
                call encode_matrix_element(t, 0.0_dp, 2)
                call encode_matrix_element(t, integral, 1)
            end if
        end if

    end subroutine calcDoubleL2R2L_stochastic

    subroutine calcDoubleRaisingStochastic(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 = "calcDoubleRaisingStochastic"

        integer :: iOrb, start2, ende1, ende2, start1
        type(WeightObj_t) :: weights
        real(dp) :: temp_pgen
        HElement_t(dp) :: integral

        ASSERT(.not. isThree(ilut, excitInfo%fullStart))
        ASSERT(.not. isThree(ilut, excitInfo%secondStart))
        ASSERT(.not. isZero(ilut, excitInfo%firstEnd))
        ASSERT(.not. isZero(ilut, excitInfo%fullEnd))

        start1 = excitInfo%fullStart
        start2 = excitInfo%secondStart
        ende1 = excitInfo%firstEnd
        ende2 = excitInfo%fullEnd

        ! : create correct weights:
        if (present(opt_weight)) then
            weights = opt_weight
        else
            weights = init_fullDoubleWeight(csf_i, start2, ende1, ende2, negSwitches(start2), &
                                            negSwitches(ende1), posSwitches(start2), posSwitches(ende1), &
                                            csf_i%B_real(start2), csf_i%B_real(ende1))
        end if

        call createStochasticStart_single(ilut, csf_i, excitInfo, weights, posSwitches, &
                                          negSwitches, t, branch_pgen)

        ! check validity
        check_abort_excit(branch_pgen, t)

        do iOrb = start1 + 1, start2 - 1
            call singleStochasticUpdate(ilut, csf_i, iOrb, excitInfo, weights, posSwitches, &
                                        negSwitches, t, temp_pgen)
            ! check validity
            branch_pgen = branch_pgen * temp_pgen

            check_abort_excit(branch_pgen, t)

        end do

        ! change weights... maybe need both single and double type weights
        ! then do lowering semi start
        ! just point to the next weight:
        weights = weights%ptr

        call calcRaisingSemiStartStochastic(ilut, csf_i, excitInfo, weights, negSwitches, &
                                            posSwitches, t, branch_pgen)

        ! check validity
        check_abort_excit(branch_pgen, t)

        do iOrb = start2 + 1, ende1 - 1
            call doubleUpdateStochastic(ilut, csf_i, iOrb, excitInfo, weights, negSwitches, &
                                        posSwitches, t, branch_pgen)
            ! check validity

            check_abort_excit(branch_pgen, t)

        end do

        ! then update weights and and to lowering semi-stop
        weights = weights%ptr

        call calcRaisingSemiStopStochastic(ilut, csf_i, excitInfo, weights, negSwitches, &
                                           posSwitches, t, branch_pgen)

        ! check validity
        check_abort_excit(branch_pgen, t)

        do iOrb = ende1 + 1, ende2 - 1
            call singleStochasticUpdate(ilut, csf_i, iOrb, excitInfo, weights, posSwitches, &
                                        negSwitches, t, temp_pgen)
            ! check validity
            branch_pgen = branch_pgen * temp_pgen

            check_abort_excit(branch_pgen, t)

        end do

        ! and finally to end step
        call singleStochasticEnd(csf_i, excitInfo, t)

        ! if we do RDMs also store the x0 and x1 coupling coeffs
        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), &
                                            x0=extract_matrix_element(t, 1), &
                                            x1=extract_matrix_element(t, 2))
        end if

        ! todo: think about the additional integral contributions and the
        ! relative sign of different order influences...
        ! and not sure yet if in this case i can use this function generally
        ! for R, RR, RR, R
        ! and L -> RL -> RL -> L
        ! since the integral/order influences are different maybe... todo!

        ! think more about the sign influence in these kind of excitations..
        ! according to my notes, if i keep the x0 and x1 elements seperate
        ! until out here i can express it neatly in form of:
        ! x0(U1 + U2) + x1(U1 - U2)
        ! but not quite sure about that anymore... hopefully yes!
        integral = (extract_matrix_element(t, 1) * (get_umat_el(start1, start2, ende1, ende2) + &
                                                   get_umat_el(start2, start1, ende2, ende1) + get_umat_el(start1, start2, ende2, ende1) + &
                                                    get_umat_el(start2, start1, ende1, ende2)) + excitInfo%order * excitInfo%order1 * &
                    extract_matrix_element(t, 2) * ( &
                    get_umat_el(start1, start2, ende1, ende2) + get_umat_el(start2, start1, ende2, ende1) - &
                    get_umat_el(start1, start2, ende2, ende1) - get_umat_el(start2, start1, ende1, ende2))) / 2.0_dp

        if (near_zero(integral)) then
            branch_pgen = 0.0_dp
            t = 0_n_int
        else
            call encode_matrix_element(t, 0.0_dp, 2)
            call encode_matrix_element(t, integral, 1)
        end if

    end subroutine calcDoubleRaisingStochastic

    subroutine calcDoubleR2L2R_stochastic(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 = "calcDoubleR2L2R_stochastic"

        integer :: iOrb, switch
        type(WeightObj_t) :: weights
        real(dp) :: temp_pgen
        HElement_t(dp) :: integral

        associate (i => excitInfo%i, j => excitInfo%j, k => excitInfo%k, &
                   l => excitInfo%l, start1 => excitInfo%fullstart, &
                   start2 => excitInfo%secondStart, ende1 => excitInfo%firstEnd, &
                   ende2 => excitInfo%fullEnd, typ => excitInfo%typ)

            ASSERT(.not. isThree(ilut, start1))
            ASSERT(.not. isZero(ilut, start2))
            ASSERT(.not. isThree(ilut, ende1))
            ASSERT(.not. isZero(ilut, ende2))

            if (present(opt_weight)) then
                weights = opt_weight
            else
                ! : create correct weights:
                weights = init_fullDoubleWeight(csf_i, start2, ende1, ende2, negSwitches(start2), &
                                                negSwitches(ende1), posSwitches(start2), posSwitches(ende1), &
                                                csf_i%B_real(start2), csf_i%B_real(ende1))
            end if

            call createStochasticStart_single(ilut, csf_i, excitInfo, weights, posSwitches, &
                                              negSwitches, t, branch_pgen)

            ! check validity
            check_abort_excit(branch_pgen, t)

            do iOrb = start1 + 1, start2 - 1
                call singleStochasticUpdate(ilut, csf_i, iOrb, excitInfo, weights, posSwitches, &
                                            negSwitches, t, temp_pgen)
                ! check validity
                branch_pgen = branch_pgen * temp_pgen

                check_abort_excit(branch_pgen, t)

            end do

            ! change weights... maybe need both single and double type weights
            ! then do lowering semi start
            weights = weights%ptr

            call calcLoweringSemiStartStochastic(ilut, csf_i, excitInfo, weights, negSwitches, &
                                                 posSwitches, t, branch_pgen)

            ! check validity
            check_abort_excit(branch_pgen, t)

            do iOrb = start2 + 1, ende1 - 1
                call doubleUpdateStochastic(ilut, csf_i, iOrb, excitInfo, weights, negSwitches, &
                                            posSwitches, t, branch_pgen)
                ! check validity

                check_abort_excit(branch_pgen, t)

            end do

            ! then update weights and and to lowering semi-stop
            weights = weights%ptr

            call calcLoweringSemiStopStochastic(ilut, csf_i, excitInfo, weights, negSwitches, &
                                                posSwitches, t, branch_pgen)

            ! check validity
            check_abort_excit(branch_pgen, t)

            do iOrb = ende1 + 1, ende2 - 1
                call singleStochasticUpdate(ilut, csf_i, iOrb, excitInfo, weights, posSwitches, &
                                            negSwitches, t, temp_pgen)

                branch_pgen = branch_pgen * temp_pgen
                ! check validity
                check_abort_excit(branch_pgen, t)

            end do

            ! and finally to end step
            call singleStochasticEnd(csf_i, excitInfo, t)

            ! if we do RDMs also store the x0 and x1 coupling coeffs
            if (tFillingStochRDMOnFly) then
                call encode_stochastic_rdm_info(GugaBits, t, rdm_ind= &
                        contract_2_rdm_ind(i, j, k, l, excit_lvl=2, excit_typ=typ), &
                                                x0=extract_matrix_element(t, 1), &
                                                x1=extract_matrix_element(t, 2))
            end if

            ! todo: think about the additional integral contributions and the
            ! relative sign of different order influences...
            ! and not sure yet if in this case i can use this function generally
            ! for L -> LL -> LL -> L
            ! and R -> RL -> RL -> R
            ! since the integral/order influences are different maybe... todo!

            ! update: on the matrix elements..
            ! i have to consider the non-overlap excitation, which can lead to
            ! the same excitation if there is no change in the stepvector
            ! in the overlap region
            ! the problem with the L2R2L and R2L2R funcitons is that the
            ! generator type changes for the non-overlap excitation so the
            ! matrix elements have to be changed more than in the R2L and L2R case

            switch = findFirstSwitch(ilut, t, start2 + 1, ende1)
            if (switch > 0) then
                integral = extract_matrix_element(t, 2) * (get_umat_el(start1, ende1, ende2, start2) + &
                                                           get_umat_el(ende1, start1, start2, ende2)) / 2.0_dp

                if (near_zero(integral)) then
                    branch_pgen = 0.0_dp
                    t = 0_n_int
                else
                    call encode_matrix_element(t, 0.0_dp, 2)
                    call encode_matrix_element(t, integral, 1)
                end if
            else
                integral = (-extract_matrix_element(t, 1) * (get_umat_el(start1, ende1, start2, ende2) + &
                                                     get_umat_el(ende1, start1, ende2, start2)) * 2.0_dp + (extract_matrix_element(t, 1) + &
                                                              extract_matrix_element(t, 2)) * (get_umat_el(start1, ende1, ende2, start2) + &
                                                                                        get_umat_el(ende1, start1, start2, ende2))) / 2.0_dp

                if (near_zero(integral)) then
                    branch_pgen = 0.0_dp
                    t = 0_n_int
                else
                    call encode_matrix_element(t, 0.0_dp, 2)
                    call encode_matrix_element(t, integral, 1)
                end if
            end if

        end associate

    end subroutine calcDoubleR2L2R_stochastic

    subroutine calcDoubleLoweringStochastic(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 = "calcDoubleLoweringStochastic"

        integer :: iOrb, start2, ende1, ende2, start1
        type(WeightObj_t) :: weights
        real(dp) :: temp_pgen
        HElement_t(dp) :: integral

        ASSERT(.not. isZero(ilut, excitInfo%fullStart))
        ASSERT(.not. isZero(ilut, excitInfo%secondStart))
        ASSERT(.not. isThree(ilut, excitInfo%firstEnd))
        ASSERT(.not. isThree(ilut, excitInfo%fullEnd))

        start1 = excitInfo%fullStart
        start2 = excitInfo%secondStart
        ende1 = excitInfo%firstEnd
        ende2 = excitInfo%fullEnd

        ! : create correct weights:
        if (present(opt_weight)) then
            weights = opt_weight
        else
            weights = init_fullDoubleWeight(csf_i, start2, ende1, ende2, negSwitches(start2), &
                                            negSwitches(ende1), posSwitches(start2), posSwitches(ende1), &
                                            csf_i%B_real(start2), csf_i%B_real(ende1))
        end if

        call createStochasticStart_single(ilut, csf_i, excitInfo, weights, posSwitches, &
                                          negSwitches, t, branch_pgen)

        ! check validity
        check_abort_excit(branch_pgen, t)

        do iOrb = start1 + 1, start2 - 1
            call singleStochasticUpdate(ilut, csf_i, iOrb, excitInfo, weights, posSwitches, &
                                        negSwitches, t, temp_pgen)
            ! check validity
            branch_pgen = branch_pgen * temp_pgen

            check_abort_excit(branch_pgen, t)
        end do

        ! change weights... maybe need both single and double type weights
        ! then do lowering semi start
        ! can i just do:
        weights = weights%ptr

        ! branch_pgen gets update insde the routine!
        call calcLoweringSemiStartStochastic(ilut, csf_i, excitInfo, weights, negSwitches, &
                                             posSwitches, t, branch_pgen)

        ! check validity
        check_abort_excit(branch_pgen, t)

        do iOrb = start2 + 1, ende1 - 1
            ! branch_pgen gets updated inside update routine
            call doubleUpdateStochastic(ilut, csf_i, iOrb, excitInfo, weights, negSwitches, &
                                        posSwitches, t, branch_pgen)
            ! here only need to have probweight, since i cant only check x1 element
            ! check validity

            check_abort_excit(branch_pgen, t)
        end do

        ! then update weights and and to lowering semi-stop
        weights = weights%ptr

        ! branch_pgen gets updated inside funciton
        call calcLoweringSemiStopStochastic(ilut, csf_i, excitInfo, weights, negSwitches, &
                                            posSwitches, t, branch_pgen)

        ! check validity
        check_abort_excit(branch_pgen, t)

        do iOrb = ende1 + 1, ende2 - 1
            call singleStochasticUpdate(ilut, csf_i, iOrb, excitInfo, weights, posSwitches, &
                                        negSwitches, t, temp_pgen)
            ! check validity
            branch_pgen = branch_pgen * temp_pgen

            check_abort_excit(branch_pgen, t)
        end do

        ! and finally to end step
        call singleStochasticEnd(csf_i, excitInfo, t)

        ! if we do RDMs also store the x0 and x1 coupling coeffs
        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), &
                                            x0=extract_matrix_element(t, 1), &
                                            x1=extract_matrix_element(t, 2))
        end if

        ! todo: think about the additional integral contributions and the
        ! relative sign of different order influences...
        ! and not sure yet if in this case i can use this function generally
        ! for L -> LL -> LL -> L
        ! and R -> RL -> RL -> R
        ! since the integral/order influences are different maybe... todo!

        ! think more about the sign influence in these kind of excitations..
        ! according to my notes, if i keep the x0 and x1 elements seperate
        ! until out here i can express it neatly in form of:
        ! x0(U1 + U2) + x1(U1 - U2)
        ! but not quite sure about that anymore... hopefully yes!
        ! try that for now!

        integral = (extract_matrix_element(t, 1) * (get_umat_el(ende1, ende2, start1, start2) + &
                                                   get_umat_el(ende2, ende1, start2, start1) + get_umat_el(ende2, ende1, start1, start2) + &
                                                    get_umat_el(ende1, ende2, start2, start1)) + excitInfo%order * excitInfo%order1 * &
                    extract_matrix_element(t, 2) * ( &
                    get_umat_el(ende1, ende2, start1, start2) + get_umat_el(ende2, ende1, start2, start1) - &
                    get_umat_el(ende2, ende1, start1, start2) - get_umat_el(ende1, ende2, start2, start1))) / 2.0_dp

        if (near_zero(integral)) then
            branch_pgen = 0.0_dp
            t = 0_n_int
        else
            call encode_matrix_element(t, 0.0_dp, 2)
            call encode_matrix_element(t, integral, 1)
        end if

    end subroutine calcDoubleLoweringStochastic

    subroutine calcFullStopL2R_stochastic(ilut, csf_i, excitInfo, t, 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(inout) :: excitInfo
        integer(n_int), intent(out) :: t(0:nifguga)
        real(dp), intent(out) :: pgen
        real(dp), intent(in) :: posSwitches(nSpatOrbs), negSwitches(nSpatOrbs)
        type(WeightObj_t), intent(in), optional :: opt_weight

        type(WeightObj_t) :: weights
        integer :: st, se, en, i, j, k, l, elecInd, holeInd
        real(dp) :: branch_pgen, &
                    temp_pgen, rdm_mat, p_orig, orb_pgen
        HElement_t(dp) :: integral

        st = excitInfo%fullStart
        se = excitInfo%secondStart
        en = excitInfo%fullEnd

        ! init weights
        if (present(opt_weight)) then
            weights = opt_weight
        else
            if (t_approx_exchange .or. (t_approx_exchange_noninits .and. (.not. is_init_guga))) then
                ! the weights should be the only necessary change to force
                ! a switch at the end, as the other branches get 0 weight..
                weights = init_forced_end_semistart_weight(csf_i, se, en, negSwitches(se), &
                                                           posSwitches(se), csf_i%B_real(se))

            else
                weights = init_semiStartWeight(csf_i, se, en, negSwitches(se), &
                                               posSwitches(se), csf_i%B_real(se))
            end if
        end if

        ! create st
        call createStochasticStart_single(ilut, csf_i, excitInfo, weights, posSwitches, &
                                          negSwitches, t, branch_pgen)

        ! in case of early access the pgen should be set to 0
        pgen = 0.0_dp

        ! check validity
        check_abort_excit(branch_pgen, t)

        do i = st + 1, se - 1
            call singleStochasticUpdate(ilut, csf_i, i, excitInfo, weights, posSwitches, &
                                        negSwitches, t, temp_pgen)
            branch_pgen = branch_pgen * temp_pgen
            ! check validity

            check_abort_excit(branch_pgen, t)
        end do

        ! do the specific se-st
        ! try the new reusing of the weights object..
        weights = weights%ptr

        call calcRaisingSemiStartStochastic(ilut, csf_i, excitInfo, weights, negSwitches, &
                                            posSwitches, t, branch_pgen)

        ! check validity
        check_abort_excit(branch_pgen, t)

        ! do the specific double update to ensure a switch
        ! although switch can also happen at end only...
        ! actually that would be, in the full-stop case, temporary measure...
        ! but would unjust favor certain types of excitations..
        do i = se + 1, en - 1
            call doubleUpdateStochastic(ilut, csf_i, i, excitInfo, &
                                        weights, negSwitches, posSwitches, t, branch_pgen)

            if (near_zero(extract_matrix_element(t, 2)) .or. near_zero(branch_pgen)) then
                t = 0_n_int
                return
            end if
        end do

        call mixedFullStopStochastic(ilut, csf_i, excitInfo, t)

        ! check if matrix element is non-zero and if a switch happened
        if (.not. near_zero(extract_matrix_element(t, 1))) then
            t = 0_n_int
            branch_pgen = 0.0_dp
            return
        end if

        if (t_approx_exchange .or. (t_approx_exchange_noninits .and. (.not. is_init_guga))) then
            if (getDeltaB(t) == 0) then
                t = 0_n_int
                branch_pgen = 0.0_dp
                return
            end if
        end if

        if (near_zero(extract_matrix_element(t, 2))) then
            branch_pgen = 0.0_dp
            t = 0_n_int
            return
        end if

        ! if we do RDMs also store the x0 and x1 coupling coeffs
        ! and I need to do it before the routines below since excitInfo
        ! gets changed there
        if (tFillingStochRDMOnFly) then
            ! i need to unbias against the total pgen later on in the
            ! RDM sampling otherwise the rdm-bias factor is not correct!
            ! encode the necessary information in the rdm-matele!
            i = excitInfo%i
            j = excitInfo%j
            k = excitInfo%k
            l = excitInfo%l
            elecInd = st
            holeInd = se
            rdm_mat = extract_matrix_element(t, 2)
            call calc_orbital_pgen_contrib_end(&
                    csf_i, [2 * elecInd, 2 * en], holeInd, orb_pgen)
            p_orig = orb_pgen * branch_pgen / real(ElecPairs, dp)
            if (csf_i%stepvector(elecInd) == 3) p_orig = p_orig * 2.0_dp
        end if

        call encode_matrix_element(t, extract_matrix_element(t, 2), 1)

        ! actually I should provide a new routine, which "just"
        ! calculates the matrix element contribution and not
        ! the modified pgen, as the spatial orbitals are now fixed
        ! this could be done in the initialisation, where i just
        ! point to a new function, which only calculates the
        ! matrix element contribution

        global_excitInfo = excitInfo

        if (t_approx_exchange .or. (t_approx_exchange_noninits .and. (.not. is_init_guga))) then
            call calc_mixed_end_contr_approx(t, csf_i, excitInfo, integral)
            pgen = branch_pgen

        else
            call calc_mixed_end_contr_integral(ilut, csf_i, t, excitInfo, &
                integral)
            call calc_mixed_end_contr_pgen(ilut, csf_i, t, excitInfo, &
                branch_pgen, pgen)

!             call calc_mixed_end_contr_sym(ilut, csf_i, t, excitInfo, branch_pgen, pgen, &
!                                           integral)
        end if

        if (tFillingStochRDMOnFly) then
            if (.not. near_zero(p_orig)) then
                call encode_stochastic_rdm_info(GugaBits, t, rdm_ind= &
                        contract_2_rdm_ind(i, j, k, l, excit_lvl=2, &
                                       excit_typ=excitInfo%typ), x0=0.0_dp, &
                                                x1=rdm_mat * pgen / p_orig)
            end if
        end if

        call encode_matrix_element(t, 0.0_dp, 2)
        call update_matrix_element(t, (get_umat_el(en, se, st, en) + &
                                       get_umat_el(se, en, en, st)) / 2.0_dp + integral, 1)

    end subroutine calcFullStopL2R_stochastic

    subroutine calc_mixed_end_contr_approx(t, csf_i, excitInfo, integral)
        ! for the approx. mixed end contribution i "just" need to
        ! calculate the correct matrix element influences
        integer(n_int), intent(in) :: t(0:nifguga)
        type(CSF_Info_t), intent(in) :: csf_i
        type(ExcitationInformation_t), intent(in) :: excitInfo
        HElement_t(dp), intent(out) :: integral
        character(*), parameter :: this_routine = "calc_mixed_end_contr_approx"

        integer :: st, se, en, elecInd, holeInd, step, sw, i
        real(dp) :: top_cont, mat_ele, stay_mat, end_mat
        logical :: above_flag

        ! do as much stuff as possible beforehand
        st = excitInfo%fullStart
        se = excitInfo%secondStart
        en = excitInfo%fullEnd
        if (excitInfo%typ == excit_type%fullstop_L_to_R) then
            elecInd = st
            holeInd = se
        else if (excitInfo%typ == excit_type%fullstop_R_to_L) then
            elecInd = se
            holeInd = st
        else
            call stop_all(this_routine, "should not be here!")
        end if

        integral = h_cast(0.0_dp)

        step = csf_i%stepvector(en)

        ! i am sure the last switch happens at the full-stop!
        sw = en

        if (en < nSpatOrbs) then
            select case (step)
            case (1)
                if (isOne(t, en)) then
                    top_cont = -Root2 * sqrt((csf_i%B_real(en) + 2.0_dp) / &
                                             csf_i%B_real(en))

                else
                    top_cont = -Root2 / sqrt(csf_i%B_real(en) * (csf_i%B_real(en) + 2.0_dp))

                end if
            case (2)
                if (isOne(t, en)) then
                    top_cont = -Root2 / sqrt(csf_i%B_real(en) * (csf_i%B_real(en) + 2.0_dp))

                else
                    top_cont = Root2 * sqrt(csf_i%B_real(en) / &
                                            (csf_i%B_real(en) + 2.0_dp))
                end if

            case default
                call stop_all(this_routine, "wrong stepvalues!")

            end select

            if (.not. near_zero(top_cont)) then

                above_flag = .false.
                mat_ele = 1.0_dp

                do i = en + 1, nSpatOrbs
                    if (csf_i%Occ_int(i) /= 1) cycle

                    ! then check if thats the last step
                    if (csf_i%stepvector(i) == 2 .and. csf_i%B_int(i) == 0) then
                        above_flag = .true.
                    end if

                    ! in the other routine i check if the orbital pgen
                    ! is 0 for the above orbitals.. do I need to do that
                    ! also here?? or is this implicit if the matrix
                    ! element will be 0??

                    step = csf_i%stepvector(i)

                    call getDoubleMatrixElement(step, step, 0, gen_type%L, gen_type%R, csf_i%B_real(i), &
                                                1.0_dp, x1_element=stay_mat)

                    call getMixedFullStop(step, step, 0, csf_i%B_real(i), &
                                          x1_element=end_mat)

                    ! this check should never be true, but just to be sure
                    if (near_zero(stay_mat)) above_flag = .true.

                    if (.not. near_zero(end_mat)) then
                        integral = integral + end_mat * mat_ele * &
                                   (get_umat_el(i, holeInd, elecInd, i) + &
                                    get_umat_el(holeInd, i, i, elecInd)) / 2.0_dp
                    end if

                    if (above_flag) exit

                    ! otherwise update your running matrix element vars
                    mat_ele = mat_ele * stay_mat

                end do

                integral = integral * top_cont
            end if
        end if

    end subroutine calc_mixed_end_contr_approx

    subroutine calcFullStopR2L_stochastic(ilut, csf_i, excitInfo, t, 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(inout) :: excitInfo
        integer(n_int), intent(out) :: t(0:nifguga)
        real(dp), intent(out) :: pgen
        real(dp), intent(in) :: posSwitches(nSpatOrbs), negSwitches(nSpatOrbs)
        type(WeightObj_t), intent(in), optional :: opt_weight

        type(WeightObj_t) :: weights
        integer :: st, se, en, i, j, k, l, elecInd, holeInd
        real(dp) :: branch_pgen, &
                    temp_pgen, p_orig, rdm_mat, orb_pgen
        HElement_t(dp) :: integral

        st = excitInfo%fullStart
        se = excitInfo%secondStart
        en = excitInfo%fullEnd

        ! init weights
        if (present(opt_weight)) then
            weights = opt_weight
        else
            if (t_approx_exchange .or. (t_approx_exchange_noninits .and. (.not. is_init_guga))) then
                weights = init_forced_end_semistart_weight(csf_i, se, en, negSwitches(se), &
                                                           posSwitches(se), csf_i%B_real(se))
            else
                weights = init_semiStartWeight(csf_i, se, en, negSwitches(se), &
                                               posSwitches(se), csf_i%B_real(se))
            end if
        end if

        ! create start
        call createStochasticStart_single(ilut, csf_i, excitInfo, weights, posSwitches, &
                                          negSwitches, t, branch_pgen)

        ! in case of early exit pgen should be set to 0
        pgen = 0.0_dp
        ! check validity

        check_abort_excit(branch_pgen, t)

        do i = st + 1, se - 1
            call singleStochasticUpdate(ilut, csf_i, i, excitInfo, weights, posSwitches, &
                                        negSwitches, t, temp_pgen)
            branch_pgen = branch_pgen * temp_pgen
            ! check validity

            check_abort_excit(branch_pgen, t)
        end do

        ! do the specific semi-start
        weights = weights%ptr

        call calcLoweringSemiStartStochastic(ilut, csf_i, excitInfo, weights, negSwitches, &
                                             posSwitches, t, branch_pgen)

        ! check validity
        check_abort_excit(branch_pgen, t)

        ! do the specific double update to ensure a switch
        ! although switch can also happen at end only...
        ! actually that would be, in the full-stop case, temporary measure...
        ! but would unjust favor certain types of excitations..
        do i = se + 1, en - 1
            call doubleUpdateStochastic(ilut, csf_i, i, excitInfo, &
                                        weights, negSwitches, posSwitches, t, branch_pgen)
            ! also here there has to be a switch at some point so check x1
            if (near_zero(extract_matrix_element(t, 2)) .or. near_zero(branch_pgen)) then
                t = 0_n_int
                return
            end if
        end do

        call mixedFullStopStochastic(ilut, csf_i, excitInfo, t)

        ! check if matrix element is non-zero and if a switch happened
        if (.not. near_zero(extract_matrix_element(t, 1))) then
            t = 0_n_int
            return
        end if
        if (t_approx_exchange .or. (t_approx_exchange_noninits .and. (.not. is_init_guga))) then
            ! make it crude for now, that we only check if the delta B value
            ! is non-zero at the end, otherwise abort this spawn..
            if (getDeltaB(t) == 0) then
                t = 0_n_int
                return
            end if
        end if

        if (near_zero(extract_matrix_element(t, 2))) then
            t = 0_n_int
            return
        end if

        ! if we do RDMs also store the x0 and x1 coupling coeffs
        ! and I need to do it before the routines below since excitInfo
        ! gets changed there
        if (tFillingStochRDMOnFly) then
            ! i need to unbias against the total pgen later on in the
            ! RDM sampling otherwise the rdm-bias factor is not correct!
            ! encode the necessary information in the rdm-matele!
            i = excitInfo%i
            j = excitInfo%j
            k = excitInfo%k
            l = excitInfo%l
            elecInd = se
            holeInd = st
            rdm_mat = extract_matrix_element(t, 2)
            call calc_orbital_pgen_contrib_end(&
                    csf_i, [2 * elecInd, 2 * en], holeInd, orb_pgen)
            p_orig = orb_pgen * branch_pgen / real(ElecPairs, dp)
            if (csf_i%stepvector(elecInd) == 3) p_orig = p_orig * 2.0_dp

        end if

        ! the x1-element is still encoded in the second entry..
        ! move it to the first elemen
        call encode_matrix_element(t, extract_matrix_element(t, 2), 1)

        global_excitInfo = excitInfo

        if (t_approx_exchange .or. (t_approx_exchange_noninits .and. (.not. is_init_guga))) then

            call calc_mixed_end_contr_approx(t, csf_i, excitInfo, integral)
            pgen = branch_pgen

        else
            call calc_mixed_end_contr_integral(ilut, csf_i, t, excitInfo, &
                integral)
            call calc_mixed_end_contr_pgen(ilut, csf_i, t, excitInfo, &
                branch_pgen, pgen)

!             call calc_mixed_end_contr_sym(ilut, csf_i, t, excitInfo, branch_pgen, pgen, integral)
        end if

        if (tFillingStochRDMOnFly) then
            if (.not. near_zero(p_orig)) then
                call encode_stochastic_rdm_info(GugaBits, t, rdm_ind= &
                        contract_2_rdm_ind(i, j, k, l, excit_lvl=2, &
                                           excit_typ=excitInfo%typ), x0=0.0_dp, &
                                                x1=rdm_mat * pgen / p_orig)
            end if
        end if

        call encode_matrix_element(t, 0.0_dp, 2)
        call update_matrix_element(t, (get_umat_el(en, st, se, en) + &
                                       get_umat_el(st, en, en, se)) / 2.0_dp + integral, 1)

    end subroutine calcFullStopR2L_stochastic

    subroutine doubleUpdateStochastic(ilut, csf_i, s, excitInfo, weights, negSwitches, &
                                      posSwitches, t, probWeight)
        integer(n_int), intent(in) :: ilut(0:nifguga)
        type(CSF_Info_t), intent(in) :: csf_i
        integer, intent(in) :: s
        type(ExcitationInformation_t), intent(in) :: excitInfo
        type(WeightObj_t), intent(in) :: weights
        real(dp), intent(in) :: negSwitches(nSpatOrbs), posSwitches(nSpatOrbs)
        integer(n_int), intent(inout) :: t(0:nifguga)
        real(dp), intent(inout) :: probWeight
        character(*), parameter :: this_routine = "doubleUpdateStochastic"

        real(dp) :: minusWeight, plusWeight, zeroWeight
        integer :: gen1, gen2, deltaB
        real(dp) :: bVal, tempWeight_0, tempWeight_1, order

        ASSERT(isProperCSF_ilut(ilut))
        ASSERT(s > 0 .and. s <= nSpatOrbs)

        if (csf_i%Occ_int(s) /= 1) then
            ! no change in stepvector or matrix element in this case
            return
        end if

        ! and for more readibility extract certain values:
        gen1 = excitInfo%gen1
        gen2 = excitInfo%gen2
        bVal = csf_i%B_real(s)
        ! stupid! only need at order at semistarts and semistops and not for
        ! the overlap region
        order = 1.0_dp

        deltaB = getDeltaB(t)

        ! new idea: make the combination stepvalue + deltaB !
        ! this give me 6 distinct integer quantities which i can choose
        ! from in a select case statement!

        select case (csf_i%stepvector(s) + deltaB)
            ! depending on the deltaB value different possibs
        case (3)
            ! d=1 + b=2 = 3
            ! only staying
            call getDoubleMatrixElement(1, 1, deltaB, gen1, gen2, bVal, &
                                        order, x1_element=tempWeight_1)

            tempWeight_0 = 0.0_dp

        case (-1)
            ! d=1 + b=-2 = -1
            ! both 0 and -2 branches are possible
            minusWeight = weights%proc%minus(negSwitches(s), &
                                             bVal, weights%dat)
            zeroWeight = weights%proc%zero(negSwitches(s), &
                                           posSwitches(s), bVal, weights%dat)

            if (near_zero(minusWeight + zeroWeight)) then
                probWeight = 0.0_dp
                t = 0
                return
            end if

            minusWeight = calcStayingProb(minusWeight, zeroWeight, bVal)

            if (genrand_real2_dSFMT() < minusWeight) then
                ! stay on -2 branch

                call getDoubleMatrixElement(1, 1, deltaB, gen1, gen2, &
                                            bVal, order, x1_element=tempWeight_1)

                tempWeight_0 = 0.0_dp

                probWeight = probWeight * minusWeight

            else
                ! switch to 0 branch
                ! 1 -> 2
                clr_orb(t, 2 * s - 1)
                set_orb(t, 2 * s)

                call setDeltaB(0, t)

                call getDoubleMatrixElement(2, 1, deltaB, gen1, gen2, &
                                            bVal, order, x1_element=tempWeight_1)

                tempWeight_0 = 0.0_dp

                probWeight = probWeight * (1.0_dp - minusWeight)

            end if

        case (1)
            ! d=1 + b=0 = 1
            ! 0 branch arrives -> have to check b value
            ! hm... is this a bug?? if d = 1, b cant be 0, its atleast
            ! 1.. and then a switch to +2 cant happen..
            ! but that should have been dealt with the weights below
            ! probably.. so thats why it didnt matter probably..
            if (csf_i%B_int(s) == 1) then
                ! only staying branch
                call getDoubleMatrixElement(1, 1, deltaB, gen1, gen2, bVal, &
                                            order, tempWeight_0, tempWeight_1)

            else
                ! both 0 and +2 branch are possible
                plusWeight = weights%proc%plus(posSwitches(s), &
                                               bVal, weights%dat)
                zeroWeight = weights%proc%zero(negSwitches(s), &
                                               posSwitches(s), bVal, weights%dat)

                if (near_zero(plusWeight + zeroWeight)) then
                    probWeight = 0.0_dp
                    t = 0
                end if

                zeroWeight = calcStayingProb(zeroWeight, plusWeight, bVal)

                if (genrand_real2_dSFMT() < zeroWeight) then
                    ! stay on 0 branch
                    call getDoubleMatrixElement(1, 1, deltaB, gen1, gen2, &
                                                bVal, order, tempWeight_0, tempWeight_1)

                    probWeight = probWeight * zeroWeight

                else
                    ! switch to +2 branch
                    ! 1 -> 2
                    clr_orb(t, 2 * s - 1)
                    set_orb(t, 2 * s)

                    call setDeltaB(2, t)

                    call getDoubleMatrixElement(2, 1, deltaB, gen1, gen2, &
                                                bVal, order, x1_element=tempWeight_1)

                    tempWeight_0 = 0.0_dp

                    probWeight = probWeight * (1.0_dp - zeroWeight)
                end if
            end if
        case (0)
            ! d=2 + b=-2 : 0
            ! only staying
            call getDoubleMatrixElement(2, 2, deltaB, gen1, gen2, bVal, &
                                        order, x1_element=tempWeight_1)

            tempWeight_0 = 0.0_dp

        case (2)
            ! d=2 + b=0 : 2
            ! always -2 and 0 branching possible
            minusWeight = weights%proc%minus(negSwitches(s), &
                                             bVal, weights%dat)
            zeroWeight = weights%proc%zero(negSwitches(s), &
                                           posSwitches(s), bVal, weights%dat)

            ! here should be the only case where b*w_0 + w_- = 0 ->
            ! have to avoid divison by 0 then
            ! but in this case only the stayin on 0 branch is valid
            ! since it should be always > 0 and the -branch has 0 weight
            if (near_zero(bVal * zeroWeight + minusWeight)) then
                zeroWeight = 1.0_dp
            else
                zeroWeight = calcStayingProb(zeroWeight, minusWeight, bVal)
            end if

            if (near_zero(zeroWeight + minusWeight)) then
                probWeight = 0.0_dp
                t = 0
                return
            end if

            if (genrand_real2_dSFMT() < zeroWeight) then
                ! stay on 0 branch
                call getDoubleMatrixElement(2, 2, deltaB, gen1, gen2, bVal, &
                                            order, tempWeight_0, tempWeight_1)

                probWeight = probWeight * zeroWeight

            else
                ! switch to -2 branch
                ! 2 -> 1
                set_orb(t, 2 * s - 1)
                clr_orb(t, 2 * s)

                call setDeltaB(-2, t)

                call getDoubleMatrixElement(1, 2, deltaB, gen1, gen2, &
                                            bVal, order, x1_element=tempWeight_1)

                tempWeight_0 = 0.0_dp

                probWeight = probWeight * (1.0_dp - zeroWeight)

            end if

        case (4)
            ! d=2 + b=2 : 4

            ! have to check b value if branching is possible
            if (csf_i%B_int(s) < 2) then

                ! only switch possible
                ! 2 -> 1
                set_orb(t, 2 * s - 1)
                clr_orb(t, 2 * s)

                call setDeltaB(0, t)

                call getDoubleMatrixElement(1, 2, deltaB, gen1, gen2, &
                                            bVal, order, x1_element=tempWeight_1)

                tempWeight_0 = 0.0_dp

            else
                ! both 0 and +2 branches possible
                plusWeight = weights%proc%plus(posSwitches(s), &
                                               bVal, weights%dat)
                zeroWeight = weights%proc%zero(negSwitches(s), &
                                               posSwitches(s), bVal, weights%dat)

                if (near_zero(plusWeight + zeroWeight)) then
                    probWeight = 0.0_dp
                    t = 0
                    return
                end if

                plusWeight = calcStayingProb(plusWeight, zeroWeight, bVal)

                if (genrand_real2_dSFMT() < plusWeight) then
                    ! stay on +2 branch
                    call getDoubleMatrixElement(2, 2, deltaB, gen1, gen2, &
                                                bVal, order, x1_element=tempWeight_1)

                    tempWeight_0 = 0.0_dp

                    probWeight = probWeight * plusWeight

                else
                    ! switch to 0 branch
                    ! 2 -> 1
                    set_orb(t, 2 * s - 1)
                    clr_orb(t, 2 * s)

                    call setDeltaB(0, t)
                    call getDoubleMatrixElement(1, 2, deltaB, gen1, gen2, &
                                                bVal, order, x1_element=tempWeight_1)

                    tempWeight_0 = 0.0_dp

                    probWeight = probWeight * (1.0_dp - plusWeight)

                end if
            end if
        end select

        call update_matrix_element(t, tempWeight_0, 1)
        call update_matrix_element(t, tempWeight_1, 2)

        if (near_zero(tempWeight_0) .and. near_zero(tempWeight_1)) then
            probWeight = 0.0_dp
            t = 0
            return
        end if

        if (t_trunc_guga_pgen .or. &
            (t_trunc_guga_pgen_noninits .and. .not. is_init_guga)) then
            if (probWeight < trunc_guga_pgen) then
                probWeight = 0.0_dp
                t = 0_n_int
            end if
        end if

        if (t_trunc_guga_matel) then
            if (abs(extract_matrix_element(t, 1)) < trunc_guga_matel .and. &
                abs(extract_matrix_element(t, 2)) < trunc_guga_matel) then
                probWeight = 0.0_dp
                t = 0_n_int
            end if
        end if
    end subroutine doubleUpdateStochastic

    subroutine mixedFullStopStochastic(ilut, csf_i, excitInfo, t)
        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(inout) :: t(0:nifguga)
        character(*), parameter :: this_routine = "mixedFullStopStochastic"

        integer :: ende, deltaB
        real(dp) :: bVal, tempWeight_0, tempWeight_1

        ASSERT(.not. isThree(ilut, excitInfo%fullEnd))
        ASSERT(.not. isZero(ilut, excitInfo%fullEnd))
        ! no 3 allowed at end or else it would be single-excitation-like

        ende = excitInfo%fullEnd
        bVal = csf_i%B_real(ende)

        deltaB = getDeltaB(t)

        ! NOTE: also remember for mixed full-stops there has to be a switch at
        ! some point of the double overlap region or else it is single-excitation
        ! like... so only consider x1 matrix element here..

        ! combine deltaB and stepvalue info here to reduce if statements
        ! asserts dont work anymore with new select case statements
        ! do it out here:
#ifdef DEBUG_
        if (csf_i%stepvector(ende) == 1) then
            ASSERT(deltaB /= 2)
        end if
        if (csf_i%stepvector(ende) == 2) then
            ASSERT(deltaB /= -2)
        end if
#endif

        select case (deltaB + csf_i%stepvector(ende))
        case (1)
            ! d=1 + b=0 : 1
            ! ! +2 branch not allowed here
            ! not sure if i can access only the x1 element down there..
            call getMixedFullStop(1, 1, 0, bVal, tempWeight_0, tempWeight_1)

        case (-1)
            ! d=1 + b=-2 : -1
            ! deltaB = -2
            ! switch 1 -> 2
            set_orb(t, 2 * ende)
            clr_orb(t, 2 * ende - 1)

            ! matrix element is 1 in this case
            tempWeight_1 = 1.0_dp
            tempWeight_0 = 0.0_dp

        case (2)
            ! d=2 + b=0 : 2

            call getMixedFullStop(2, 2, 0, bVal, tempWeight_0, tempWeight_1)

        case (4)
            ! d=2 + b=2 : 4
            ! deltab = 2
            ! switch 2 -> 1
            set_orb(t, 2 * ende - 1)
            clr_orb(t, 2 * ende)

            tempWeight_1 = 1.0_dp
            tempWeight_0 = 0.0_dp
        end select

        ! thats kind of stupid what i did here...
        ! check if x0-matrix element is non-zero which is an indicator that
        ! no switch happened in the double overlap region

        ! just for completion reasons still update the x0 matrix element here
        ! although it should be 0 anyway..
        call update_matrix_element(t, tempWeight_0, 1)
        call update_matrix_element(t, tempWeight_1, 2)

        ! todo... have to write some excitation cancellation function...
        ! to get back to the start of an excitation or somehow ensure that
        ! switch happens...
    end subroutine mixedFullStopStochastic

    subroutine calcRaisingSemiStartStochastic(ilut, csf_i, excitInfo, weights, negSwitches, &
                                              posSwitches, t, probWeight)
        integer(n_int), intent(in) :: ilut(0:nifguga)
        type(CSF_Info_t), intent(in) :: csf_i
        type(ExcitationInformation_t), intent(in) :: excitInfo
        type(WeightObj_t), intent(in) :: weights
        real(dp), intent(in) :: negSwitches(nSpatOrbs), posSwitches(nSpatOrbs)
        integer(n_int), intent(inout) :: t(0:nifguga)
        real(dp), intent(inout) :: probWeight
        character(*), parameter :: this_routine = "calcRaisingSemiStartStochastic"

        integer :: se, deltaB
        real(dp) :: tempWeight_0, tempWeight_1, minusWeight, &
                    plusWeight, zeroWeight, bVal

        ASSERT(.not. isThree(ilut, excitInfo%secondStart))

        ! i can be sure that there is no 3 or 0 at the fullEnd, or otherwise
        ! this would be single-excitation like.
        se = excitInfo%secondStart
        bVal = csf_i%B_real(se)

        deltaB = getDeltaB(t)

        ! do non-choosing possibs first
        ! why does this cause a segfault on compilation with gfortran??
        ! do some debugging:

        ! fix for gfortran compilation for some reasono
        ! i can probably fix it when i finally get to this point in
        ! test running

        select case (csf_i%stepvector(se))
        case (1)
            ! 1 -> 3
            set_orb(t, 2 * se)

            call setDeltaB(deltaB + 1, t)

            call getDoubleMatrixElement(3, 1, deltaB, excitInfo%gen1, &
                                        excitInfo%gen2, bVal, excitInfo%order, &
                                        tempWeight_0, tempWeight_1)

        case (2)
            ! 2 -> 3
            set_orb(t, 2 * se - 1)

            call setDeltaB(deltaB - 1, t)

            call getDoubleMatrixElement(3, 2, deltaB, excitInfo%gen1, &
                                        excitInfo%gen2, bVal, excitInfo%order, &
                                        tempWeight_0, tempWeight_1)

        case (0)
            ! 0:
            ! for arriving -1 branch branching is always possible
            if (deltaB == -1) then
                ! here the choice is between 0 and -2 branch
                minusWeight = weights%proc%minus(negSwitches(se), bVal, weights%dat)
                zeroWeight = weights%proc%zero(negSwitches(se), posSwitches(se), &
                                               bVal, weights%dat)

                if (near_zero(minusWeight + zeroWeight)) then
                    probWeight = 0.0_dp
                    t = 0
                    return
                end if

                ! cant directly cant assign it to probWeight
                zeroWeight = calcStartProb(zeroWeight, minusWeight)

                if (genrand_real2_dSFMT() < zeroWeight) then
                    ! go to 0 branch
                    ! 0 -> 2
                    set_orb(t, 2 * se)

                    call setDeltaB(0, t)

                    call getDoubleMatrixElement(2, 0, deltaB, excitInfo%gen1, &
                                                excitInfo%gen2, bVal, excitInfo%order, &
                                                tempWeight_0, tempWeight_1)

                    probWeight = probWeight * zeroWeight

                else
                    ! to to -2 branch
                    ! 0 -> 1
                    set_orb(t, 2 * se - 1)

                    call setDeltaB(-2, t)

                    call getDoubleMatrixElement(1, 0, deltaB, excitInfo%gen1, &
                                                excitInfo%gen2, bVal, excitInfo%order, &
                                                x1_element=tempWeight_1)

                    tempWeight_0 = 0.0_dp

                    probWeight = probWeight * (1.0_dp - zeroWeight)
                end if
            else
                ! +1 branch arriving -> have to check b values
                ! UPDATE: include b value check into probWeight calculation
                ! todo
                if (csf_i%B_int(se) < 2) then
                    ! only 0 branch possible
                    ! todo: in this forced cases due to the b value, have to
                    ! think about, how that influences the probWeight...
                    ! 0 -> 1
                    set_orb(t, 2 * se - 1)

                    call setDeltaB(0, t)

                    call getDoubleMatrixElement(1, 0, deltaB, excitInfo%gen1, &
                                                excitInfo%gen2, bVal, excitInfo%order, &
                                                tempWeight_0, tempWeight_1)

                else
                    ! need probs
                    plusWeight = weights%proc%plus(posSwitches(se), &
                                                   bVal, weights%dat)
                    zeroWeight = weights%proc%zero(negSwitches(se), posSwitches(se), &
                                                   bVal, weights%dat)

                    if (near_zero(plusWeight + zeroWeight)) then
                        probWeight = 0.0_dp
                        t = 0
                        return
                    end if

                    ! cant directly cant assign it to probWeight
                    zeroWeight = calcStartProb(zeroWeight, plusWeight)

                    if (genrand_real2_dSFMT() < zeroWeight) then
                        ! do 0 branch
                        ! 0 -> 1
                        set_orb(t, 2 * se - 1)

                        call setDeltaB(0, t)

                        call getDoubleMatrixElement(1, 0, deltaB, excitInfo%gen1, &
                                                    excitInfo%gen2, bVal, excitInfo%order, &
                                                    tempWeight_0, tempWeight_1)

                        probWeight = probWeight * zeroWeight

                    else
                        ! do +2 branch
                        ! 0 -> 2
                        set_orb(t, 2 * se)

                        call setDeltaB(2, t)

                        call getDoubleMatrixElement(2, 0, deltaB, excitInfo%gen1, &
                                                    excitInfo%gen2, bVal, excitInfo%order, &
                                                    x1_element=tempWeight_1)

                        tempWeight_0 = 0.0_dp

                        probWeight = probWeight * (1.0_dp - zeroWeight)
                    end if
                end if
            end if
        end select

        call encode_matrix_element(t, extract_matrix_element(t, 1) * tempWeight_1, 2)
        call update_matrix_element(t, tempWeight_0, 1)

        if (near_zero(tempWeight_0) .and. near_zero(tempWeight_1)) then
            probWeight = 0.0_dp
            t = 0
        end if

    end subroutine calcRaisingSemiStartStochastic

    subroutine calcLoweringSemiStartStochastic(ilut, csf_i, excitInfo, weights, negSwitches, &
                                               posSwitches, t, probWeight)
        integer(n_int), intent(in) :: ilut(0:nifguga)
        type(CSF_Info_t), intent(in) :: csf_i
        type(ExcitationInformation_t), intent(in) :: excitInfo
        type(WeightObj_t), intent(in) :: weights
        real(dp), intent(in) :: negSwitches(nSpatOrbs), posSwitches(nSpatOrbs)
        integer(n_int), intent(inout) :: t(0:nifguga)
        real(dp), intent(inout) :: probWeight
        character(*), parameter :: this_routine = "calcLoweringSemiStartStochastic"

        integer :: se, deltaB
        real(dp) :: tempWeight_0, tempWeight_1, minusWeight, &
                    plusWeight, zeroWeight, bVal

        ASSERT(.not. isZero(ilut, excitInfo%secondStart))

        ! i can be sure that there is no 3 or 0 at the fullEnd, or otherwise
        ! this would be single-excitation like.

        se = excitInfo%secondStart
        bVal = csf_i%B_real(se)

        deltaB = getDeltaB(t)

        ! do non-choosing possibs first
        ! why does this cause a segfault on compilation with gfortran??
        ! do some debugging:
        ! same gfortran compilex issue fix as above

        select case (csf_i%stepvector(se))
        case (1)
            ! 1 -> 0
            clr_orb(t, 2 * se - 1)

            call setDeltaB(deltaB + 1, t)

            call getDoubleMatrixElement(0, 1, deltaB, excitInfo%gen1, &
                                        excitInfo%gen2, bVal, excitInfo%order, &
                                        tempWeight_0, tempWeight_1)

        case (2)
            ! 2 -> 0
            clr_orb(t, 2 * se)

            call setDeltaB(deltaB - 1, t)

            call getDoubleMatrixElement(0, 2, deltaB, excitInfo%gen1, &
                                        excitInfo%gen2, bVal, excitInfo%order, &
                                        tempWeight_0, tempWeight_1)

        case (3)
            ! 3:
            ! for arriving -1 branch branching is always possible
            if (deltaB == -1) then
                ! here the choice is between 0 and -2 branch
                minusWeight = weights%proc%minus(negSwitches(se), &
                                                 bVal, weights%dat)
                zeroWeight = weights%proc%zero(negSwitches(se), posSwitches(se), &
                                               bVal, weights%dat)

                if (near_zero(minusWeight + zeroWeight)) then
                    probWeight = 0.0_dp
                    t = 0
                    return
                end if
                ! cant directly cant assign it to probWeight
                zeroWeight = calcStartProb(zeroWeight, minusWeight)

                if (genrand_real2_dSFMT() < zeroWeight) then
                    ! go to 0 branch
                    ! 3 -> 2
                    clr_orb(t, 2 * se - 1)

                    call setDeltaB(0, t)

                    call getDoubleMatrixElement(2, 3, deltaB, excitInfo%gen1, &
                                                excitInfo%gen2, bVal, excitInfo%order, &
                                                tempWeight_0, tempWeight_1)

                    probWeight = probWeight * zeroWeight

                else
                    ! to to -2 branch
                    ! 3 -> 1
                    clr_orb(t, 2 * se)

                    call setDeltaB(-2, t)

                    call getDoubleMatrixElement(1, 3, deltaB, excitInfo%gen1, &
                                                excitInfo%gen2, bVal, excitInfo%order, &
                                                x1_element=tempWeight_1)

                    tempWeight_0 = 0.0_dp

                    probWeight = probWeight * (1.0_dp - zeroWeight)
                end if
            else
                ! +1 branch arriving -> have to check b values
                if (csf_i%B_int(se) < 2) then
                    ! only 0 branch possible
                    ! todo: in this forced cases due to the b value, have to
                    ! think about, how that influences the probWeight...
                    ! 3 -> 1
                    clr_orb(t, 2 * se)

                    call setDeltaB(0, t)

                    call getDoubleMatrixElement(1, 3, deltaB, excitInfo%gen1, &
                                                excitInfo%gen2, bVal, excitInfo%order, &
                                                tempWeight_0, tempWeight_1)

                else
                    ! need probs
                    plusWeight = weights%proc%plus(posSwitches(se), &
                                                   bVal, weights%dat)
                    zeroWeight = weights%proc%zero(negSwitches(se), posSwitches(se), &
                                                   bVal, weights%dat)

                    if (near_zero(plusWeight + zeroWeight)) then
                        probWeight = 0.0_dp
                        t = 0
                        return
                    end if

                    ! cant directly cant assign it to probWeight
                    zeroWeight = calcStartProb(zeroWeight, plusWeight)

                    if (genrand_real2_dSFMT() < zeroWeight) then
                        ! do 0 branch
                        ! 3 -> 1
                        clr_orb(t, 2 * se)

                        call setDeltaB(0, t)

                        call getDoubleMatrixElement(1, 3, deltaB, excitInfo%gen1, &
                                                    excitInfo%gen2, bVal, excitInfo%order, &
                                                    tempWeight_0, tempWeight_1)

                        probWeight = probWeight * zeroWeight

                    else
                        ! do +2 branch
                        ! 3 -> 2
                        clr_orb(t, 2 * se - 1)

                        call setDeltaB(2, t)

                        call getDoubleMatrixElement(2, 3, deltaB, excitInfo%gen1, &
                                                    excitInfo%gen2, bVal, excitInfo%order, &
                                                    x1_element=tempWeight_1)

                        tempWeight_0 = 0.0_dp

                        probWeight = probWeight * (1.0_dp - zeroWeight)
                    end if
                end if
            end if
        end select

        call encode_matrix_element(t, extract_matrix_element(t, 1) * tempWeight_1, 2)
        call update_matrix_element(t, tempWeight_0, 1)

        if (near_zero(tempWeight_0) .and. near_zero(tempWeight_1)) then
            probWeight = 0.0_dp
            t = 0
        end if

    end subroutine calcLoweringSemiStartStochastic

    subroutine calcRaisingSemiStopStochastic(ilut, csf_i, excitInfo, weights, negSwitches, &
                                             posSwitches, t, probWeight)
        integer(n_int), intent(in) :: ilut(0:nifguga)
        type(CSF_Info_t), intent(in) :: csf_i
        type(ExcitationInformation_t), intent(in) :: excitInfo
        type(WeightObj_t), intent(in) :: weights
        real(dp), intent(in) :: negSwitches(nSpatOrbs), posSwitches(nSpatOrbs)
        integer(n_int), intent(inout) :: t(0:nifguga)
        real(dp), intent(inout) :: probWeight
        character(*), parameter :: this_routine = "calcRaisingSemiStopStochastic"

        integer :: semi, deltaB
        real(dp) :: tempWeight_0, tempWeight_1, minusWeight, &
                    plusWeight, bVal

        ASSERT(.not. isZero(ilut, excitInfo%firstEnd))

        semi = excitInfo%firstEnd
        bVal = csf_i%B_real(semi)

        ! in the stochastic case i am sure that at there is no 3 or 0 at the
        ! full start... so definetly all deltaB branches can arrive here
        ! first deal with the forced ones
        deltaB = getDeltaB(t)

        select case (csf_i%stepvector(semi))
        case (1)
            if (t_approx_exchange .or. (t_approx_exchange_noninits .and. (.not. is_init_guga))) then
                if (getDeltaB(t) == 2) then
                    t = 0_n_int
                    probWeight = 0.0_dp
                    return
                end if
            end if

            ASSERT(getDeltaB(t) /= 2)

            ! do the 1 -> 0 switch
            clr_orb(t, 2 * semi - 1)

            call setDeltaB(deltaB + 1, t)

            call getDoubleMatrixElement(0, 1, deltaB, excitInfo%gen1, &
                                        excitInfo%gen2, bVal, excitInfo%order1, &
                                        tempWeight_0, tempWeight_1)

        case (2)
            if (t_approx_exchange .or. (t_approx_exchange_noninits .and. (.not. is_init_guga))) then
                if (getDeltaB(t) == -2) then
                    t = 0_n_int
                    probWeight = 0.0_dp
                    return
                end if
            end if

            ASSERT(getDeltaB(t) /= -2)

            ! do 2 -> 0
            clr_orb(t, 2 * semi)

            call setDeltaB(deltaB - 1, t)

            call getDoubleMatrixElement(0, 2, deltaB, excitInfo%gen1, &
                                        excitInfo%gen2, bVal, excitInfo%order1, &
                                        tempWeight_0, tempWeight_1)

        case (3)
            ! its a 3 -> have to check b if a 0 branch arrives
            select case (deltaB)
            case (-2)
                ! do 3 -> 2
                clr_orb(t, 2 * semi - 1)

                call getDoubleMatrixElement(2, 3, deltaB, excitInfo%gen1, &
                                            excitInfo%gen2, bVal, excitInfo%order1, &
                                            x1_element=tempWeight_1)

                tempWeight_0 = 0.0_dp

                call setDeltaB(-1, t)

            case (2)
                ! do 3 -> 1
                clr_orb(t, 2 * semi)

                call setDeltaB(+1, t)

                call getDoubleMatrixElement(1, 3, deltaB, excitInfo%gen1, &
                                            excitInfo%gen2, bVal, excitInfo%order1, &
                                            x1_element=tempWeight_1)

                tempWeight_0 = 0.0_dp

            case (0)
                ! deltaB = 0 branch arrives -> check b
                if (csf_i%B_int(semi) == 0) then
                    ! only 3 -> 1 possble
                    clr_orb(t, 2 * semi)

                    call setDeltaB(-1, t)

                    call getDoubleMatrixElement(1, 3, deltaB, excitInfo%gen1, &
                                                excitInfo%gen2, bVal, excitInfo%order1, &
                                                tempWeight_0, tempWeight_1)

                else
                    ! finally have to check the probs
                    minusWeight = weights%proc%minus(negSwitches(semi), &
                                                     bVal, weights%dat)
                    plusWeight = weights%proc%plus(posSwitches(semi), &
                                                   bVal, weights%dat)

                    if (near_zero(minusWeight + plusWeight)) then
                        probWeight = 0.0_dp
                        t = 0
                        return
                    end if

                    ! here i cant directly save it to probWeight ...
                    minusWeight = calcStartProb(minusWeight, plusWeight)

                    if (genrand_real2_dSFMT() < minusWeight) then
                        ! do -1 branch:
                        ! set 3 -> 1
                        clr_orb(t, 2 * semi)

                        call setDeltaB(-1, t)

                        call getDoubleMatrixElement(1, 3, deltaB, excitInfo%gen1, &
                                                    excitInfo%gen2, bVal, excitInfo%order1, &
                                                    tempWeight_0, tempWeight_1)

                        probWeight = probWeight * minusWeight

                    else
                        ! do +1 branch
                        ! set 3 -> 2
                        clr_orb(t, 2 * semi - 1)

                        call setDeltaB(1, t)

                        call getDoubleMatrixElement(2, 3, deltaB, excitInfo%gen1, &
                                                    excitInfo%gen2, bVal, excitInfo%order1, &
                                                    tempWeight_0, tempWeight_1)

                        probWeight = probWeight * (1 - minusWeight)
                    end if
                end if
            end select
        end select

        ! if its a mixed full-start raising into lowering -> check i a
        ! switch happened in the double overlap region
        if (excitInfo%typ == excit_type%fullstart_R_to_L) then
            ! this is indicated by a non-zero x0-matrix element
            if (.not. near_zero(extract_matrix_element(t, 1) * tempWeight_0)) then
                probWeight = 0.0_dp
                t = 0
                return
            end if
        end if

        ! do not combine the x0 and x1 elements at this point, because its
        ! easier to deal with the contributing integrals if we keep them
        ! seperate until the end!

        if (near_zero(extract_matrix_element(t, 1) * tempWeight_0) .and. &
            near_zero(extract_matrix_element(t, 2) * tempWeight_1)) then
            probWeight = 0.0_dp
            t = 0_n_int
            return
        end if

        call update_matrix_element(t, tempWeight_0, 1)
        call update_matrix_element(t, tempWeight_1, 2)

    end subroutine calcRaisingSemiStopStochastic

    subroutine calcLoweringSemiStopStochastic(ilut, csf_i, excitInfo, weights, negSwitches, &
                                              posSwitches, t, probWeight)
        integer(n_int), intent(in) :: ilut(0:nifguga)
        type(CSF_Info_t), intent(in) :: csf_i
        type(ExcitationInformation_t), intent(in) :: excitInfo
        type(WeightObj_t), intent(in) :: weights
        real(dp), intent(in) :: negSwitches(nSpatOrbs), posSwitches(nSpatOrbs)
        integer(n_int), intent(inout) :: t(0:nifguga)
        real(dp), intent(inout) :: probWeight
        character(*), parameter :: this_routine = "calcLoweringSemiStopStochastic"

        integer :: semi, deltaB
        real(dp) :: tempWeight_0, tempWeight_1, minusWeight, &
                    plusWeight, bVal

        ASSERT(.not. isThree(ilut, excitInfo%firstEnd))

        semi = excitInfo%firstEnd
        bVal = csf_i%B_real(semi)

        ! in the stochastic case i am sure that at there is no 3 or 0 at the
        ! full start... so definetly all deltaB branches can arrive here
        ! first deal with the forced ones
        deltaB = getDeltaB(t)

        select case (csf_i%stepvector(semi))
        case (1)
            if (t_approx_exchange .or. (t_approx_exchange_noninits .and. (.not. is_init_guga))) then
                if (getDeltaB(t) == 2) then
                    t = 0
                    return
                end if
            end if

            ASSERT(getDeltaB(t) /= 2)

            ! do the 1 -> 3 switch
            set_orb(t, 2 * semi)

            call setDeltaB(deltaB + 1, t)

            call getDoubleMatrixElement(3, 1, deltaB, excitInfo%gen1, &
                                        excitInfo%gen2, bVal, excitInfo%order1, &
                                        tempWeight_0, tempWeight_1)

        case (2)
            if (t_approx_exchange .or. (t_approx_exchange_noninits .and. (.not. is_init_guga))) then
                if (getDeltaB(t) == -2) then
                    t = 0
                    return
                end if
            end if

            ASSERT(getDeltaB(t) /= -2)

            ! do 2 -> 3
            set_orb(t, 2 * semi - 1)

            call setDeltaB(deltaB - 1, t)

            call getDoubleMatrixElement(3, 2, deltaB, excitInfo%gen1, &
                                        excitInfo%gen2, bVal, excitInfo%order1, &
                                        tempWeight_0, tempWeight_1)

        case (0)
            ! its a 0 -> have to check b if a 0 branch arrives
            select case (deltaB)
            case (-2)
                ! do 0 -> 2
                set_orb(t, 2 * semi)

                call getDoubleMatrixElement(2, 0, deltaB, excitInfo%gen1, &
                                            excitInfo%gen2, bVal, excitInfo%order1, &
                                            x1_element=tempWeight_1)

                tempWeight_0 = 0.0_dp

                call setDeltaB(-1, t)

            case (2)
                ! do 0 -> 1
                set_orb(t, 2 * semi - 1)

                call setDeltaB(+1, t)

                call getDoubleMatrixElement(1, 0, deltaB, excitInfo%gen1, &
                                            excitInfo%gen2, bVal, excitInfo%order1, &
                                            x1_element=tempWeight_1)

                tempWeight_0 = 0.0_dp

            case (0)
                ! deltaB=0 branch arrives -> check b
                if (csf_i%B_int(semi) == 0) then
                    ! only 0 -> 1 possble
                    set_orb(t, 2 * semi - 1)

                    call setDeltaB(-1, t)

                    call getDoubleMatrixElement(1, 0, deltaB, excitInfo%gen1, &
                                                excitInfo%gen2, bVal, excitInfo%order1, &
                                                tempWeight_0, tempWeight_1)

                else
                    ! finally have to check the probs
                    minusWeight = weights%proc%minus(negSwitches(semi), &
                                                     bVal, weights%dat)
                    plusWeight = weights%proc%plus(posSwitches(semi), &
                                                   bVal, weights%dat)

                    if (near_zero(minusWeight + plusWeight)) then
                        probWeight = 0.0_dp
                        t = 0
                        return
                    end if

                    ! here i cant directly save it to probWeight ...
                    minusWeight = calcStartProb(minusWeight, plusWeight)

                    if (genrand_real2_dSFMT() < minusWeight) then
                        ! do -1 branch:
                        ! set 0 -> 1
                        set_orb(t, 2 * semi - 1)

                        call setDeltaB(-1, t)

                        call getDoubleMatrixElement(1, 0, deltaB, excitInfo%gen1, &
                                                    excitInfo%gen2, bVal, excitInfo%order1, &
                                                    tempWeight_0, tempWeight_1)

                        probWeight = probWeight * minusWeight

                    else
                        ! do +1 branch
                        ! set 0 -> 2
                        set_orb(t, 2 * semi)

                        call setDeltaB(1, t)

                        call getDoubleMatrixElement(2, 0, deltaB, excitInfo%gen1, &
                                                    excitInfo%gen2, bVal, excitInfo%order1, &
                                                    tempWeight_0, tempWeight_1)

                        probWeight = probWeight * (1.0_dp - minusWeight)
                    end if
                end if
            end select
        end select

        ! for mixed fullstart check if no switch happened in the double
        ! overlap region, indicated by a non-zero-x0 matrix element
        if (excitInfo%typ == excit_type%fullStart_L_to_R) then
            if (.not. near_zero(extract_matrix_element(t, 1) * tempWeight_0)) then
                probWeight = 0.0_dp
                t = 0
                return
            end if
        end if

        if (near_zero(extract_matrix_element(t, 1) * tempWeight_0) .and. &
            near_zero(extract_matrix_element(t, 2) * tempWeight_1)) then
            probWeight = 0.0_dp
            t = 0
            return
        end if

        call update_matrix_element(t, tempWeight_0, 1)
        call update_matrix_element(t, tempWeight_1, 2)

    end subroutine calcLoweringSemiStopStochastic

    subroutine calcFullStartR2L_stochastic(ilut, csf_i, excitInfo, t, 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(inout) :: excitInfo
        integer(n_int), intent(out) :: t(0:nifguga)
        real(dp), intent(out) :: pgen
        real(dp), intent(in) :: posSwitches(nSpatOrbs), negSwitches(nSpatOrbs)
        type(WeightObj_t), intent(in), optional :: opt_weight
        character(*), parameter :: this_routine = "calcFullStartR2L_stochastic"

        integer :: i, st, en, se, gen, j, k, l, holeInd, elecInd
        type(WeightObj_t) :: weights
        real(dp) :: branch_pgen, temp_pgen, &
                    p_orig, rdm_mat, orb_pgen
        HElement_t(dp) :: integral
        procedure(calc_pgen_general), pointer :: calc_pgen_yix_start

        ASSERT(isProperCSF_ilut(ilut))
        ASSERT(.not. isZero(ilut, excitInfo%fullStart))
        ASSERT(.not. isThree(ilut, excitInfo%fullStart))

        ! create the fullStart
        st = excitInfo%fullStart
        en = excitInfo%fullEnd
        se = excitInfo%firstEnd
        gen = excitInfo%lastGen

        ! create correct weights:
        if (present(opt_weight)) then
            weights = opt_weight
        else
            weights = init_fullStartWeight(csf_i, se, en, negSwitches(se), &
                                           posSwitches(se), csf_i%B_real(se))
        end if

        if (t_approx_exchange .or. (t_approx_exchange_noninits .and. (.not. is_init_guga))) then
            ! todo the switch
            call forced_mixed_start(ilut, csf_i, excitInfo, t, branch_pgen)
        else
            call mixedFullStartStochastic(ilut, csf_i, excitInfo, weights, posSwitches, &
                                          negSwitches, t, branch_pgen)
        end if

        ! set pgen to 0 in case of early exit
        pgen = 0.0_dp

        ! check validity
        check_abort_excit(branch_pgen, t)

        ! then for the overlap region i need a double update routine, which
        ! somehow garuantees a switch happens at some point, to avoid
        ! single excitations

        do i = st + 1, se - 1
            call doubleUpdateStochastic(ilut, csf_i, i, excitInfo, &
                                        weights, negSwitches, posSwitches, t, branch_pgen)

            ! check validity
            if (near_zero(extract_matrix_element(t, 2)) .or. near_zero(branch_pgen)) then
                t = 0_n_int
                return
            end if
        end do

        ! then deal with specific semi-stop
        ! and update weights here
        ! i also could use the fact that the single weights are already
        ! initialized within the fullstart weights or?
        ! and then use smth like
        weights = weights%ptr

        call calcRaisingSemiStopStochastic(ilut, csf_i, excitInfo, weights, negSwitches, &
                                           posSwitches, t, branch_pgen)

        ! check validity
        check_abort_excit(branch_pgen, t)

        excitInfo%currentGen = excitInfo%lastGen

        do i = se + 1, en - 1
            call singleStochasticUpdate(ilut, csf_i, i, excitInfo, weights, posSwitches, &
                                        negSwitches, t, temp_pgen)
            branch_pgen = branch_pgen * temp_pgen

            if (t_trunc_guga_pgen .or. &
                (t_trunc_guga_pgen_noninits .and. .not. is_init_guga)) then
                if (branch_pgen < trunc_guga_pgen) then
                    t = 0_n_int
                    return
                end if
            end if

            ! check validity
            check_abort_excit(branch_pgen, t)

        end do

        call singleStochasticEnd(csf_i, excitInfo, t)

        ! if we do RDMs also store the x0 and x1 coupling coeffs
        ! and I need to do it before the routines below since excitInfo
        ! gets changed there
        if (tFillingStochRDMOnFly) then
            ! i need to unbias against the total pgen later on in the
            ! RDM sampling otherwise the rdm-bias factor is not correct!
            ! encode the necessary information in the rdm-matele!
            i = excitInfo%i
            j = excitInfo%j
            k = excitInfo%k
            l = excitInfo%l
            elecInd = se
            holeInd = en
            rdm_mat = extract_matrix_element(t, 2)
            call calc_orbital_pgen_contrib_start(&
                    csf_i, [2 * st, 2 * elecInd], holeInd, orb_pgen)
            p_orig = orb_pgen * branch_pgen / real(ElecPairs, dp)
            if (csf_i%stepvector(elecInd) == 3) p_orig = p_orig * 2.0_dp
        end if

        call encode_matrix_element(t, extract_matrix_element(t, 1) + &
                                   extract_matrix_element(t, 2), 1)
        call encode_matrix_element(t, 0.0_dp, 2)

        ! update: since i can formulate everything in terms of the already
        ! calculated matrix element i can abort here if it is zero
        if (near_zero(extract_matrix_element(t, 1))) then
            t = 0_n_int
            return
        end if

        global_excitInfo = excitInfo

        if (t_approx_exchange .or. (t_approx_exchange_noninits .and. (.not. is_init_guga))) then
            call calc_mixed_start_contr_approx(t, csf_i, excitInfo, integral)
            pgen = branch_pgen

        else
            call calc_mixed_start_contr_integral(ilut, csf_i, t, excitInfo, &
                integral)
            call calc_mixed_start_contr_pgen(ilut, csf_i, t, excitInfo, &
                branch_pgen, pgen)

!             call calc_mixed_start_contr_sym(ilut, csf_i, t, excitInfo, branch_pgen, pgen, &
!                                             integral)
        end if

        if (tFillingStochRDMOnFly) then
            if (.not. near_zero(p_orig)) then
                call encode_stochastic_rdm_info(GugaBits, t, rdm_ind= &
                                    contract_2_rdm_ind(i, j, k, l, excit_lvl=2, &
                                       excit_typ=excitInfo%typ), x0=0.0_dp, &
                                                x1=rdm_mat * pgen / p_orig)
            end if
        end if

        ! and finally update the matrix element with all contributions
        call update_matrix_element(t, (get_umat_el(st, en, se, st) + &
                                       get_umat_el(en, st, st, se)) / 2.0_dp + integral, 1)

    end subroutine calcFullStartR2L_stochastic

    subroutine calc_mixed_start_contr_approx(t, csf_i, excitInfo, integral)
        integer(n_int), intent(in) :: t(0:nifguga)
        type(CSF_Info_t), intent(in) :: csf_i
        type(ExcitationInformation_t), intent(in) :: excitInfo
        HElement_t(dp), intent(out) :: integral
        character(*), parameter :: this_routine = "calc_mixed_start_contr_approx"

        integer :: st, se, en, elecInd, holeInd, sw, step, i
        real(dp) :: bot_cont, mat_ele, stay_mat, start_mat
        logical :: below_flag

        st = excitInfo%fullStart
        se = excitInfo%firstEnd
        en = excitInfo%fullEnd
        ! depending on the type of excitaiton, calculation of orbital pgens
        ! change
        if (excitInfo%typ == excit_type%fullStart_L_to_R) then
            elecInd = en
            holeInd = se
        else if (excitInfo%typ == excit_type%fullstart_R_to_L) then
            elecInd = se
            holeInd = en
        else
            call stop_all(this_routine, "should not be here!")
        end if

        ! I am sure first switch is at full-start
        sw = st

        ! what can i precalculate beforehand?
        step = csf_i%stepvector(st)

        integral = h_cast(0.0_dp)

        if (step == 1) then

            ASSERT(isTwo(t, st))

            bot_cont = -sqrt(2.0_dp / ((csf_i%B_real(st) - 1.0_dp) * &
                                       (csf_i%B_real(st) + 1.0_dp)))

        else

            ASSERT(isOne(t, st))

            bot_cont = -sqrt(2.0_dp / ((csf_i%B_real(st) + 1.0_dp) * &
                                       (csf_i%B_real(st) + 3.0_dp)))

        end if

        if (.not. near_zero(bot_cont)) then

            mat_ele = 1.0_dp
            below_flag = .false.

            do i = st - 1, 1, -1
                if (csf_i%Occ_int(i) /= 1) cycle

                ! then check if thats the last stepvalue to consider
                if (csf_i%stepvector(i) == 1 .and. csf_i%B_int(i) == 1) then
                    below_flag = .true.
                end if

                ! then deal with the matrix element and branching probabilities
                step = csf_i%stepvector(i)

                ! get both start and staying matrix elements -> and update
                ! matrix element contributions on the fly to avoid second loop!
                call getDoubleMatrixElement(step, step, -1, gen_type%R, gen_type%L, csf_i%B_real(i), &
                                            1.0_dp, x1_element=start_mat)

                call getDoubleMatrixElement(step, step, 0, gen_type%R, gen_type%L, csf_i%B_real(i), &
                                            1.0_dp, x1_element=stay_mat)

                if (near_zero(stay_mat)) below_flag = .true.
                ! "normally" matrix element shouldnt be 0 anymore... still check

                if (.not. near_zero(start_mat)) then
                    integral = integral + start_mat * mat_ele * (get_umat_el(i, holeInd, elecInd, i) &
                                                                 + get_umat_el(holeInd, i, i, elecInd)) / 2.0_dp
                end if

                ! also update matrix element on the fly
                mat_ele = stay_mat * mat_ele

            end do
        end if

        ! maybe I need to deal with the pgens here, if they are not
        ! correctly considered outside..

    end subroutine calc_mixed_start_contr_approx

    subroutine calcFullStartL2R_stochastic(ilut, csf_i, excitInfo, t, 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(inout) :: excitInfo
        integer(n_int), intent(out) :: t(0:nifguga)
        real(dp), intent(out) :: pgen
        real(dp), intent(in) :: posSwitches(nSpatOrbs), negSwitches(nSpatOrbs)
        type(WeightObj_t), intent(in), optional :: opt_weight
        character(*), parameter :: this_routine = "calcFullStartL2R_stochastic"

        integer :: i, st, en, se, gen, j, k, l, elecInd, holeInd
        type(WeightObj_t) :: weights
        real(dp) :: branch_pgen, temp_pgen, p_orig, orb_pgen, rdm_mat
        HElement_t(dp) :: integral
        procedure(calc_pgen_general), pointer :: calc_pgen_yix_start

        ! the integral contributions to this mixed excitaiton full starts are
        ! more involved then i have imagined. since all the orbitals from below
        ! the first stepvector change(which has to happen at some point of the
        ! excitation, or else it is single-like,(which is also still unclear
        ! how i should implement that)) contribute to the same excitation.
        ! similar to the full-start full-stop mixed generator case.
        ! to ensure a stepvector change i could write some kind of restricted
        ! double excitiation update function, which accumulates switch
        ! probability with decreasing switch possibilities...
        ! and then during the excitation creation i have to save the first
        ! stepvector change and depending on that calculate the remaining
        ! integral contribution. todo
        ! and same with the fullstartr2l, and fullstoprl2/l2r
        ! for both full starts and full stop routines, the contributions are
        ! the same atleast.
        ASSERT(isProperCSF_ilut(ilut))
        ASSERT(.not. isZero(ilut, excitInfo%fullStart))
        ASSERT(.not. isThree(ilut, excitInfo%fullStart))

        ! create the fullStart
        st = excitInfo%fullStart
        en = excitInfo%fullEnd
        se = excitInfo%firstEnd
        gen = excitInfo%lastGen

        ! create correct weights:
        if (present(opt_weight)) then
            weights = opt_weight
        else
            weights = init_fullStartWeight(csf_i, se, en, negSwitches(se), &
                                           posSwitches(se), csf_i%B_real(se))
        end if

        ! in the case of the approximate exchange excitations I need to
        ! force a switch at the beginning
        if (t_approx_exchange .or. (t_approx_exchange_noninits .and. (.not. is_init_guga))) then
            ! do the switch
            call forced_mixed_start(ilut, csf_i, excitInfo, t, branch_pgen)

        else
            call mixedFullStartStochastic(ilut, csf_i, excitInfo, weights, posSwitches, &
                                          negSwitches, t, branch_pgen)
        end if

        ! in case of early exit pgen should be set to 0
        pgen = 0.0_dp

        check_abort_excit(branch_pgen, t)

        ! in the mixed fullstart case there has to be a switch at some point
        ! in the double overlap region, or else it would be a single-like
        ! excitation -> so check if x1 matrix element is 0

        ! then for the overlap region i need a double update routine, which
        ! somehow garuantees a switch happens at some point, to avoid
        ! single excitations
        do i = st + 1, se - 1
            call doubleUpdateStochastic(ilut, csf_i, i, excitInfo, &
                                        weights, negSwitches, posSwitches, t, branch_pgen)

            ! to keep it general, i cant only check weights in doubleUpdate
            if (near_zero(extract_matrix_element(t, 2)) .or. near_zero(branch_pgen)) then
                t = 0_n_int
                return
            end if
        end do

        ! then deal with specific semi-stop
        ! and update weights here
        weights = weights%ptr

        call calcLoweringSemiStopStochastic(ilut, csf_i, excitInfo, weights, negSwitches, &
                                            posSwitches, t, branch_pgen)

        ! check validity
        check_abort_excit(branch_pgen, t)

        excitInfo%currentGen = excitInfo%lastGen

        do i = se + 1, en - 1
            call singleStochasticUpdate(ilut, csf_i, i, excitInfo, weights, posSwitches, &
                                        negSwitches, t, temp_pgen)
            branch_pgen = branch_pgen * temp_pgen

            if (t_trunc_guga_pgen .or. &
                (t_trunc_guga_pgen_noninits .and. .not. is_init_guga)) then
                if (branch_pgen < trunc_guga_pgen) then
                    t = 0_n_int
                    return
                end if
            end if

            ! check validity
            check_abort_excit(branch_pgen, t)
        end do

        call singleStochasticEnd(csf_i, excitInfo, t)

        ! if we do RDMs also store the x0 and x1 coupling coeffs
        ! and I need to do it before the routines below since excitInfo
        ! gets changed there
        if (tFillingStochRDMOnFly) then
            ! i need to unbias against the total pgen later on in the
            ! RDM sampling otherwise the rdm-bias factor is not correct!
            ! encode the necessary information in the rdm-matele!
            i = excitInfo%i
            j = excitInfo%j
            k = excitInfo%k
            l = excitInfo%l
            elecInd = en
            holeInd = se
            rdm_mat = extract_matrix_element(t, 2)
            call calc_orbital_pgen_contrib_start(&
                    csf_i, [2 * st, 2 * elecInd], holeInd, orb_pgen)
            p_orig = orb_pgen * branch_pgen / real(ElecPairs, dp)
            if (csf_i%stepvector(elecInd) == 3) p_orig = p_orig * 2.0_dp
        end if

        ! put everything in first entry
        call encode_matrix_element(t, extract_matrix_element(t, 1) + &
                                   extract_matrix_element(t, 2), 1)
        call encode_matrix_element(t, 0.0_dp, 2)

        ! now have to check if there has be a change in the double overlap
        ! region, or else its just a single-like excitation:
        ! todo write a routine for that
        ! switchFlag = checkForSwitch(ilut, t, st, se-1)

        ! if a switch happened i have to calculated the additional contributing
        ! matrix elements from all indices below the fullstart
        ! can i formulate that in terms of the already calulated matrix
        ! element, as in the diagonal matrix element calculation?
        ! probably... but todo!

        ! update: since i can formulate everything in terms of the already
        ! calculated matrix element i can abort here if it is zero
        if (near_zero(extract_matrix_element(t, 1))) then
            t = 0_n_int
            return
        end if

        global_excitInfo = excitInfo

        if (t_approx_exchange .or. (t_approx_exchange_noninits .and. (.not. is_init_guga))) then
            call calc_mixed_start_contr_approx(t, csf_i, excitInfo, integral)
            pgen = branch_pgen
        else
            call calc_mixed_start_contr_integral(ilut, csf_i, t, excitInfo, &
                integral)
            call calc_mixed_start_contr_pgen(ilut, csf_i, t, excitInfo, &
                branch_pgen, pgen)

!             call calc_mixed_start_contr_sym(ilut, csf_i, t, excitInfo, branch_pgen, pgen, integral)
        end if

        if (tFillingStochRDMOnFly) then
            if (.not. near_zero(p_orig)) then
                call encode_stochastic_rdm_info(GugaBits, t, rdm_ind= &
                            contract_2_rdm_ind(i, j, k, l, excit_lvl=2, &
                                       excit_typ=excitInfo%typ), x0=0.0_dp, &
                                                x1=rdm_mat * pgen / p_orig)
            end if
        end if

        ! and finally update the matrix element with all contributions
        call update_matrix_element(t, (get_umat_el(st, se, en, st) + &
                                       get_umat_el(se, st, st, en)) / 2.0_dp + integral, 1)

    end subroutine calcFullStartL2R_stochastic

    subroutine calc_mixed_x2x_ueg(ilut, csf_i, t, excitInfo, branch_pgen, pgen, &
                                  integral, rdm_ind, rdm_mat)
        integer(n_int), intent(in) :: ilut(0:nifguga), t(0:nifguga)
        type(CSF_Info_t), intent(in) :: csf_i
        type(ExcitationInformation_t), intent(inout) :: excitInfo
        real(dp), intent(inout) :: branch_pgen
        real(dp), intent(out) :: pgen
        HElement_t(dp), intent(out) :: integral
        integer(int_rdm), intent(out), allocatable, optional :: rdm_ind(:)
        real(dp), intent(out), allocatable, optional :: rdm_mat(:)
        character(*), parameter :: this_routine = "calc_mixed_x2x_ueg"

        pgen = 0.0_dp
        integral = 0.0_dp
        unused_var(ilut); unused_var(t); unused_var(excitInfo); unused_var(branch_pgen);
        unused_var(csf_i)
        if (present(rdm_ind)) then
            allocate(rdm_ind(0), source=0_int_rdm)
        end if
        if (present(rdm_mat)) then
            allocate(rdm_mat(0), source=0.0_dp)
        end if
        call stop_all(this_routine, &
                      "in Hubbard/UEG calculations with full k-point symmetry, this excitation shouldnt be reached!")
    end subroutine calc_mixed_x2x_ueg

    subroutine calc_orbital_pgen_contrib_end_def(csf_i, occ_orbs, orb_a, orb_pgen)
        ! write a combined function for both r2l and l2r since its only
        ! one difference -> only one if condition to adjust for both!
        type(CSF_Info_t), intent(in) :: csf_i
        integer, intent(in) :: occ_orbs(2), orb_a
        real(dp), intent(out) :: orb_pgen

        integer :: i, j, orb
        real(dp) :: cum_sum, cpt_a, cpt_b, cpt_ba, cpt_ab, ba_sum, ab_sum

        ! electron indices
        i = gtID(occ_orbs(1))
        j = gtID(occ_orbs(2))

        cum_sum = 0.0_dp
        if (tGen_guga_weighted) then
            do orb = 1, orb_a - 1
                ! calc. the p(a)
                if (csf_i%stepvector(orb) /= 3) then
                    cum_sum = cum_sum + get_guga_integral_contrib(occ_orbs, orb, -1)

                end if

            end do
        end if

        ! deal with orb (a) in specific way:
        cpt_a = get_guga_integral_contrib(occ_orbs, orb_a, -1)

        cum_sum = cum_sum + cpt_a

        ! also get p(b|a)
        ! depending if its a r2l or l2r full-stop:
        if (i < orb_a) then
            ! its a L2R -> so no restrictions
            call pgen_select_orb_guga_mol(csf_i, occ_orbs, orb_a, j, cpt_ba, ba_sum)
        else
            ! its a R2L so orbital i is off-limits
            call pgen_select_orb_guga_mol(csf_i, occ_orbs, orb_a, j, cpt_ba, ba_sum, i)
        end if

        ! change to the fullstart into fullstop: loop until orbital j for the
        ! fullstop implementattion
        if (tGen_guga_weighted) then
            do orb = orb_a + 1, j - 1
                if (csf_i%stepvector(orb) /= 3) then
                    cum_sum = cum_sum + get_guga_integral_contrib(occ_orbs, orb, -1)
                end if
            end do
        end if

        ! deal with j also speciallly
        cpt_b = get_guga_integral_contrib(occ_orbs, j, -1)

        cum_sum = cum_sum + cpt_b

        ! and get p(a|b)
        ! only orbitals below j are allowed!
        call pgen_select_orb_guga_mol(csf_i, occ_orbs, j, orb_a, cpt_ab, ab_sum, -j, .true.)

        ! and deal with rest:
        if (tGen_guga_weighted) then
            do orb = j + 1, nSpatOrbs
                if (csf_i%stepvector(orb) /= 3) then
                    cum_sum = cum_sum + get_guga_integral_contrib(occ_orbs, orb, -1)
                end if
            end do
        end if

        if (.not. tGen_guga_weighted) then
            cum_sum = csf_i%cum_list(nSpatOrbs)
        end if

        if (near_zero(cum_sum) .or. near_zero(ab_sum) .or. near_zero(ba_sum)) then
            orb_pgen = 0.0_dp
        else
            ! and get hopefully correct final values:
            cpt_a = cpt_a / cum_sum * cpt_ba / ba_sum
            cpt_b = cpt_b / cum_sum * cpt_ab / ab_sum

            ! and add them up to the final orbital pgen
            orb_pgen = cpt_a + cpt_b
        end if

    end subroutine calc_orbital_pgen_contrib_end_def

    subroutine calc_orbital_pgen_contrib_start_def(csf_i, occ_orbs, orb_a, orb_pgen)
        ! write a combined function for both r2l and l2r since its only
        ! one difference -> only one if condition to adjust for both!
        type(CSF_Info_t), intent(in) :: csf_i
        integer, intent(in) :: occ_orbs(2), orb_a
        real(dp), intent(out) :: orb_pgen

        integer :: i, j, orb
        real(dp) :: cum_sum, cpt_a, cpt_b, cpt_ba, cpt_ab, ba_sum, ab_sum

        ! electron indices
        i = gtID(occ_orbs(1))
        j = gtID(occ_orbs(2))

        ! damn... i need the probability of the elec-pair picking too or?

        cum_sum = 0.0_dp
        if (tGen_guga_weighted) then
            do orb = 1, i - 1
                ! calc. the p(a)
                if (csf_i%stepvector(orb) /= 3) then
                    cum_sum = cum_sum + get_guga_integral_contrib(occ_orbs, orb, -1)
                end if

            end do
        end if

        ! deal with orb (i) in specific way:
        cpt_a = get_guga_integral_contrib(occ_orbs, i, -1)

        cum_sum = cum_sum + cpt_a

        ! also get p(b|a)
        ! did i mix up i and j here below.. what do i assume picked already
        ! here? if its really p(b|a) i should switch up i and j..
        ! hm.. the inputted j is the holeInd i which is fixed.. i gets looped
        ! over on the outside and is assumed picked first or 2nd here??
        ! taking i and j here is wrong! i is the open orbital, but j
        ! is the already picked electron! it has to be orb_a here or?
        call pgen_select_orb_guga_mol(csf_i, occ_orbs, i, orb_a, cpt_ba, ba_sum, i, .true.)

        ! change to the fullstart into fullstop: loop until orbital a
        if (tGen_guga_weighted) then
            do orb = i + 1, orb_a - 1
                if (csf_i%stepvector(orb) /= 3) then
                    cum_sum = cum_sum + get_guga_integral_contrib(occ_orbs, orb, -1)
                end if
            end do
        end if

        ! deal with a also speciallly
        cpt_b = get_guga_integral_contrib(occ_orbs, orb_a, -1)

        cum_sum = cum_sum + cpt_b

        ! and get p(a|b)
        ! here the only difference between r2l and l2r fullstarts come into
        ! play!
        if (orb_a > j) then
            ! then orb_j is off-limits
            call pgen_select_orb_guga_mol(csf_i, occ_orbs, orb_a, i, cpt_ab, ab_sum, j)
        else
            ! in this case there is no restriction guga-wise..
            call pgen_select_orb_guga_mol(csf_i, occ_orbs, orb_a, i, cpt_ab, ab_sum)
        end if

        ! and deal with rest:
        if (tGen_guga_weighted) then
            do orb = orb_a + 1, nSpatOrbs
                if (csf_i%stepvector(orb) /= 3) then
                    cum_sum = cum_sum + get_guga_integral_contrib(occ_orbs, orb, -1)
                end if
            end do
        end if

        if (.not. tGen_guga_weighted) then
            cum_sum = csf_i%cum_list(nSpatOrbs)
        end if

        if (near_zero(cum_sum) .or. near_zero(ab_sum) .or. near_zero(ba_sum)) then
            orb_pgen = 0.0_dp
        else
            ! and get hopefully correct final values:
            cpt_a = cpt_a / cum_sum * cpt_ba / ba_sum
            cpt_b = cpt_b / cum_sum * cpt_ab / ab_sum

            ! and add them up to the final orbital pgen
            orb_pgen = cpt_a + cpt_b
        end if

    end subroutine calc_orbital_pgen_contrib_start_def

    subroutine forced_mixed_start(ilut, csf_i, excitInfo, t, probWeight)
        ! NOTE: mixed full-start matrix elements are stores in the same row
        ! as delta B = -1 ones -> so access getDoubleMatrixElement with
        ! db = -1 below!
        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) :: probWeight
        character(*), parameter :: this_routine = "forced_mixed_start"

        real(dp) :: tempWeight, bVal, tempWeight_1
        integer :: st

        ASSERT(isProperCSF_ilut(ilut))
        ASSERT(.not. isZero(ilut, excitInfo%fullStart))
        ASSERT(.not. isThree(ilut, excitInfo%fullStart))

        ! cant be 0 or zero matrix element, but cant be 3 either or only
        ! deltaB=0 branch would have non-zero matrix element and that would
        ! to a single-excitation-like DE

        st = excitInfo%fullStart
        bVal = csf_i%B_real(st)

        t = ilut

        select case (csf_i%stepvector(st))
        case (1)
            if (csf_i%B_int(st) < 2) then
                ! actually not possible in this case, as only 0 branch valid..
                probWeight = 0.0_dp
                t = 0_n_int
                return
            end if

            ! otherwise switch 1 -> 2
            clr_one(t, st)
            set_two(t, st)

            call getDoubleMatrixElement(2, 1, -1, gen_type%L, gen_type%R, &
                                        bVal, 1.0_dp, x1_element=tempWeight_1)

            ! x0 matrix element:
            tempWeight = 0.0_dp

            call setDeltaB(2, t)

        case (2)
            ! choose -2 branch 2 -> 1
            clr_two(t, st)
            set_one(t, st)

            call getDoubleMatrixElement(1, 2, -1, gen_type%L, gen_type%R, &
                                        bVal, 1.0_dp, x1_element=tempWeight_1)

            tempWeight = 0.0_dp

            call setDeltaB(-2, t)
#ifdef DEBUG_
        case default
            call stop_all(this_routine, "wrong stepvalue!")
#endif
        end select

        call encode_matrix_element(t, tempWeight_1, 2)
        call encode_matrix_element(t, tempWeight, 1)

        if (near_zero(abs(tempWeight) + abs(tempWeight_1))) then
            probWeight = 0.0_dp
            t = 0_n_int
            return
        end if

        ! and since there is no choice: branch_pgen is 1
        probWeight = 1.0_dp

    end subroutine forced_mixed_start

    subroutine mixedFullStartStochastic(ilut, csf_i, excitInfo, weights, posSwitches, &
                                        negSwitches, t, probWeight)
        ! NOTE: mixed full-start matrix elements are stores in the same row
        ! as delta B = -1 ones -> so access getDoubleMatrixElement with
        ! db = -1 below!
        integer(n_int), intent(in) :: ilut(0:nifguga)
        type(CSF_Info_t), intent(in) :: csf_i
        type(ExcitationInformation_t), intent(in) :: excitInfo
        type(WeightObj_t), intent(in) :: weights
        real(dp), intent(in) :: posSwitches(nSpatOrbs), negSwitches(nSpatOrbs)
        integer(n_int), intent(out) :: t(0:nifguga)
        real(dp), intent(out) :: probWeight
        character(*), parameter :: this_routine = "mixedFullStartStochastic"

        real(dp) :: minusWeight, plusWeight, zeroWeight, tempWeight, bVal, &
                    tempWeight_1
        integer :: st

        ASSERT(isProperCSF_ilut(ilut))
        ASSERT(.not. isZero(ilut, excitInfo%fullStart))
        ASSERT(.not. isThree(ilut, excitInfo%fullStart))

        ! cant be 0 or zero matrix element, but cant be 3 either or only
        ! deltaB=0 branch would have non-zero matrix element and that would
        ! to a single-excitation-like DE

        st = excitInfo%fullStart
        bVal = csf_i%B_real(st)

        t = ilut

        ! another note on the matrix elements... for a mixed full-start
        ! i know that the branch has to switch to +-2 at some point to
        ! yield a non-single excitation -> so i could just ignore the x0
        ! matrix element, as it definetly gets 0 at some point...
        ! but i also then have to write a specific double update function,
        ! which also deals with this kind of "restricted" or better enforced
        ! double excitation regime!

        ! no do not do that, as i have to check probweights inbetween if
        ! excitation yields non-zero matrix element
        select case (csf_i%stepvector(st))
        case (1)
            if (csf_i%B_int(st) < 2) then
                ! only 0-branch start possible, which always has weight > 0
                ! no change in stepvector, so just calc matrix element
                call getDoubleMatrixElement(1, 1, -1, gen_type%L, gen_type%R, &
                                            bVal, 1.0_dp, tempWeight, tempWeight_1)

                call setDeltaB(0, t)

                probWeight = 1.0_dp
            else
                ! both branches possible
                plusWeight = weights%proc%plus(posSwitches(st), &
                                               bVal, weights%dat)
                zeroWeight = weights%proc%zero(negSwitches(st), &
                                               posSwitches(st), bVal, weights%dat)

                probWeight = calcStartProb(zeroWeight, plusWeight)

                if (genrand_real2_dSFMT() < probWeight) then
                    ! choose 0 branch

                    call getDoubleMatrixElement(1, 1, -1, gen_type%L, gen_type%R, &
                                                bVal, 1.0_dp, tempWeight, tempWeight_1)

                    call setDeltaB(0, t)

                else
                    ! +2 branch:
                    ! switch 1 -> 2
                    clr_orb(t, 2 * st - 1)
                    set_orb(t, 2 * st)

                    call getDoubleMatrixElement(2, 1, -1, gen_type%L, gen_type%R, &
                                                bVal, 1.0_dp, x1_element=tempWeight_1)

                    tempWeight = 0.0_dp

                    call setDeltaB(2, t)

                    probWeight = 1.0_dp - probWeight
                end if

            end if

        case (2)
            ! here always both branches are possible
            minusWeight = weights%proc%minus(negSwitches(st), &
                                             bVal, weights%dat)
            zeroWeight = weights%proc%zero(negSwitches(st), &
                                           posSwitches(st), bVal, weights%dat)

            probWeight = calcStartProb(zeroWeight, minusWeight)

            if (genrand_real2_dSFMT() < probWeight) then
                ! choose 0 branch
                call getDoubleMatrixElement(2, 2, -1, gen_type%L, gen_type%R, &
                                            bVal, 1.0_dp, tempWeight, tempWeight_1)

                call setDeltaB(0, t)

            else
                ! choose -2 branch:
                ! change 2 -> 1
                set_orb(t, 2 * st - 1)
                clr_orb(t, 2 * st)

                call getDoubleMatrixElement(1, 2, -1, gen_type%L, gen_type%R, &
                                            bVal, 1.0_dp, x1_element=tempWeight_1)

                tempWeight = 0.0_dp

                call setDeltaB(-2, t)

                probWeight = 1.0_dp - probWeight

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

        call encode_matrix_element(t, tempWeight_1, 2)
        call encode_matrix_element(t, tempWeight, 1)

        ! since i only need this routine in excitations, where there has to be
        ! a switch in the double overlap part i could check if it is 0 and then
        ! set probWeight to 0 and return
        if (near_zero(abs(tempWeight) + abs(tempWeight_1))) then
            probWeight = 0.0_dp
            t = 0
            return
        end if

    end subroutine mixedFullStartStochastic

    subroutine calcSingleOverlapMixedStochastic(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(inout) :: excitInfo
        integer(n_int),