diff --git a/python/triqs_dft_tools/sumk_dft_tools.py b/python/triqs_dft_tools/sumk_dft_tools.py
index fdcae04f..c40e9838 100644
--- a/python/triqs_dft_tools/sumk_dft_tools.py
+++ b/python/triqs_dft_tools/sumk_dft_tools.py
@@ -153,9 +153,7 @@ class SumkDFTTools(SumkDFT):
         # G_loc can now also be used to look at orbitally-resolved quantities
         for ish in range(self.n_inequiv_shells):
             for bname, gf in G_loc[self.inequiv_to_corr[ish]]:  # loop over spins
-                for iom in range(n_om):
-                    DOSproj[ish][bname][iom] -= gf.data[iom,
-                                                        :, :].imag.trace() / numpy.pi
+                DOSproj[ish][bname] = -gf.data.imag.trace(axis1=1, axis2=2) / numpy.pi
                 DOSproj_orb[ish][bname][
                     :, :, :] += (1.0j*(gf-gf.conjugate().transpose())/2.0/numpy.pi).data[:,:,:]
 
@@ -279,8 +277,7 @@ class SumkDFTTools(SumkDFT):
 
         # G_loc can now also be used to look at orbitally-resolved quantities
         for bname, gf in G_loc_all:  # loop over spins
-            for iom in range(n_om):
-                DOSproj[bname][iom] -= gf.data[iom,:,:].imag.trace() / numpy.pi
+            DOSproj[bname] = -gf.data.imag.trace(axis1=1, axis2=2) / numpy.pi
             DOSproj_orb[bname][:,:,:] += (1.0j*(gf-gf.conjugate().transpose())/2.0/numpy.pi).data[:,:,:]
         # Write to files
         if save_to_file and mpi.is_master_node():
@@ -392,10 +389,8 @@ class SumkDFTTools(SumkDFT):
             G_latt_w *= self.bz_weights[ik]
 
             # Non-projected DOS
-            for iom in range(n_om):
-                for bname, gf in G_latt_w:
-                    DOS[bname][iom] -= gf.data[iom, :, :].imag.trace() / \
-                        numpy.pi
+            for bname, gf in G_latt_w:
+                DOS[bname] -= gf.data.imag.trace(axis1=1, axis2=2) / numpy.pi
 
             # Projected DOS:
             for ish in range(self.n_shells):
@@ -427,9 +422,7 @@ class SumkDFTTools(SumkDFT):
         # G_loc can now also be used to look at orbitally-resolved quantities
         for ish in range(self.n_shells):
             for bname, gf in G_loc[ish]:
-                for iom in range(n_om):
-                    DOSproj[ish][bname][iom] -= gf.data[iom,
-                                                        :, :].imag.trace() / numpy.pi
+                DOSproj[ish][bname] = -gf.data.imag.trace(axis1=1, axis2=2) / numpy.pi
                 DOSproj_orb[ish][bname][
                     :, :, :] += (1.0j*(gf-gf.conjugate().transpose())/2.0/numpy.pi).data[:,:,:]
 
diff --git a/test/python/SrVO3_spectral.h5 b/test/python/SrVO3_spectral.h5
new file mode 100644
index 00000000..455e5962
Binary files /dev/null and b/test/python/SrVO3_spectral.h5 differ
diff --git a/test/python/srvo3_spectral.py b/test/python/srvo3_spectral.py
new file mode 100644
index 00000000..63b833a7
--- /dev/null
+++ b/test/python/srvo3_spectral.py
@@ -0,0 +1,36 @@
+from triqs_dft_tools.sumk_dft_tools import *
+from triqs.utility.h5diff import h5diff 
+from h5 import HDFArchive
+
+beta = 40
+
+SK = SumkDFTTools(hdf_file='SrVO3_spectral.h5', use_dft_blocks=True)
+
+if mpi.is_master_node():
+    with HDFArchive('SrVO3_Sigma.h5', 'a') as ar:
+        Sigma = ar['dmft_transp_input']['Sigma_w']
+        SK.chemical_potential = ar['dmft_transp_input']['chemical_potential']
+        SK.dc_imp = ar['dmft_transp_input']['dc_imp']
+
+Sigma = mpi.bcast(Sigma)
+SK.chemical_potential = mpi.bcast(SK.chemical_potential)
+SK.dc_imp = mpi.bcast(SK.dc_imp)
+SK.set_Sigma([Sigma])
+
+dos_wannier = SK.dos_wannier_basis(broadening=0.01, with_Sigma=True, with_dc=True, save_to_file=False)
+dos_parproj = SK.dos_parproj_basis(broadening=0.01, with_Sigma=True, with_dc=True, save_to_file=False)
+spaghetti = SK.spaghettis(broadening=0.01, plot_shift=0.0, plot_range=(-1,1), ishell=None, save_to_file=False)
+
+if mpi.is_master_node():
+
+    # with HDFArchive('srvo3_spectral.ref.h5', 'a') as ar:
+    #     ar['dos_wannier'] = dos_wannier 
+    #     ar['dos_parproj'] = dos_parproj
+    #     ar['spaghetti'] = spaghetti
+
+    with HDFArchive('srvo3_spectral.out.h5', 'a') as ar:
+        ar['dos_wannier'] = dos_wannier 
+        ar['dos_parproj'] = dos_parproj
+        ar['spaghetti'] = spaghetti
+
+    h5diff('srvo3_spectral.out.h5', 'srvo3_spectral.ref.h5')
diff --git a/test/python/srvo3_spectral.ref.h5 b/test/python/srvo3_spectral.ref.h5
new file mode 100644
index 00000000..9d12fd40
Binary files /dev/null and b/test/python/srvo3_spectral.ref.h5 differ