fill_spawn_rdm_diag_guga Subroutine

    subroutine fill_spawn_rdm_diag_guga(spawn, nI, full_sign)
        ! i have to write a routine, which correctly takes into
        ! account all the diagonal contributions from double excitations
        type(rdm_spawn_t), intent(inout) :: spawn
        integer, intent(in) :: nI(nel)
        real(dp), intent(in) :: full_sign(spawn%rdm_send%sign_length)

        integer :: i, s_orbs(nel), s, inc_i, j, inc_j, p

        real(dp) :: occ_i, occ_j, x0, x1

        ! i have to figure out what exactly contributes to here..
        ! according to the paper about the two-RDM labels the
        ! diagonal elements of the 2-RDM are given by
        ! D_ij,kl = <psi| e_{ki,lj} | psi>
        ! which for the diagonal contribution kl = ij yields
        ! D_{ij,ij} = <psi| e_{ii,jj} |psi>
        ! which is definetly purely diagonal!

        ! i could try to use the add_to_rdm_spawn_t routine in
        ! rdm_data_utils with the according modified matrix elements.
        i = 1
        s_orbs = gtID(nI)

        ! TODO: I could also try to add the diagonal exchange
        ! contributions directly here! so this routine would be
        ! applicable in the stochastic spawning already
        do while (i <= nel)

            s = s_orbs(i)

            if (isDouble(nI, i)) then
                occ_i = 2.0_dp
                inc_i = 2

                ! i think here i have to put the full diagonal
                ! D_{ii,ii} entry, which counts double occupancies
                ! and i will try for now to reuse the following routine:
                ! put with the spatial orbital s
                ! but maybe I have to also add the 'switched' indices
                ! contribution.. which would mean a factor of 2..
                call add_to_rdm_spawn_t(spawn, s, s, s, s, occ_i * full_sign, .true.)

            else
                occ_i = 1.0_dp
                inc_i = 1
            end if

            j = i + inc_i

            do while (j <= nel)

                p = s_orbs(j)

                if (isDouble(nI, j)) then
                    occ_j = 2.0_dp
                    inc_j = 2
                else
                    occ_j = 1.0_dp
                    inc_j = 1
                end if

                ! i could also just multiply by 2 here, since this will
                ! get strored in the same D_{ij,ij} RDM element!
                ! but for now I decided to fill in all disting 2-RDM
                ! elements
                ! since here we do not have spawning this should also
                ! get stored symmetrically! as otherwise this will
                ! junk up the hermiticity error calculation!
                call add_to_rdm_spawn_t(spawn, s, s, p, p, &
                                        occ_i * occ_j * full_sign, .true.)
                call add_to_rdm_spawn_t(spawn, p, p, s, s, &
                                        occ_i * occ_j * full_sign, .true.)

                ! i can also add the fully diagonal exchange contributions
                ! here. this is also necessary to do, if I want to use
                ! this routine in the stochastic sampling

                ! but for open-shell to open-shell exchange excitations
                ! I have to calculate the correct x1 matrix element..

                ! x0 matrix element is easy
                x0 = -occ_i * occ_j / 2.0_dp

                if (inc_i == 1 .and. inc_j == 1) then
                    ! if we have open-shell to open chell
                    x1 = calcDiagExchangeGUGA_nI(i, j, nI) / 2.0_dp

                    call add_to_rdm_spawn_t(spawn, s, p, p, s, &
                                            (x0 - x1) * full_sign, .true.)
                    call add_to_rdm_spawn_t(spawn, p, s, s, p, &
                                            (x0 - x1) * full_sign, .true.)

                else
                    call add_to_rdm_spawn_t(spawn, s, p, p, s, &
                                            x0 * full_sign, .true.)
                    ! and the symmetric version:
                    call add_to_rdm_spawn_t(spawn, p, s, s, p, &
                                            x0 * full_sign, .true.)
                end if

                j = j + inc_j
            end do
            i = i + inc_i
        end do

    end subroutine fill_spawn_rdm_diag_guga