diff --git a/doc/Ce-gamma.py b/doc/Ce-gamma.py
index 8d42d848..104753a7 100644
--- a/doc/Ce-gamma.py
+++ b/doc/Ce-gamma.py
@@ -70,7 +70,7 @@ for Iteration_Number in range(1,Loops+1):
 
         # Compute the SumK, possibly fixing mu by dichotomy
         if SK.density_required and (Iteration_Number > 0):
-            Chemical_potential = SK.find_mu( precision = 0.01 )
+            Chemical_potential = SK.calc_mu( precision = 0.01 )
         else:
             mpi.report("No adjustment of chemical potential\nTotal density  = %.3f"%SK.total_density(mu=Chemical_potential))
 
@@ -80,7 +80,7 @@ for Iteration_Number in range(1,Loops+1):
         dm = S.G.density()
 
         if ((Iteration_Number==1)and(previous_present==False)):
-	    SK.set_dc( dens_mat=dm, U_interact = U_int, J_hund = J_hund, orb = 0, use_dc_formula = DC_type, use_val=use_val)
+	    SK.calc_dc( dens_mat=dm, U_interact = U_int, J_hund = J_hund, orb = 0, use_dc_formula = DC_type, use_val=use_val)
 
         # set atomic levels:
         eal = SK.eff_atomic_levels()[0]
@@ -119,7 +119,7 @@ for Iteration_Number in range(1,Loops+1):
         
         # after the Solver has finished, set new double counting: 
         dm = S.G.density()
-        SK.set_dc( dm, U_interact = U_int, J_hund = J_hund, orb = 0, use_dc_formula = DC_type , use_val=use_val)
+        SK.calc_dc( dm, U_interact = U_int, J_hund = J_hund, orb = 0, use_dc_formula = DC_type , use_val=use_val)
         # correlation energy calculations:
         correnerg = 0.5 * (S.G * S.Sigma).total_density()
         mpi.report("Corr. energy = %s"%correnerg)
@@ -153,7 +153,7 @@ for Iteration_Number in range(1,Loops+1):
 
 # find exact chemical potential
 if (SK.density_required):
-    SK.chemical_potential = SK.find_mu( precision = 0.000001 )
+    SK.chemical_potential = SK.calc_mu( precision = 0.000001 )
 dN,d = SK.calc_density_correction(filename = dft_filename+'.qdmft')
 
 mpi.report("Trace of Density Matrix: %s"%d)
diff --git a/doc/DFTDMFTmain.rst b/doc/DFTDMFTmain.rst
index 4727a2ca..e3e08b89 100644
--- a/doc/DFTDMFTmain.rst
+++ b/doc/DFTDMFTmain.rst
@@ -92,22 +92,22 @@ set up the loop over DMFT iterations and the self-consistency condition::
   for iteration_number in range(n_loops) :            # start the DMFT loop
 
           SK.put_Sigma(Sigma_imp = [ S.Sigma ])      # Put self energy to the SumK class
-          chemical_potential = SK.find_mu()          # find the chemical potential for the given density
-          S.G << SK.extract_G_loc()[0]              # extract the local Green function
-          S.G0 << inverse(S.Sigma + inverse(S.G))   # finally get G0, the input for the Solver
+          chemical_potential = SK.calc_mu()          # calculate the chemical potential for the given density
+          S.G << SK.extract_G_loc()[0]               # extract the local Green function
+          S.G0 << inverse(S.Sigma + inverse(S.G))    # finally get G0, the input for the Solver
 
           S.solve(U_interact,J_hund,use_spinflip=False,use_matrix=True,     # now solve the impurity problem
                            l=2,T=None, dim_reps=None, irep=None, n_cycles =10000,
                            length_cycle=200,n_warmup_cycles=1000)
 
 	  dm = S.G.density()                                                 # Density matrix of the impurity problem  
-          SK.set_dc( dm, U_interact = U, J_hund = J, use_dc_formula = 0)     # Set the double counting term
+          SK.calc_dc( dm, U_interact = U, J_hund = J, use_dc_formula = 0)     # Set the double counting term
           SK.save(['chemical_potential','dc_imp','dc_energ'])                # Save data in the hdf5 archive
 
 These basic steps are enough to set up the basic DMFT Loop. For a detailed description of the :class:`SumkDFT` routines,
 see the reference manual. After the self-consistency steps, the solution of the Anderson impurity problem is calculation by CTQMC. 
 Different to model calculations, we have to do a few more steps after this, because of the double-counting correction. We first 
-calculate the density of the impurity problem. Then, the routine `set_dc` takes as parameters this density matrix, the 
+calculate the density of the impurity problem. Then, the routine `calc_dc` takes as parameters this density matrix, the 
 Coulomb interaction, Hund's rule coupling, and the type of double-counting that should be used. Possible values for `use_dc_formula` are:
 
   * `0`: Full-localised limit
diff --git a/doc/advanced.rst b/doc/advanced.rst
index 0ed99ea2..475bef98 100644
--- a/doc/advanced.rst
+++ b/doc/advanced.rst
@@ -89,14 +89,14 @@ previous section, with some additional refinement::
         SK.symm_deg_gf(S.Sigma,orb=0)                           # symmetrise Sigma
         SK.put_Sigma(Sigma_imp = [ S.Sigma ])                   # put Sigma into the SumK class:
   
-        chemical_potential = SK.find_mu( precision = prec_mu )  # find the chemical potential
+        chemical_potential = SK.calc_mu( precision = prec_mu )  # find the chemical potential
         S.G << SK.extract_G_loc()[0]                           # calculation of the local Green function
         mpi.report("Total charge of Gloc : %.6f"%S.G.total_density())
   
         if ((iteration_number==1)and(previous_present==False)):
             # Init the DC term and the real part of Sigma, if no previous run was found:
             dm = S.G.density()
-            SK.set_dc( dm, U_interact = U, J_hund = J, orb = 0, use_dc_formula = dc_type)
+            SK.calc_dc( dm, U_interact = U, J_hund = J, orb = 0, use_dc_formula = dc_type)
             S.Sigma << SK.dc_imp[0]['up'][0,0]
   
         S.G0 << inverse(S.Sigma + inverse(S.G))
@@ -141,7 +141,7 @@ previous section, with some additional refinement::
 
         # Now set new double counting:
         dm = S.G.density()
-        SK.set_dc( dm, U_interact = U, J_hund = J, orb = 0, use_dc_formula = dc_type)
+        SK.calc_dc( dm, U_interact = U, J_hund = J, orb = 0, use_dc_formula = dc_type)
   
         #Save stuff:
         SK.save(['chemical_potential','dc_imp','dc_energ'])
diff --git a/doc/selfcons.rst b/doc/selfcons.rst
index 255399fb..51af05f2 100644
--- a/doc/selfcons.rst
+++ b/doc/selfcons.rst
@@ -34,7 +34,7 @@ Now we calculate the modified charge density::
 
   # find exact chemical potential
   SK.put_Sigma(Sigma_imp = [ S.Sigma ])
-  chemical_potential = SK.find_mu( precision = 0.000001 )
+  chemical_potential = SK.calc_mu( precision = 0.000001 )
   dN, d = SK.calc_density_correction(filename = dft_filename+'.qdmft')
   SK.save(['chemical_potential','dc_imp','dc_energ'])
 
diff --git a/lda_dmft_cthyb.py b/lda_dmft_cthyb.py
index e3ead752..d332ac0d 100644
--- a/lda_dmft_cthyb.py
+++ b/lda_dmft_cthyb.py
@@ -73,14 +73,14 @@ for iteration_number in range(1,loops+1):
 
       SK.symm_deg_gf(S.Sigma_iw,orb=0)                        # symmetrise Sigma
       SK.put_Sigma(Sigma_imp = [ S.Sigma_iw ])                # put Sigma into the SumK class
-      chemical_potential = SK.find_mu( precision = prec_mu )  # find the chemical potential for the given density
+      chemical_potential = SK.calc_mu( precision = prec_mu )  # find the chemical potential for the given density
       S.G_iw << SK.extract_G_loc()[0]                           # extract the local Green function
       mpi.report("Total charge of Gloc : %.6f"%S.G_iw.total_density())
 
       if ((iteration_number==1)and(previous_present==False)):
           # Init the DC term and the real part of Sigma, if no previous run was found:
           dm = S.G_iw.density()
-          SK.set_dc(dm, U_interact = U, J_hund = J, orb = 0, use_dc_formula = dc_type)
+          SK.calc_dc(dm, U_interact = U, J_hund = J, orb = 0, use_dc_formula = dc_type)
           S.Sigma_iw << SK.dc_imp[0]['up'][0,0]
 
       # now calculate new G0_iw to input into the solver:
@@ -126,7 +126,7 @@ for iteration_number in range(1,loops+1):
 
       dm = S.G_iw.density() # compute the density matrix of the impurity problem
       # Set the double counting
-      SK.set_dc( dm, U_interact = U, J_hund = J, orb = 0, use_dc_formula = dc_type)
+      SK.calc_dc( dm, U_interact = U, J_hund = J, orb = 0, use_dc_formula = dc_type)
 
       # Save stuff into the dft_output group of hdf5 archive in case of rerun:
       SK.save(['chemical_potential','dc_imp','dc_energ'])
diff --git a/python/U_matrix.py b/python/U_matrix.py
index 67156cb4..ebb203d8 100644
--- a/python/U_matrix.py
+++ b/python/U_matrix.py
@@ -3,7 +3,6 @@ from scipy.misc import factorial as fact
 from itertools import product
 import numpy as np
 
-
 # The interaction matrix in desired basis
 # U^{spherical}_{m1 m2 m3 m4} = \sum_{k=0}^{2l} F_k angular_matrix_element(l, k, m1, m2, m3, m4)
 def U_matrix(l, radial_integrals=None, U_int=None, J_hund=None, basis="spherical", T=None):
@@ -207,7 +206,6 @@ def clebsch_gordan(jm1, jm2, jm3):
 #  column 0 for 1st dim,
 #  columns 2 and 3 for 2nd dim,
 #  columns 0,1,2 and 3 for 3rd dim.
-#def subarray(a,idxlist,n=len(a.shape)-1) :
 def subarray(a,idxlist,n=None) :
     if n is None: n = len(a.shape)-1
     sa = a[tuple(slice(x) for x in a.shape[:n]) + (idxlist[n],)]
diff --git a/python/clear_h5_output.py b/python/clear_h5_output.py
index 17eb0ccb..8edb1775 100644
--- a/python/clear_h5_output.py
+++ b/python/clear_h5_output.py
@@ -13,7 +13,7 @@ and to restore it to the original post-converter state.
 
 filename = sys.argv[1]
 A = h5py.File(filename)
-for group in ['dft_output','dmft_output']:
+for group in ['dft_output','user_data']:
     if group in A: del(A[group])
 A.close()
 
diff --git a/python/sumk_dft.py b/python/sumk_dft.py
index c7c8fd8b..eba5f267 100644
--- a/python/sumk_dft.py
+++ b/python/sumk_dft.py
@@ -32,9 +32,9 @@ class SumkDFT:
     """This class provides a general SumK method for combining ab-initio code and pytriqs."""
 
 
-    def __init__(self, hdf_file, mu = 0.0, h_field = 0.0, use_dft_blocks = False, 
+    def __init__(self, hdf_file, h_field = 0.0, use_dft_blocks = False, 
 	         dft_data = 'dft_input', symmcorr_data = 'dft_symmcorr_input', parproj_data = 'dft_parproj_input', 
-                 symmpar_data = 'dft_symmpar_input', bands_data = 'dft_bands_input', dft_output = 'dft_output'):
+                 symmpar_data = 'dft_symmpar_input', bands_data = 'dft_bands_input'):
         """
         Initialises the class from data previously stored into an HDF5
         """
@@ -48,16 +48,16 @@ class SumkDFT:
             self.parproj_data = parproj_data
             self.symmpar_data = symmpar_data
             self.bands_data = bands_data
-            self.dft_output = dft_output
             self.G_upfold = None
             self.h_field = h_field
 
-            # read input from HDF:
+            # Read input from HDF:
             things_to_read = ['energy_unit','n_k','k_dep_projection','SP','SO','charge_below','density_required',
                               'symm_op','n_shells','shells','n_corr_shells','corr_shells','use_rotations','rot_mat',
                               'rot_mat_time_inv','n_reps','dim_reps','T','n_orbitals','proj_mat','bz_weights','hopping',
                               'n_inequiv_shells', 'corr_to_inequiv', 'inequiv_to_corr']
             self.subgroup_present, self.value_read = self.read_input_from_hdf(subgrp = self.dft_data, things_to_read = things_to_read)
+            if self.symm_op: self.symmcorr = Symmetry(hdf_file,subgroup=self.symmcorr_data)
 
             if self.SO and (abs(self.h_field) > 0.000001):
                 self.h_field = 0.0
@@ -65,7 +65,7 @@ class SumkDFT:
 
             self.spin_block_names = [ ['up','down'], ['ud'] ]
             self.n_spin_blocks = [2,1]
-            # convert spin_block_names to indices -- if spin polarized, differentiate up and down blocks 
+            # Convert spin_block_names to indices -- if spin polarized, differentiate up and down blocks 
             self.spin_names_to_ind = [{}, {}]
             for iso in range(2): # SO = 0 or 1
                 for isp in range(self.n_spin_blocks[iso]):
@@ -75,50 +75,33 @@ class SumkDFT:
             # Most general form allowing for all hybridisation, i.e. largest blocks possible
             self.gf_struct_sumk = [ [ (sp, range( self.corr_shells[icrsh]['dim'])) for sp in self.spin_block_names[self.corr_shells[icrsh]['SO']] ]
                                    for icrsh in range(self.n_corr_shells) ]
+            # First set a standard gf_struct solver:
+            self.gf_struct_solver = [ dict([ (sp, range(self.corr_shells[self.inequiv_to_corr[ish]]['dim']) )
+                                                    for sp in self.spin_block_names[self.corr_shells[self.inequiv_to_corr[ish]]['SO']] ])
+                                      for ish in range(self.n_inequiv_shells) ]
+            # Set standard (identity) maps from gf_struct_sumk <-> gf_struct_solver
+            self.sumk_to_solver = [ {} for ish in range(self.n_inequiv_shells) ]
+            self.solver_to_sumk = [ {} for ish in range(self.n_inequiv_shells) ]
+            self.solver_to_sumk_block = [ {} for ish in range(self.n_inequiv_shells) ]
+            for ish in range(self.n_inequiv_shells):
+                for block,inner_list in self.gf_struct_sumk[self.inequiv_to_corr[ish]]:
+                    self.solver_to_sumk_block[ish][block] = block
+                    for inner in inner_list:
+                        self.sumk_to_solver[ish][(block,inner)] = (block,inner)
+                        self.solver_to_sumk[ish][(block,inner)] = (block,inner)
+            self.deg_shells = [ [] for ish in range(self.n_inequiv_shells) ] # assume no shells are degenerate
 
-            #-----
-            # If these quantities are not in HDF, set them up
-            optional_things = ['gf_struct_solver','sumk_to_solver','solver_to_sumk','solver_to_sumk_block','chemical_potential','dc_imp','dc_energ','deg_shells']
-            self.subgroup_present, self.value_read = self.read_input_from_hdf(subgrp = self.dft_output, things_to_read = [], 
-                                                                              optional_things = optional_things)
-            if (not self.subgroup_present) or (not self.value_read['gf_struct_solver']):
-                # No gf_struct was stored in HDF, so first set a standard one:
-                self.gf_struct_solver = [ dict([ (sp, range(self.corr_shells[self.inequiv_to_corr[ish]]['dim']) )
-                                                        for sp in self.spin_block_names[self.corr_shells[self.inequiv_to_corr[ish]]['SO']] ])
-                                          for ish in range(self.n_inequiv_shells)
-                                        ]
-                # Set standard (identity) maps from gf_struct_sumk <-> gf_struct_solver
-                self.sumk_to_solver = [ {} for ish in range(self.n_inequiv_shells) ]
-                self.solver_to_sumk = [ {} for ish in range(self.n_inequiv_shells) ]
-                self.solver_to_sumk_block = [ {} for ish in range(self.n_inequiv_shells) ]
-                for ish in range(self.n_inequiv_shells):
-                    for block,inner_list in self.gf_struct_sumk[self.inequiv_to_corr[ish]]:
-                        self.solver_to_sumk_block[ish][block] = block
-                        for inner in inner_list:
-                            self.sumk_to_solver[ish][(block,inner)] = (block,inner)
-                            self.solver_to_sumk[ish][(block,inner)] = (block,inner)
+            self.chemical_potential = 0.0 # initialise mu
+            self.init_dc() # initialise the double counting
 
-            if (not self.subgroup_present) or (not self.value_read['dc_imp']):
-                self.__init_dc() # initialise the double counting
-
-            if (not self.subgroup_present) or (not self.value_read['chemical_potential']):
-                self.chemical_potential = mu
-
-            if (not self.subgroup_present) or (not self.value_read['deg_shells']):
-                self.deg_shells = [ [] for ish in range(self.n_inequiv_shells)]
-            #-----
-
-            if self.symm_op:
-                self.symmcorr = Symmetry(hdf_file,subgroup=self.symmcorr_data)
-
-            # Analyse the block structure and determine the smallest blocks, if desired
-            if use_dft_blocks: dm = self.analyse_block_structure()
+            # Analyse the block structure and determine the smallest gf_struct blocks and maps, if desired
+            if use_dft_blocks: self.analyse_block_structure()
 
 ################
 # HDF5 FUNCTIONS
 ################
 
-    def read_input_from_hdf(self, subgrp, things_to_read=[], optional_things=[]):
+    def read_input_from_hdf(self, subgrp, things_to_read):
         """
         Reads data from the HDF file
         """
@@ -126,7 +109,6 @@ class SumkDFT:
         value_read = True
         # initialise variables on all nodes to ensure mpi broadcast works at the end
         for it in things_to_read: setattr(self,it,0)
-        for it in optional_things: setattr(self,it,0)
         subgroup_present = 0
 
         if mpi.is_master_node():
@@ -140,16 +122,6 @@ class SumkDFT:
                     else:
                         mpi.report("Loading %s failed!"%it)
                         value_read = False
-
-                if value_read and (len(optional_things) > 0):
-                    # if successfully read necessary items, read optional things:
-                    value_read = {}
-                    for it in optional_things:
-                        if it in ar[subgrp]:
-                            setattr(self,it,ar[subgrp][it])
-                            value_read['%s'%it] = True
-                        else:
-                            value_read['%s'%it] = False
             else:
                 if (len(things_to_read) != 0): mpi.report("Loading failed: No %s subgroup in HDF5!"%subgrp)
                 subgroup_present = False
@@ -158,26 +130,41 @@ class SumkDFT:
             del ar
         # now do the broadcasting:
         for it in things_to_read: setattr(self,it,mpi.bcast(getattr(self,it)))
-        for it in optional_things: setattr(self,it,mpi.bcast(getattr(self,it)))
         subgroup_present = mpi.bcast(subgroup_present)
         value_read = mpi.bcast(value_read)
 
         return subgroup_present, value_read
 
 
-    def save(self,things_to_save):
-        """Saves given quantities into the 'dft_output' subgroup of the HDF5 archive"""
+    def save(self, things_to_save, subgrp='user_data'):
+        """Saves given quantities into the subgroup ('user_data' by default) of the HDF5 archive"""
 
         if not (mpi.is_master_node()): return # do nothing on nodes
         ar = HDFArchive(self.hdf_file,'a')
-        if not self.dft_output in ar: ar.create_group(self.dft_output)
+        if not subgrp in ar: ar.create_group(subgrp)
         for it in things_to_save: 
             try:
-                ar[self.dft_output][it] = getattr(self,it)
+                ar[subgrp][it] = getattr(self,it)
             except:
-                mpi.report("%s not found, and so not stored."%it)
+                mpi.report("%s not found, and so not saved."%it)
         del ar
 
+
+    def load(self, things_to_load, subgrp='user_data'):
+        """Loads given quantities from the subgroup ('user_data' by default) of the HDF5 archive"""
+
+        if not (mpi.is_master_node()): return # do nothing on nodes
+        ar = HDFArchive(self.hdf_file,'a')
+        if not subgrp in ar: mpi.report("Loading %s failed!"%subgrp)
+        list_to_return = []
+        for it in things_to_load: 
+            try:
+                list_to_return.append(ar[subgrp][it])
+            except:
+                raise ValueError, "%s not found, and so not loaded."%it 
+        del ar
+        return list_to_return
+
 ################
 # CORE FUNCTIONS
 ################
@@ -431,7 +418,7 @@ class SumkDFT:
                         self.solver_to_sumk[ish][(block_solv,inner_solv)] = (block_sumk,inner_sumk)
                         self.solver_to_sumk_block[ish][block_solv] = block_sumk
 
-            # now calculate degeneracies of orbitals:
+            # Now calculate degeneracies of orbitals
             dm = {}
             for block,inner in self.gf_struct_solver[ish].iteritems():
                 # get dm for the blocks:
@@ -446,9 +433,9 @@ class SumkDFT:
                 for block2 in self.gf_struct_solver[ish].iterkeys(): 
                     if dm[block1].shape == dm[block2].shape:
                         if ( (abs(dm[block1] - dm[block2]) < threshold).all() ) and (block1 != block2):
-                            # check if it was already there:
                             ind1 = -1
                             ind2 = -2
+                            # check if it was already there:
                             for n,ind in enumerate(self.deg_shells[ish]):
                                 if block1 in ind: ind1 = n
                                 if block2 in ind: ind2 = n
@@ -459,11 +446,6 @@ class SumkDFT:
                             elif (ind1 < 0) and (ind2 < 0):
                                 self.deg_shells[ish].append([block1,block2])
 
-        things_to_save = ['gf_struct_solver','sumk_to_solver','solver_to_sumk','solver_to_sumk_block','deg_shells']
-        self.save(things_to_save)
-
-        return dens_mat
-
 
     def density_matrix(self, method = 'using_gf', beta = 40.0):
         """Calculate density matrices in one of two ways:
@@ -506,6 +488,8 @@ class SumkDFT:
                         else:
                             MMat[isp][inu,inu] = 0.0
 
+            else: raise ValueError, "density_matrix: the method '%s' is not supported."%method
+
             for icrsh in range(self.n_corr_shells):
                 for isp, sp in enumerate(self.spin_block_names[self.corr_shells[icrsh]['SO']]):
                     isp = self.spin_names_to_ind[self.corr_shells[icrsh]['SO']][sp]
@@ -550,7 +534,8 @@ class SumkDFT:
 
         # Chemical Potential:
         for ish in range(self.n_inequiv_shells):
-            for ii in eff_atlevels[ish]: eff_atlevels[ish][ii] *= -self.chemical_potential
+            for ii in eff_atlevels[ish]: 
+                eff_atlevels[ish][ii] *= -self.chemical_potential
 
         # double counting term:
         for ish in range(self.n_inequiv_shells):
@@ -599,9 +584,8 @@ class SumkDFT:
         return eff_atlevels
 
 
-    def __init_dc(self):
-
-        # construct the density matrix dm_imp and double counting arrays
+    def init_dc(self):
+        """ Initialise the double counting terms to have the correct shape."""
         self.dc_imp = [ {} for icrsh in range(self.n_corr_shells)]
         for icrsh in range(self.n_corr_shells):
             dim = self.corr_shells[icrsh]['dim']
@@ -610,46 +594,45 @@ class SumkDFT:
         self.dc_energ = [0.0 for icrsh in range(self.n_corr_shells)]
 
 
-    def set_dc(self,dens_mat,U_interact,J_hund,orb=0,use_dc_formula=0,use_val=None):
-        """Sets the double counting term for inequiv orbital orb:
+    def set_dc(self,dc_imp,dc_energ):
+        """Sets double counting terms dc_imp and dc_energ to known values."""
+
+        self.dc_imp = dc_imp
+        self.dc_energ = dc_energ
+
+
+    def calc_dc(self,dens_mat,orb=0,U_interact=None,J_hund=None,use_dc_formula=0,use_dc_value=None):
+        """Sets the double counting corrections in the correct form for inequiv orbital orb:
+           1) either using U_interact, J_hund and 
            use_dc_formula = 0: fully-localised limit (FLL),
            use_dc_formula = 1: Held's formula,
            use_dc_formula = 2: around mean-field (AMF).
+           2) or using a given dc value in use_dc_value.
            Be sure that you are using the correct interaction Hamiltonian!"""
 
         for icrsh in range(self.n_corr_shells):
 
-            iorb = self.corr_to_inequiv[icrsh]    # iorb is the index of the inequivalent shell corresponding to icrsh
-
-            if iorb != orb: continue # ignore this orbital
-
-            Ncr = {}
+            ish = self.corr_to_inequiv[icrsh]    # ish is the index of the inequivalent shell corresponding to icrsh
+            if ish != orb: continue # ignore this orbital
             dim = self.corr_shells[icrsh]['dim'] #*(1+self.corr_shells[icrsh]['SO'])
-
             spn = self.spin_block_names[self.corr_shells[icrsh]['SO']]
-            for sp in spn:
-                self.dc_imp[icrsh][sp] = numpy.identity(dim,numpy.float_)
-                Ncr[sp] = 0.0
 
-            for block,inner in self.gf_struct_solver[iorb].iteritems():
-                bl = self.solver_to_sumk_block[iorb][block]
+            Ncr = { sp: 0.0 for sp in spn }
+            for block,inner in self.gf_struct_solver[ish].iteritems():
+                bl = self.solver_to_sumk_block[ish][block]
                 Ncr[bl] += dens_mat[block].real.trace()
-
-            Ncrtot = 0.0
-
-            spn = self.spin_block_names[self.corr_shells[icrsh]['SO']]
+            Ncrtot = sum(Ncr.itervalues())
             for sp in spn: 
-                Ncrtot += Ncr[sp]
-
-            # average the densities if there is no SP:
-            if self.SP == 0:
-                for sp in spn: Ncr[sp] = Ncrtot / len(spn)
-            # correction for SO: we have only one block in this case, but in DC we need N/2
-            elif self.SP == 1 and self.SO == 1:
-                for sp in spn: Ncr[sp] = Ncrtot / 2.0
+                self.dc_imp[icrsh][sp] = numpy.identity(dim,numpy.float_)
+                if self.SP == 0: # average the densities if there is no SP:
+                    Ncr[sp] = Ncrtot / len(spn)
+                elif self.SP == 1 and self.SO == 1: # correction for SO: we have only one block in this case, but in DC we need N/2
+                    Ncr[sp] = Ncrtot / 2.0
 
             if use_val is None:
 
+                if U_interact is None and J_hund is None: raise ValueError, "set_dc: either provide U_interact and J_hund or set use_val to dc value."
+
                 if use_dc_formula == 0: # FLL
 
                     self.dc_energ[icrsh] = U_interact / 2.0 * Ncrtot * (Ncrtot-1.0)
@@ -676,16 +659,13 @@ class SumkDFT:
                         self.dc_energ[icrsh] -= (U_interact + (dim-1)*J_hund)/dim * 0.5 * Ncr[sp] * Ncr[sp]
                         mpi.report("DC for shell %(icrsh)i and block %(sp)s = %(Uav)f"%locals())
 
-                # output:
                 mpi.report("DC energy for shell %s = %s"%(icrsh,self.dc_energ[icrsh]))
 
-            else:
+            else: # use value provided for user to determine dc_energ and dc_imp
 
                 self.dc_energ[icrsh] = use_val * Ncrtot
-                for sp in spn:
-                    self.dc_imp[icrsh][sp] *= use_val
+                for sp in spn: self.dc_imp[icrsh][sp] *= use_val
 
-                # output:
                 mpi.report("DC for shell %(icrsh)i = %(use_val)f"%locals())
                 mpi.report("DC energy = %s"%self.dc_energ[icrsh])
 
@@ -745,7 +725,7 @@ class SumkDFT:
         self.chemical_potential = mu
 
 
-    def find_mu(self, precision = 0.01):
+    def calc_mu(self, precision = 0.01):
         """
         Searches for mu in order to give the desired charge
         A desired precision can be specified in precision.
@@ -840,8 +820,8 @@ class SumkDFT:
 # FIXME LEAVE UNDOCUMENTED
 ################
 
-    # FIXME Merge with find_mu?
-    def find_mu_nonint(self, dens_req, orb = None, precision = 0.01):
+    # FIXME Merge with calc_mu?
+    def calc_mu_nonint(self, dens_req, orb = None, precision = 0.01):
 
         def F(mu):
             gnonint = self.extract_G_loc(mu = mu, with_Sigma = False)
@@ -873,7 +853,7 @@ class SumkDFT:
         mu = self.chemical_potential
 
         def F(dc):
-            self.set_dc(dens_mat = dens_mat, U_interact = 0, J_hund = 0, orb = orb, use_val = dc)
+            self.calc_dc(dens_mat = dens_mat, U_interact = 0, J_hund = 0, orb = orb, use_val = dc)
             if dens_req is None:
                 return self.total_density(mu = mu)
             else:
diff --git a/python/sumk_dft_tools.py b/python/sumk_dft_tools.py
index bb5dffb0..a9b0c6f8 100644
--- a/python/sumk_dft_tools.py
+++ b/python/sumk_dft_tools.py
@@ -32,11 +32,11 @@ class SumkDFTTools(SumkDFT):
     """Extends the SumkDFT class with some tools for analysing the data."""
 
 
-    def __init__(self, hdf_file, mu = 0.0, h_field = 0.0, use_dft_blocks = False, dft_data = 'dft_input', symmcorr_data = 'dft_symmcorr_input',
+    def __init__(self, hdf_file, h_field = 0.0, use_dft_blocks = False, dft_data = 'dft_input', symmcorr_data = 'dft_symmcorr_input',
                  parproj_data = 'dft_parproj_input', symmpar_data = 'dft_symmpar_input', bands_data = 'dft_bands_input'):
 
         self.G_upfold_refreq = None
-        SumkDFT.__init__(self, hdf_file=hdf_file, mu=mu, h_field=h_field, use_dft_blocks=use_dft_blocks,
+        SumkDFT.__init__(self, hdf_file=hdf_file, h_field=h_field, use_dft_blocks=use_dft_blocks,
                           dft_data=dft_data, symmcorr_data=symmcorr_data, parproj_data=parproj_data, 
                           symmpar_data=symmpar_data, bands_data=bands_data)
 
diff --git a/python/update_archive.py b/python/update_archive.py
index 936a3783..e25f7071 100644
--- a/python/update_archive.py
+++ b/python/update_archive.py
@@ -62,15 +62,18 @@ for old, new in old_to_new.iteritems():
     A.copy(old,new)
     del(A[old])
 
-# Move output items from dft_input to dft_output
-move_to_output = ['gf_struct_solver','map_inv','map',
-                  'chemical_potential','dc_imp','dc_energ','deg_shells',
-                  'h_field']
+# Move output items from dft_input to user_data
+move_to_output = ['chemical_potential','dc_imp','dc_energ']
 for obj in move_to_output:
     if obj in A['dft_input'].keys():
-       if not 'dft_output' in A: A.create_group('dft_output')
-       print "Moving %s to dft_output ..."%obj
-       A.copy('dft_input/'+obj,'dft_output/'+obj)
+       if 'user_data' not in A: A.create_group('user_data')
+       print "Moving %s to user_data ..."%obj
+       A.copy('dft_input/'+obj,'user_data/'+obj)
+       del(A['dft_input'][obj])
+# Delete obsolete quantities
+to_delete = ['gf_struct_solver','map_inv','map','deg_shells','h_field']
+for obj in to_delete:
+    if obj in A['dft_input'].keys():
        del(A['dft_input'][obj])
 
 # Update shells and corr_shells to list of dicts
@@ -87,16 +90,18 @@ HDFArchive(filename,'a')['dft_input']['corr_shells'] = corr_shells
 A = h5py.File(filename)
 
 # Add shell equivalency quantities
-equiv_shell_info = det_shell_equivalence(corr_shells)
-A['dft_input']['n_inequiv_shells'] = equiv_shell_info[0]
-A['dft_input']['corr_to_inequiv'] = equiv_shell_info[1]
-A['dft_input']['inequiv_to_corr'] = equiv_shell_info[2]
+if 'n_inequiv_shells' not in A['dft_input']:
+    equiv_shell_info = det_shell_equivalence(corr_shells)
+    A['dft_input']['n_inequiv_shells'] = equiv_shell_info[0]
+    A['dft_input']['corr_to_inequiv'] = equiv_shell_info[1]
+    A['dft_input']['inequiv_to_corr'] = equiv_shell_info[2]
 
 # Rename variables
 groups = ['dft_symmcorr_input','dft_symmpar_input']
 
 for group in groups:
     if group not in A.keys(): continue
+    if 'n_s' not in group: continue
     print "Changing n_s to n_symm ..."
     A[group].move('n_s','n_symm')
     # Convert orbits to list of dicts
diff --git a/test/srvo3_Gloc.output.h5 b/test/srvo3_Gloc.output.h5
index c33d0ab4..9d30e79c 100644
Binary files a/test/srvo3_Gloc.output.h5 and b/test/srvo3_Gloc.output.h5 differ