calc_explicit_2_rdm_guga Subroutine

    subroutine calc_explicit_2_rdm_guga(ilut, csf_i, n_tot, excitations)
        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_2_rdm_guga"

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

        nMax = 6 + 4 * (nSpatOrbs)**3 * (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
                do k = 1, nSpatOrbs
                    do l = 1, nSpatOrbs

                        ! pure diagonal contribution:
                        if (i == j .and. k == l) cycle

                        call calc_all_excits_guga_rdm_doubles(ilut, csf_i, i, j, k, l, &
                                                              temp_excits, n_excits)

#ifdef DEBUG_
                        do n = 1, n_excits
                            if (.not. isProperCSF_ilut(temp_excits(:, n), .true.)) then
                                print *, "===="
                                call write_det_guga(stdout, ilut, .true.)
                                call write_det_guga(stdout, temp_excits(:, n), .true.)
                                print *, i, j, k, l
                            end if
                            ASSERT(isProperCSF_ilut(temp_excits(:, n), .true.))
                        end do
#endif
                        if (i == l .and. j == k) then
                            ! exclude the diagonal exchange here,
                            ! as it is already accounted for in the
                            ! diagonal contribution routine
#ifdef DEBUG_
                            if (n_excits > 0) then
                                if (.not. DetBitEQ(ilut, temp_excits(:, 1), nifd)) then
                                    print *, "not equal!"
                                end if
                            end if
#endif

                            if (n_excits > 1) then
                                call add_guga_lists_rdm(n_tot, n_excits - 1, &
                                                        tmp_all_excits, temp_excits(:, 2:))
                            end if
                        else
                            if (n_excits > 0) then
                                call add_guga_lists_rdm(n_tot, n_excits, tmp_all_excits, temp_excits)
                            end if
                        end if

                        deallocate(temp_excits)

                    end do
                end do
            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)

        deallocate(tmp_all_excits)
        call LogMemDealloc(this_routine, tag_tmp_excits)

    end subroutine calc_explicit_2_rdm_guga