add_rdm_from_ij_pair_guga_exact Subroutine

    subroutine add_rdm_from_ij_pair_guga_exact(spawn, one_rdms, ilutI, csf_i, ilutJ, csf_j, &
                                               sign_i, sign_j)
        ! this routine is called for RDM sampling within the semi-stochastic
        ! space or for the connection to the 'true' HF det!
        ! i also need the calc-type input, as in the semi-stochastic space
        ! we initialize the csf info of nI, while for the HF connections we dont
        type(rdm_spawn_t), intent(inout) :: spawn
        type(one_rdm_t), intent(inout) :: one_rdms(:)
        integer(n_int), intent(in) :: ilutI(0:IlutBits%len_tot), &
                                      ilutJ(0:IlutBits%len_tot)
        type(CSF_Info_t), intent(in) :: csf_i, csf_j
        real(dp), intent(in) :: sign_i(:), sign_j(:)
        character(*), parameter :: this_routine = "add_rdm_from_ij_pair_guga_exact"
        integer(int_rdm), allocatable :: rdm_ind(:)
        real(dp), allocatable :: rdm_mat(:)
        type(ExcitationInformation_t) :: excitInfo
        HElement_t(dp) :: mat_ele
        integer :: p, q, r, s, n
        real(dp) :: full_sign(spawn%rdm_send%sign_length)

        call calc_guga_matrix_element(ilutI, csf_i, ilutJ, csf_j, excitInfo, mat_ele, &
                                      t_hamil=.false., rdm_ind=rdm_ind, rdm_mat=rdm_mat)

        ! i assume sign_i and sign_j are not 0 if we end up here..
        if (allocated(rdm_ind)) then
            ASSERT(allocated(rdm_mat))
            ASSERT(size(rdm_ind) == size(rdm_mat))
            do n = 1, size(rdm_ind)
                if (.not. near_zero(rdm_mat(n))) then
                    if (excitInfo%excitLvl == 1) then
                        if (RDMExcitLevel == 1) then
                            call fill_sings_1rdm_guga(one_rdms, sign_i, sign_j, &
                                                      rdm_mat(n), rdm_ind(n))
                        else
                            call fill_sings_2rdm_guga(spawn, ilutI, &
                                                      ilutJ, sign_i, sign_j, &
                                                      rdm_mat(n), rdm_ind(n))
                        end if
                    else if (excitInfo%excitLvl == 2 .and. RDMExcitLevel /= 1) then
                        call extract_2_rdm_ind(rdm_ind(n), p, q, r, s)
                        full_sign = sign_i * sign_j * rdm_mat(n)
                        ! here in the 'exact' filling (coming from HF or
                        ! within the semistochastic space I think it makes
                        ! sense to fill symmetrically.. since here no
                        ! stochastic spawning is happening and this does not
                        ! give us information about the hermiticity error!
                        call add_to_rdm_spawn_t(spawn, p, q, r, s, &
                                                full_sign, .true.)
                        if (.not. &
                            (excitInfo%typ == excit_type%fullstart_stop_alike)) then
                            call add_to_rdm_spawn_t(spawn, r, s, p, q, &
                                                    full_sign, .true.)
                        end if
                    end if
                end if
            end do
        end if

    end subroutine add_rdm_from_ij_pair_guga_exact