mixedFullStartStochastic Subroutine

    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