create_all_rdm_contribs Subroutine

    subroutine create_all_rdm_contribs(rdm_ind, x0, x1, rdm_ind_out, rdm_mat)
        ! I also need a function which creates me the remaining
        ! index combination of not directly sampled excitation in the
        ! GUGA excitation generation
        integer(int_rdm), intent(in) :: rdm_ind
        real(dp) :: x0, x1
        integer(int_rdm), intent(out), allocatable :: rdm_ind_out(:)
        real(dp), intent(out), allocatable :: rdm_mat(:)
        debug_function_name("create_all_rdm_contribs")
        integer :: ex_lvl, ex_typ, i, j, k, l

        ex_lvl = extract_excit_lvl_rdm(rdm_ind)
        ex_typ = extract_excit_type_rdm(rdm_ind)

        ASSERT(ex_lvl == 1 .or. ex_lvl == 2)
        ASSERT(ex_typ /= excit_type%invalid)

        if (ex_lvl == 1) then
            ! for singles we 'only' have the original index
            allocate(rdm_ind_out(1), source=rdm_ind)
            ! and the element is only the x0 contrib
            ASSERT(near_zero(x1))
            allocate(rdm_mat(1), source=x0)

        else if (ex_lvl == 2) then

            select case (ex_typ)

                ! we only take stuff from stochastic double excitations here..
                ! so this excludes some excitation types.. assert this below!
            case (excit_type%single_overlap_R_to_L, &
                  excit_type%single_overlap_L_to_R, &
                  excit_type%fullstop_lowering, &
                  excit_type%fullstop_raising, &
                  excit_type%fullstart_raising, &
                  excit_type%fullstart_lowering, &
                  excit_type%fullstart_stop_alike)

                ! here we also only have the ij <-> kl conjugation which
                ! can (and should for now) be handled outside of this function!
                allocate(rdm_ind_out(1), source=rdm_ind)
                ! also x1 must be 0 in these cases
                ASSERT(near_zero(x1))
                allocate(rdm_mat(1), source=x0)

                ! group everything together where x0 is 0:
                ! the conjugated elements of those are actually dealt with
                ! in the single excitations (here they would actually be 0
                ! since we enforce a spin-flip in the overlap range!
            case (excit_type%fullstop_L_to_R, &
                  excit_type%fullstop_R_to_L, &
                  excit_type%fullstart_L_to_R, &
                  excit_type%fullstart_R_to_L, &
                  excit_type%fullstart_stop_mixed)

                allocate(rdm_ind_out(1), source=rdm_ind)
                ASSERT(near_zero(x0))
                allocate(rdm_mat(1), source=x1)

            case (excit_type%double_lowering, &
                  excit_type%double_raising)

                ! here we need to do something
                call extract_2_rdm_ind(rdm_ind, i, j, k, l)
                ! a switch of j <-> l will induce a sign change between x0 and x1!
                allocate(rdm_ind_out(2), source=0_int_rdm)
                allocate(rdm_mat(2), source=0.0_dp)

                rdm_ind_out(1) = rdm_ind
                rdm_mat(1) = x0 + generator_sign(i, j, k, l) * x1

                rdm_ind_out(2) = contract_2_rdm_ind(i, l, k, j, &
                                                    excit_lvl=2, excit_typ=ex_typ)

                rdm_mat(2) = x0 + generator_sign(i, l, k, j) * x1

            case (excit_type%double_L_to_R_to_L, &
                  excit_type%double_R_to_L_to_R, &
                  excit_type%double_L_to_R, &
                  excit_type%double_R_to_L)

                ! these cases contain a non-overlap excitaion if no spin-flip
                ! in the overlap range happened (i hope for now that a
                ! zero x0 element indicates such a spin-flip.. although
                ! it should be fine, since when there is a spin-flip it
                ! is 0 definitely and even if there was, x0 is still 0 and this
                ! is all what is encoded for this flipped version!

                allocate(rdm_ind_out(2), source=0_int_rdm)
                allocate(rdm_mat(2), source=0.0_dp)

                ! the first is the original
                rdm_ind_out(1) = rdm_ind
                rdm_mat(1) = x0 + x1

                call extract_2_rdm_ind(rdm_ind, i, j, k, l)
                ! the second one does hopefully not need a type encoded
                ! in the rdm! since this is not available here!!
                rdm_ind_out(2) = contract_2_rdm_ind(i, l, k, j, excit_lvl=2)
                rdm_mat(2) = -2.0_dp * x0

            case default
                ! should not be here for now!
                ASSERT(.false.)

            end select
        else
            ! but should actually not be here now!
            ASSERT(.false.)
            ! other wise no contributions
            allocate(rdm_ind_out(0))
            allocate(rdm_mat(0))
        end if

    end subroutine create_all_rdm_contribs