diff --git a/python/triqs_dft_tools/sumk_dft_tools.py b/python/triqs_dft_tools/sumk_dft_tools.py
index 20221abd..44fcd6f8 100644
--- a/python/triqs_dft_tools/sumk_dft_tools.py
+++ b/python/triqs_dft_tools/sumk_dft_tools.py
@@ -53,12 +53,11 @@ class SumkDFTTools(SumkDFT):
                         parproj_data=parproj_data, symmpar_data=symmpar_data, bands_data=bands_data,
                         transp_data=transp_data, misc_data=misc_data, cont_data=cont_data)
 
-    # Uses .data of only GfReFreq objects.
-    # Uses .data of only GfReFreq objects.
     def density_of_states(self, mu=None, broadening=None, mesh=None, with_Sigma=True, with_dc=True, proj_type=None, dosocc=False, save_to_file=True):
         """
-        Calculates the density of states. The basis of the projected density of states is 
-        specified by proj_type.
+        Calculates the density of states and the projected density of states. 
+        The basis of the projected density of states is specified by proj_type.
+
         Parameters
         ----------
         mu           : double, optional
@@ -96,51 +95,53 @@ class SumkDFTTools(SumkDFT):
                        DOS projected to atoms and resolved into orbital contributions.
                        Empty if proj_type = None
         """
-        if(proj_type!=None):
-            #assert proj_type in ('wann', 'vasp','wien2k','elk'), "'proj_type' must be either 'wann', 'vasp', 'wien2k', or 'elk'"
-            assert proj_type in ('wann', 'vasp','wien2k',), "'proj_type' must be either 'wann', 'vasp', 'wien2k'"
-            if(proj_type!='wann'):
-                assert proj_type==self.dft_code, "proj_type must be from the corresponding dft inputs."
+        if (proj_type != None):
+            # assert proj_type in ('wann', 'vasp','wien2k','elk'), "'proj_type' must be either 'wann', 'vasp', 'wien2k', or 'elk'"
+            assert proj_type in ('wann', 'vasp', 'wien2k',
+                                 ), "'proj_type' must be either 'wann', 'vasp', 'wien2k'"
+            if (proj_type != 'wann'):
+                assert proj_type == self.dft_code, "proj_type must be from the corresponding dft inputs."
 
         if (with_Sigma):
-            assert isinstance(self.Sigma_imp[0].mesh, MeshReFreq), "SumkDFT.mesh must be real if with_Sigma is True"
-            mesh=self.Sigma_imp[0].mesh
+            assert isinstance(
+                self.Sigma_imp[0].mesh, MeshReFreq), "SumkDFT.mesh must be real if with_Sigma is True"
+            mesh = self.Sigma_imp[0].mesh
         elif mesh is not None:
             assert isinstance(mesh, MeshReFreq), "mesh must be of form MeshReFreq"
             if broadening is None:
-                broadening=0.001
+                broadening = 0.001
         elif self.mesh is not None:
             assert isinstance(self.mesh, MeshReFreq), "self.mesh must be of form MeshReFreq"
-            mesh=self.mesh
+            mesh = self.mesh
             if broadening is None:
-                broadening=0.001
+                broadening = 0.001
         else:
             assert 0, "ReFreqMesh input required for calculations without real frequency self-energy"
-        mesh_val = numpy.linspace(mesh.w_min,mesh.w_max,len(mesh))
+        mesh_val = numpy.linspace(mesh.w_min, mesh.w_max, len(mesh))
         n_om = len(mesh)
         om_minplot = mesh_val[0] - 0.001
         om_maxplot = mesh_val[-1] + 0.001
 
-        #Read in occupations from HDF5 file if required
-        if(dosocc):
+        # Read in occupations from HDF5 file if required
+        if (dosocc):
             mpi.report('Reading occupations generated by self.occupations().')
             thingstoread = ['occik']
             subgroup_present, values_not_read = self.read_input_from_hdf(
                 subgrp=self.misc_data, things_to_read=thingstoread)
             if len(values_not_read) > 0 and mpi.is_master_node:
-            	raise ValueError(
-                'ERROR: One or more necessary SumK input properties have not been found in the given h5 archive:', self.values_not_read)    
+                raise ValueError(
+                       'ERROR: One or more necessary SumK input properties have not been found in the given h5 archive:', self.values_not_read)
 
-        #initialise projected DOS type if required
+        # initialise projected DOS type if required
         spn = self.spin_block_names[self.SO]
-        n_shells=1
+        n_shells = 1
         if (proj_type == 'wann'):
             n_shells = self.n_corr_shells
             gf_struct = self.gf_struct_sumk.copy()
             dims = [self.corr_shells[ish]['dim'] for ish in range(n_shells)]
             shells_type = 'corr'
         elif (proj_type == 'vasp'):
-            n_shells=1
+            n_shells = 1
             gf_struct = [[(sp, list(range(self.proj_mat_csc.shape[2]))) for sp in spn]]
             dims = [self.proj_mat_csc.shape[2]]
             shells_type = 'csc'
@@ -148,10 +149,10 @@ class SumkDFTTools(SumkDFT):
             self.load_parproj()
             n_shells = self.n_shells
             gf_struct = [[(sp, self.shells[ish]['dim']) for sp in spn]
-                             for ish in range(n_shells)]
+                         for ish in range(n_shells)]
             dims = [self.shells[ish]['dim'] for ish in range(n_shells)]
             shells_type = 'all'
-#	#commented out for now - unsure this produces DFT+DMFT PDOS
+# #commented out for now - unsure this produces DFT+DMFT PDOS
 #        elif (proj_type == 'elk'):
 #            n_shells = self.n_shells
 #            dims = [self.shells[ish]['dim'] for ish in range(n_shells)]
@@ -164,16 +165,16 @@ class SumkDFTTools(SumkDFT):
 #            	raise ValueError(
 #                'ERROR: One or more necessary SumK input properties have not been found in the given h5 archive:', self.values_not_read)
 
-        #set-up output arrays
-        DOS = {sp: numpy.zeros([n_om],float) for sp in spn}
+        # set-up output arrays
+        DOS = {sp: numpy.zeros([n_om], float) for sp in spn}
         DOSproj = [{} for ish in range(n_shells)]
         DOSproj_orb = [{} for ish in range(n_shells)]
-        #set-up Green's function object
+        # set-up Green's function object
         if (proj_type != None):
           G_loc = []
           for ish in range(n_shells):
               glist = [GfReFreq(target_shape=(block_dim, block_dim), mesh=mesh)
-                     for block, block_dim in gf_struct[ish]]
+                       for block, block_dim in gf_struct[ish]]
               G_loc.append(
                   BlockGf(name_list=spn, block_list=glist, make_copies=False))
               G_loc[ish].zero()
@@ -183,17 +184,17 @@ class SumkDFTTools(SumkDFT):
                   DOSproj_orb[ish][sp] = numpy.zeros(
                       [n_om, dim, dim], complex)
 
-        #calculate the DOS
+        # calculate the DOS
         ikarray = numpy.array(list(range(self.n_k)))
         for ik in mpi.slice_array(ikarray):
             G_latt_w = self.lattice_gf(
                 ik=ik, mu=mu, broadening=broadening, mesh=mesh, with_Sigma=with_Sigma, with_dc=with_dc)
             G_latt_w *= self.bz_weights[ik]
-            #output occupied DOS if nk inputted
-            if(dosocc):
+            # output occupied DOS if nk inputted
+            if (dosocc):
               for bname, gf in G_latt_w:
-                G_latt_w[bname].data[:,:,:] *= self.occik[bname][ik]
-            # DOS 
+                G_latt_w[bname].data[:, :, :] *= self.occik[bname][ik]
+            # DOS
             for bname, gf in G_latt_w:
                 DOS[bname] -= gf.data.imag.trace(axis1=1, axis2=2)/numpy.pi
             # Projected DOS:
@@ -209,13 +210,13 @@ class SumkDFTTools(SumkDFT):
         for bname in DOS:
             DOS[bname] = mpi.all_reduce(DOS[bname])
         # Collect data from mpi and put in projected arrays
-        if(proj_type != None):
+        if (proj_type != None):
           for ish in range(n_shells):
               G_loc[ish] << mpi.all_reduce(G_loc[ish])
         # Symmetrize and rotate to local coord. system if needed:
-          if((proj_type!='vasp') and (proj_type!='elk')):
+          if ((proj_type != 'vasp') and (proj_type != 'elk')):
             if self.symm_op != 0:
-                if proj_type=='wann':
+                if proj_type == 'wann':
                     G_loc = self.symmcorr.symmetrize(G_loc)
                 else:
                     G_loc = self.symmpar.symmetrize(G_loc)
@@ -223,13 +224,13 @@ class SumkDFTTools(SumkDFT):
                 for ish in range(n_shells):
                     for bname, gf in G_loc[ish]:
                         G_loc[ish][bname] << self.rotloc(
-                            ish, gf, direction='toLocal',shells=shells_type)
+                            ish, gf, direction='toLocal', shells=shells_type)
         # G_loc can now also be used to look at orbitally-resolved quantities
           for ish in range(n_shells):
             for bname, gf in G_loc[ish]:  # loop over spins
                 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[:,:,:]
+                    :, :, :] += (1.0j*(gf-gf.conjugate().transpose())/2.0/numpy.pi).data[:, :, :]
 
         # Write to files
         if save_to_file and mpi.is_master_node():
@@ -239,9 +240,9 @@ class SumkDFTTools(SumkDFT):
                     f.write("%s    %s\n" % (mesh_val[iom], DOS[sp][iom]))
                 f.close()
                 # Partial
-                if(proj_type!=None):
+                if (proj_type != None):
                   for ish in range(n_shells):
-                    f = open('DOS_' + proj_type +  '_%s_proj%s.dat' % (sp, ish), 'w')
+                    f = open('DOS_' + proj_type + '_%s_proj%s.dat' % (sp, ish), 'w')
                     for iom in range(n_om):
                         f.write("%s    %s\n" %
                                 (mesh_val[iom], DOSproj[ish][sp][iom]))
@@ -249,18 +250,17 @@ class SumkDFTTools(SumkDFT):
                     # Orbitally-resolved
                     for i in range(dims[ish]):
                         for j in range(dims[ish]):
-                            #For Elk with parproj - skip off-diagonal elements
-                            #if(proj_type=='elk') and (i!=j): continue
+                            # For Elk with parproj - skip off-diagonal elements
+                            # if(proj_type=='elk') and (i!=j): continue
                             f = open('DOS_' + proj_type + '_' + sp + '_proj' + str(ish) +
                                      '_' + str(i) + '_' + str(j) + '.dat', 'w')
                             for iom in range(n_om):
                                 f.write("%s    %s    %s\n" % (
-                                    mesh_val[iom], DOSproj_orb[ish][sp][iom, i, j].real,DOSproj_orb[ish][sp][iom, i, j].imag))
+                                    mesh_val[iom], DOSproj_orb[ish][sp][iom, i, j].real, DOSproj_orb[ish][sp][iom, i, j].imag))
                             f.close()
 
         return DOS, DOSproj, DOSproj_orb
 
-
     def proj_type_G_loc(self, G_latt, G_inp, ik, ish, proj_type=None):
         """
         Internal routine which calculates the project Green's function subject to the 
@@ -292,7 +292,7 @@ class SumkDFTTools(SumkDFT):
         if (proj_type == 'wann'):
            for bname, gf in G_proj:
               G_proj[bname] << self.downfold(ik, ish, bname,
-                                G_latt[bname], gf)  # downfolding G
+                                             G_latt[bname], gf)  # downfolding G
         elif (proj_type == 'vasp'):
            for bname, gf in G_latt:
               G_proj[bname] << self.downfold(ik, ish, bname, gf, G_proj[bname], shells='csc')
@@ -302,7 +302,7 @@ class SumkDFTTools(SumkDFT):
               tmp.zero()
               for bname, gf in tmp:
                  tmp[bname] << self.downfold(ik, ish, bname,
-                                G_latt[bname], gf, shells='all', ir=ir)
+                                             G_latt[bname], gf, shells='all', ir=ir)
               G_proj += tmp
 #        elif (proj_type == 'elk'):
 #           dim = self.shells[ish]['dim']
@@ -328,7 +328,7 @@ class SumkDFTTools(SumkDFT):
 
         return G_proj
 
-    def load_parproj(self,data_type=None):
+    def load_parproj(self, data_type=None):
         """
         Internal routine which loads the n_parproj, proj_mat_all, rot_mat_all and 
         rot_mat_all_time_inv from parproj data from .h5 file.
@@ -339,20 +339,20 @@ class SumkDFTTools(SumkDFT):
                     'band' - reads data converted by bands_convert()        
                     None - reads data converted by parproj_convert()
         """
-        #read in the projectors
+        # read in the projectors
         things_to_read = ['n_parproj', 'proj_mat_all']
         if data_type == 'band':
             subgroup_present, values_not_read = self.read_input_from_hdf(
-              subgrp=self.bands_data, things_to_read=things_to_read)
+                subgrp=self.bands_data, things_to_read=things_to_read)
         else:
             subgroup_present, values_not_read = self.read_input_from_hdf(
-              subgrp=self.parproj_data, things_to_read=things_to_read)
+                subgrp=self.parproj_data, things_to_read=things_to_read)
             if self.symm_op:
                 self.symmpar = Symmetry(self.hdf_file, subgroup=self.symmpar_data)
         if len(values_not_read) > 0 and mpi.is_master_node:
             raise ValueError(
                 'ERROR: One or more necessary SumK input properties have not been found in the given h5 archive:', self.values_not_read)
-        #read general data
+        # read general data
         things_to_read = ['rot_mat_all', 'rot_mat_all_time_inv']
         subgroup_present, values_not_read = self.read_input_from_hdf(
             subgrp=self.parproj_data, things_to_read=things_to_read)
@@ -364,6 +364,7 @@ class SumkDFTTools(SumkDFT):
         """
         Calculates the band resolved density matrices (occupations) from the Matsubara 
         frequency self-energy.
+
         Parameters
         ----------
         mu           : double, optional
@@ -376,7 +377,8 @@ class SumkDFTTools(SumkDFT):
         save_occ     : boolean, optional
                        If True, saves the band resolved density matrix in misc_data.
         save_to_file : boolean, optional
-                       If True, text files with the calculated data will be created.
+                       If True, text files with the calculated data will be created.\
+
         Returns
         -------
         occik        : Dict of numpy arrays
@@ -387,7 +389,8 @@ class SumkDFTTools(SumkDFT):
             mesh = self.Sigma_imp[0].mesh
         else:
             mesh = self.mesh
-        assert isinstance(mesh, MeshImFreq), "SumkDFT.mesh must be real if with_Sigma is True or mesh is not given"
+        assert isinstance(
+            mesh, MeshImFreq), "SumkDFT.mesh must be real if with_Sigma is True or mesh is not given"
 
         if mu is None:
             mu = self.chemical_potential
@@ -395,20 +398,21 @@ class SumkDFTTools(SumkDFT):
         spn = self.spin_block_names[self.SO]
         occik = {}
         for sp in spn:
-          #same format as gf.data ndarray  
-          occik[sp] = [numpy.zeros([1, self.n_orbitals[ik, ntoi[sp]], self.n_orbitals[ik, ntoi[sp]]], numpy.double) for ik in range(self.n_k)]
-        #calculate the occupations
+          # same format as gf.data ndarray
+          occik[sp] = [numpy.zeros([1, self.n_orbitals[ik, ntoi[sp]],
+                                   self.n_orbitals[ik, ntoi[sp]]], numpy.double) for ik in range(self.n_k)]
+        # calculate the occupations
         ikarray = numpy.array(range(self.n_k))
         for ik in range(self.n_k):
           G_latt = self.lattice_gf(
               ik=ik, mu=mu, with_Sigma=with_Sigma, with_dc=with_dc)
           for bname, gf in G_latt:
-            occik[bname][ik][0,:,:] = gf.density().real
+            occik[bname][ik][0, :, :] = gf.density().real
         # Collect data from mpi:
         for sp in spn:
           occik[sp] = mpi.all_reduce(occik[sp])
         mpi.barrier()
-        #save to HDF5 file (if specified)
+        # save to HDF5 file (if specified)
         if save_occ and mpi.is_master_node():
           things_to_save_misc = ['occik']
           # Save it to the HDF:
@@ -420,7 +424,6 @@ class SumkDFTTools(SumkDFT):
           del ar
         return occik
 
-    # Uses .data of only GfReFreq objects.
     def spectral_contours(self, mu=None, broadening=None, mesh=None, plot_range=None, FS=True, with_Sigma=True, with_dc=True, proj_type=None, save_to_file=True):
         """
         Calculates the correlated spectral function at the Fermi level (relating to the Fermi
@@ -465,13 +468,13 @@ class SumkDFTTools(SumkDFT):
                        (Correlated) k-resolved spectral function projected to atoms and 
                        resolved into orbital contributions. Empty if proj_type = None
         """
-        if(proj_type!=None):
+        if (proj_type != None):
             assert proj_type in ('wann'), "'proj_type' must be 'wann' if not None"
-        #read in the energy contour energies and projectors
-        things_to_read = ['n_k','bmat','BZ_n_k','BZ_iknr','BZ_vkl',
+        # read in the energy contour energies and projectors
+        things_to_read = ['n_k', 'bmat', 'BZ_n_k', 'BZ_iknr', 'BZ_vkl',
                           'n_orbitals', 'proj_mat', 'hopping']
         subgroup_present, values_not_read = self.read_input_from_hdf(
-          subgrp=self.cont_data, things_to_read=things_to_read)
+            subgrp=self.cont_data, things_to_read=things_to_read)
         if len(values_not_read) > 0 and mpi.is_master_node:
             raise ValueError(
                 'ERROR: One or more necessary SumK input properties have not been found in the given h5 archive:', self.values_not_read)
@@ -479,107 +482,110 @@ class SumkDFTTools(SumkDFT):
         if mu is None:
             mu = self.chemical_potential
         if (with_Sigma):
-            assert isinstance(self.Sigma_imp[0].mesh, MeshReFreq), "SumkDFT.mesh must be real if with_Sigma is True"
-            mesh=self.Sigma_imp[0].mesh
+            assert isinstance(
+                self.Sigma_imp[0].mesh, MeshReFreq), "SumkDFT.mesh must be real if with_Sigma is True"
+            mesh = self.Sigma_imp[0].mesh
         elif mesh is not None:
             assert isinstance(mesh, MeshReFreq), "mesh must be of form MeshReFreq"
             if broadening is None:
-                broadening=0.001
+                broadening = 0.001
         elif self.mesh is not None:
             assert isinstance(self.mesh, MeshReFreq), "self.mesh must be of form MeshReFreq"
-            mesh=self.mesh
+            mesh = self.mesh
             if broadening is None:
-                broadening=0.001
+                broadening = 0.001
         else:
             assert 0, "ReFreqMesh input required for calculations without real frequency self-energy"
-        mesh_val = numpy.linspace(mesh.w_min,mesh.w_max,len(mesh))
+        mesh_val = numpy.linspace(mesh.w_min, mesh.w_max, len(mesh))
         n_om = len(mesh)
         om_minplot = mesh_val[0] - 0.001
         om_maxplot = mesh_val[-1] + 0.001
-        #for Fermi Surface calculations
+        # for Fermi Surface calculations
         if FS:
             dw = abs(mesh_val[1]-mesh_val[0])
-            #ensure that a few frequencies around the Fermi level are included
+            # ensure that a few frequencies around the Fermi level are included
             plot_range = [-2*dw, 2*dw]
             mpi.report('Generated A(k,w) will be evaluted at closest frequency to 0.0 in given mesh ')
         if plot_range is None:
-            n_om = len(mesh_val[(mesh_val > om_minplot)&(mesh_val < om_maxplot)])
-            mesh_val2 = mesh_val[(mesh_val > om_minplot)&(mesh_val < om_maxplot)]
+            n_om = len(mesh_val[(mesh_val > om_minplot) & (mesh_val < om_maxplot)])
+            mesh_val2 = mesh_val[(mesh_val > om_minplot) & (mesh_val < om_maxplot)]
         else:
             om_minplot = plot_range[0]
             om_maxplot = plot_range[1]
-            n_om = len(mesh_val[(mesh_val > om_minplot)&(mesh_val < om_maxplot)])
-            mesh_val2 = mesh_val[(mesh_val > om_minplot)&(mesh_val < om_maxplot)]
-            #\omega ~= 0.0 index for FS file
+            n_om = len(mesh_val[(mesh_val > om_minplot) & (mesh_val < om_maxplot)])
+            mesh_val2 = mesh_val[(mesh_val > om_minplot) & (mesh_val < om_maxplot)]
+            # \omega ~= 0.0 index for FS file
             abs_mesh_val = [abs(i) for i in mesh_val2]
-            jw=[i for i in range(len(abs_mesh_val)) if abs_mesh_val[i] == numpy.min(abs_mesh_val[:])]
+            jw = [i for i in range(len(abs_mesh_val)) if abs_mesh_val[i]
+                  == numpy.min(abs_mesh_val[:])]
+
+        # calculate the spectral functions for the irreducible set of k-points
+        [Akw, pAkw, pAkw_orb] = self.gen_Akw(mu=mu, broadening=broadening, mesh=mesh,
+                                             plot_shift=0.0, plot_range=plot_range,
+                                             shell_list=None, with_Sigma=with_Sigma, with_dc=with_dc,
+                                             proj_type=proj_type)
 
-        #calculate the spectral functions for the irreducible set of k-points
-        [Akw, pAkw, pAkw_orb] = self.gen_Akw(mu=mu, broadening=broadening, mesh=mesh, \
-                                plot_shift=0.0, plot_range=plot_range, \
-                                shell_list=None, with_Sigma=with_Sigma, with_dc=with_dc, \
-                                proj_type=proj_type)
-   
         if save_to_file and mpi.is_master_node():
            spn = self.spin_block_names[self.SO]
            vkc = numpy.zeros(3, float)
-           mesh_val2 = mesh_val[(mesh_val > om_minplot)&(mesh_val < om_maxplot)]
+           mesh_val2 = mesh_val[(mesh_val > om_minplot) & (mesh_val < om_maxplot)]
            if FS:
              n_om = 1
            else:
              n_om = len(mesh_val2)
            for sp in spn:
-                # Open file for storage:
-                for iom in range(n_om):
+               # Open file for storage:
+               for iom in range(n_om):
                   if FS:
                     f = open('Akw_FS_' + sp + '.dat', 'w')
-                    jom=jw[0]
+                    jom = jw[0]
                   else:
                     f = open('Akw_omega_%s_%s.dat' % (iom, sp), 'w')
-                    jom=iom
-                  f.write("#Spectral function evaluated at frequency = %s\n" %mesh_val2[jom])
+                    jom = iom
+                  f.write("#Spectral function evaluated at frequency = %s\n" % mesh_val2[jom])
                   for ik in range(self.BZ_n_k):
-                     jk=self.BZ_iknr[ik]
-                     vkc[:] = numpy.matmul(self.bmat,self.BZ_vkl[ik,:])
+                     jk = self.BZ_iknr[ik]
+                     vkc[:] = numpy.matmul(self.bmat, self.BZ_vkl[ik, :])
                      f.write("%s    %s    %s    %s\n" % (vkc[0], vkc[1], vkc[2], Akw[sp][jk, jom]))
                   f.close()
-           if (proj_type!=None):
-                n_shells = len(pAkw[:])
-                for iom in range(n_om):
+           if (proj_type != None):
+               n_shells = len(pAkw[:])
+               for iom in range(n_om):
                   for sp in spn:
                     for ish in range(n_shells):
                       if FS:
                         strng = 'Akw_FS' + '_' + proj_type + '_' + sp + '_proj' + str(ish)
-                        jom=jw[0]
+                        jom = jw[0]
                       else:
-                        strng = 'Akw_omega_' + str(iom) + '_' + proj_type + '_' + sp + '_proj' + str(ish)
-                        jom=iom
+                        strng = 'Akw_omega_' + str(iom) + '_' + proj_type + \
+                            '_' + sp + '_proj' + str(ish)
+                        jom = iom
                       f = open(strng + '.dat', 'w')
-                      f.write("#Spectral function evaluated at frequency = %s\n" %mesh_val2[jom])
+                      f.write("#Spectral function evaluated at frequency = %s\n" % mesh_val2[jom])
                       for ik in range(self.BZ_n_k):
-                        jk=self.BZ_iknr[ik]
-                        vkc[:] = numpy.matmul(self.bmat,self.BZ_vkl[ik,:])
+                        jk = self.BZ_iknr[ik]
+                        vkc[:] = numpy.matmul(self.bmat, self.BZ_vkl[ik, :])
                         f.write("%s    %s    %s    %s\n" % (vkc[0], vkc[1], vkc[2],
-                                 pAkw[ish][sp][jk, jom]))
+                                                            pAkw[ish][sp][jk, jom]))
                       f.close()
-                      dim=len(pAkw_orb[ish][sp][0, 0, 0, :])
+                      dim = len(pAkw_orb[ish][sp][0, 0, 0, :])
                       for i in range(dim):
                         for j in range(dim):
-                        #For Elk with parproj - skip off-diagonal elements
-                          if(proj_type=='elk') and (i!=j): continue
+                            # For Elk with parproj - skip off-diagonal elements
+                          if (proj_type == 'elk') and (i != j):
+                              continue
                           strng2 = strng + '_' + str(i) + '_' + str(j)
                         # Open file for storage:
                           f = open(strng2 + '.dat', 'w')
                           for ik in range(self.BZ_n_k):
-                             jk=self.BZ_iknr[ik]
-                             vkc[:] = numpy.matmul(self.bmat,self.BZ_vkl[ik,:])
+                             jk = self.BZ_iknr[ik]
+                             vkc[:] = numpy.matmul(self.bmat, self.BZ_vkl[ik, :])
                              f.write("%s    %s    %s    %s\n" % (vkc[0], vkc[1], vkc[2],
-                                    pAkw_orb[ish][sp][jk, jom, i, j]))
+                                                                 pAkw_orb[ish][sp][jk, jom, i, j]))
                           f.close()
 
         return Akw, pAkw, pAkw_orb
 
-    # Uses .data of only GfReFreq objects.
     def spaghettis(self, mu=None, broadening=None, mesh=None, plot_shift=0.0, plot_range=None, shell_list=None, with_Sigma=True, with_dc=True, proj_type=None, save_to_file=True):
         """
         Calculates the k-resolved spectral function A(k,w) (band structure)
@@ -597,12 +603,12 @@ class SumkDFTTools(SumkDFT):
         plot_shift   : double, optional
                        Offset for each A(k,w) for stacked plotting of spectra.
         plot_range   : list of double, optional
-                       Sets the energy window for plotting to (plot_range[0],plot_range[1]). 
-                       If not provided, the min and max values of the energy mesh is used.           
+                       Sets the energy window for plotting to (plot_range[0],plot_range[1]).
+                       If not provided, the min and max values of the energy mesh is used.
         shell_list   : list of integers, optional
-                       Contains the indices of the shells of which the projected spectral function 
-                       is calculated for. 
-                       If shell_list = None and proj_type is not None, then the projected spectral 
+                       Contains the indices of the shells of which the projected spectral function
+                       is calculated for.
+                       If shell_list = None and proj_type is not None, then the projected spectral
                        function is calculated for all shells.
         with_Sigma   : boolean, optional
                        If True, the self energy is used for the calculation.
@@ -620,44 +626,45 @@ class SumkDFTTools(SumkDFT):
         Akw          : Dict of numpy arrays
                        (Correlated) k-resolved spectral function
         pAkw      : Dict of numpy arrays
-                       (Correlated) k-resolved spectral function projected to atoms. 
+                       (Correlated) k-resolved spectral function projected to atoms.
                        Empty if proj_type = None
         pAkw_orb  : Dict of numpy arrays
-                       (Correlated) k-resolved spectral function projected to atoms and 
+                       (Correlated) k-resolved spectral function projected to atoms and
                        resolved into orbital contributions. Empty if proj_type = None
         """
 
-        #initialisation
-        if(proj_type!=None):
+        # initialisation
+        if (proj_type != None):
             assert proj_type in ('wann', 'wien2k'), "'proj_type' must be either 'wann', 'wien2k'"
-            if(proj_type!='wann'):
-                assert proj_type==self.dft_code, "proj_type must be from the corresponding dft inputs."
+            if (proj_type != 'wann'):
+                assert proj_type == self.dft_code, "proj_type must be from the corresponding dft inputs."
         things_to_read = ['n_k', 'n_orbitals', 'proj_mat', 'hopping']
         subgroup_present, values_not_read = self.read_input_from_hdf(
             subgrp=self.bands_data, things_to_read=things_to_read)
         if len(values_not_read) > 0 and mpi.is_master_node:
             raise ValueError(
                 'ERROR: One or more necessary SumK input properties have not been found in the given h5 archive:', self.values_not_read)
-        if(proj_type=='wien2k'):
+        if (proj_type == 'wien2k'):
            self.load_parproj(data_type='band')
 
         if mu is None:
             mu = self.chemical_potential
         if (with_Sigma):
-            assert isinstance(self.Sigma_imp[0].mesh, MeshReFreq), "SumkDFT.mesh must be real if with_Sigma is True"
-            mesh=self.Sigma_imp[0].mesh
+            assert isinstance(
+                self.Sigma_imp[0].mesh, MeshReFreq), "SumkDFT.mesh must be real if with_Sigma is True"
+            mesh = self.Sigma_imp[0].mesh
         elif mesh is not None:
             assert isinstance(mesh, MeshReFreq), "mesh must be of form MeshReFreq"
             if broadening is None:
-                broadening=0.001
+                broadening = 0.001
         elif self.mesh is not None:
             assert isinstance(self.mesh, MeshReFreq), "self.mesh must be of form MeshReFreq"
-            mesh=self.mesh
+            mesh = self.mesh
             if broadening is None:
-                broadening=0.001
+                broadening = 0.001
         else:
             assert 0, "ReFreqMesh input required for calculations without real frequency self-energy"
-        mesh_val = numpy.linspace(mesh.w_min,mesh.w_max,len(mesh))
+        mesh_val = numpy.linspace(mesh.w_min, mesh.w_max, len(mesh))
         n_om = len(mesh)
         om_minplot = mesh_val[0] - 0.001
         om_maxplot = mesh_val[-1] + 0.001
@@ -667,52 +674,53 @@ class SumkDFTTools(SumkDFT):
         else:
             om_minplot = plot_range[0]
             om_maxplot = plot_range[1]
-        n_om = len(mesh_val[(mesh_val > om_minplot)&(mesh_val < om_maxplot)])
+        n_om = len(mesh_val[(mesh_val > om_minplot) & (mesh_val < om_maxplot)])
 
-        [Akw, pAkw, pAkw_orb] = self.gen_Akw(mu=mu, broadening=broadening, mesh=mesh, \
-                                plot_shift=plot_shift, plot_range=plot_range, \
-                                shell_list=shell_list, with_Sigma=with_Sigma, with_dc=with_dc, \
-                                proj_type=proj_type)
+        [Akw, pAkw, pAkw_orb] = self.gen_Akw(mu=mu, broadening=broadening, mesh=mesh,
+                                             plot_shift=plot_shift, plot_range=plot_range,
+                                             shell_list=shell_list, with_Sigma=with_Sigma, with_dc=with_dc,
+                                             proj_type=proj_type)
 
         if save_to_file and mpi.is_master_node():
-           mesh_val2 = mesh_val[(mesh_val > om_minplot)&(mesh_val < om_maxplot)]
-           spn = self.spin_block_names[self.SO]
-           for sp in spn:
+            mesh_val2 = mesh_val[(mesh_val > om_minplot) & (mesh_val < om_maxplot)]
+            spn = self.spin_block_names[self.SO]
+            for sp in spn:
                 # Open file for storage:
                 f = open('Akw_' + sp + '.dat', 'w')
                 for ik in range(self.n_k):
                     for iom in range(n_om):
-                        f.write('%s     %s      %s\n' %(ik, mesh_val2[iom], Akw[sp][ik, iom]))
+                        f.write('%s     %s      %s\n' % (ik, mesh_val2[iom], Akw[sp][ik, iom]))
                     f.write('\n')
                 f.close()
-           if (proj_type!=None):
+            if (proj_type != None):
                 n_shells = len(pAkw[:])
-                if shell_list==None:
-                  shell_list=[ish for ish in range(n_shells)]
+                if shell_list == None:
+                  shell_list = [ish for ish in range(n_shells)]
                 for sp in spn:
                   for ish in range(n_shells):
-                    jsh=shell_list[ish]
+                    jsh = shell_list[ish]
                     f = open('Akw_' + proj_type + '_' +
-                                 sp + '_proj' + str(jsh) + '.dat', 'w')
+                             sp + '_proj' + str(jsh) + '.dat', 'w')
                     for ik in range(self.n_k):
                        for iom in range(n_om):
                           f.write('%s     %s      %s\n' % (
-                                 ik, mesh_val2[iom], pAkw[ish][sp][ik, iom]))
+                              ik, mesh_val2[iom], pAkw[ish][sp][ik, iom]))
                        f.write('\n')
                     f.close()
-                    #get orbital dimension from the length of dimension of the array
-                    dim=len(pAkw_orb[ish][sp][0, 0, 0, :])
+                    # get orbital dimension from the length of dimension of the array
+                    dim = len(pAkw_orb[ish][sp][0, 0, 0, :])
                     for i in range(dim):
                       for j in range(dim):
-                        #For Elk with parproj - skip off-diagonal elements
-                        if(proj_type=='elk') and (i!=j): continue
+                        # For Elk with parproj - skip off-diagonal elements
+                        if(proj_type =='elk') and (i!=j):
+                            continue
                         # Open file for storage:
                         f = open('Akw_' + proj_type + '_' + sp + '_proj' + str(jsh)
-                                + '_' + str(i) + '_' + str(j) + '.dat', 'w')
+                                 + '_' + str(i) + '_' + str(j) + '.dat', 'w')
                         for ik in range(self.n_k):
                             for iom in range(n_om):
                                 f.write('%s     %s      %s\n' % (
-                                            ik, mesh_val2[iom], pAkw_orb[ish][sp][ik, iom, i, j]))
+                                    ik, mesh_val2[iom], pAkw_orb[ish][sp][ik, iom, i, j]))
                             f.write('\n')
                         f.close()