calc_explicit_1_rdm_guga Subroutine

    subroutine calc_explicit_1_rdm_guga(ilut, csf_i, n_tot, excitations)
        ! routine to calculate the one-RDM explicitly in the GUGA formalism
        ! the one RDM is given by rho_ij = <Psi|E_ij|Psi>
        ! so, similar to the actHamiltonian routine I need to loop over all
        ! the orbital indices (i,j) and calculate the action of a given E_ij
        ! on the sampled wavefunction |Psi> and the resulting overlap with <Psi|
        ! make this routine similar to GenExcDjs() in rdm_explicit.F90
        ! to insert it there to calculate the GUGA RDMs in this case
        integer(n_int), intent(in) :: ilut(0:GugaBits%len_tot)
        type(CSF_Info_t), intent(in) :: csf_i
        integer, intent(out) :: n_tot
        integer(n_int), intent(out), allocatable :: excitations(:, :)
        character(*), parameter :: this_routine = "calc_explicit_1_rdm_guga"

        integer :: i, j, nMax, ierr, n, n_excits
        integer(n_int), allocatable :: temp_excits(:, :), tmp_all_excits(:, :)

        nMax = 6 + 4 * (nSpatOrbs)**2 * (count_open_orbs(ilut) + 1)
        allocate(tmp_all_excits(0:GugaBits%len_tot, nMax), stat=ierr)
        call LogMemAlloc('tmp_all_excits', (GugaBits%len_tot + 1) * nMax, 8, this_routine, tag_tmp_excits)

        n_tot = 0

        do i = 1, nSpatOrbs
            do j = 1, nSpatOrbs
                if (i == j) cycle

                call calc_all_excits_guga_rdm_singles(ilut, csf_i, i, j, temp_excits, &
                                                      n_excits)

#ifdef DEBUG_
                do n = 1, n_excits
                    ASSERT(isProperCSF_ilut(temp_excits(:, n), .true.))
                end do
#endif

                if (n_excits > 0) then
                    call add_guga_lists_rdm(n_tot, n_excits, tmp_all_excits, temp_excits)
                end if

                deallocate(temp_excits)

            end do
        end do

        j = 1
        do i = 1, n_tot
            if (near_zero(extract_matrix_element(tmp_all_excits(:, i), 1))) cycle

            tmp_all_excits(:, j) = tmp_all_excits(:, i)

            j = j + 1

        end do

        n_tot = j - 1

        allocate(excitations(0:GugaBits%len_tot, n_tot), &
                  source=tmp_all_excits(0:GugaBits%len_tot, 1:n_tot), stat=ierr)
        ! hm to log that does not make so much sense.. since it gets called
        ! more than once and is only a temporary array..
        call LogMemAlloc('excitations', n_tot, 8, this_routine, tag_excitations)

        call sort(excitations(:, 1:n_tot), ilut_lt, ilut_gt)

        deallocate(tmp_all_excits)
        call LogMemDealloc(this_routine, tag_tmp_excits)

    end subroutine calc_explicit_1_rdm_guga