diff --git a/doc/Ce-gamma.py b/doc/Ce-gamma.py
index a60a999e..8d42d848 100644
--- a/doc/Ce-gamma.py
+++ b/doc/Ce-gamma.py
@@ -44,8 +44,8 @@ previous_present = mpi.bcast(previous_present)
 # Init the SumK class
 SK=SumkDFT(hdf_file=dft_filename+'.h5',use_dft_blocks=False)
 
-Norb = SK.corr_shells[0][3]
-l    = SK.corr_shells[0][2]
+Norb = SK.corr_shells[0]['dim']
+l    = SK.corr_shells[0]['l']
 
 # Init the Solver:
 S = Solver(beta = beta, l = l)
diff --git a/doc/Ce-gamma_DOS.py b/doc/Ce-gamma_DOS.py
index 573084c3..b2975361 100644
--- a/doc/Ce-gamma_DOS.py
+++ b/doc/Ce-gamma_DOS.py
@@ -54,8 +54,8 @@ SK = SumkDFTTools(hdf_file=dft_filename+'.h5',use_dft_blocks=False)
 if (mpi.is_master_node()):
     print 'DC after reading SK: ',SK.dc_imp[SK.invshellmap[0]]
 
-N = SK.corr_shells[0][3]
-l = SK.corr_shells[0][2]
+N = SK.corr_shells[0]['dim']
+l = SK.corr_shells[0]['l']
 
 # Init the Solver:
 S = Solver(beta = Beta, l = l)
diff --git a/doc/advanced.rst b/doc/advanced.rst
index 10725193..0ed99ea2 100644
--- a/doc/advanced.rst
+++ b/doc/advanced.rst
@@ -63,8 +63,8 @@ Now we can use all this information to initialise the :class:`SumkDFT` class::
 If there was a previous run, we know already about the block structure, and therefore `UseDFTBlocs` is set to `False`.
 The next step is to initialise the Solver::
 
-  Norb = SK.corr_shells[0][3]
-  l = SK.corr_shells[0][2]
+  Norb = SK.corr_shells[0]['dim']
+  l = SK.corr_shells[0]['l']
   S = SolverMultiBand(beta=beta,n_orb=Norb,gf_struct=SK.gf_struct_solver[0],map=SK.map[0])
 
 As we can see, many options of the solver are set by properties of the :class:`SumkDFT` class, so we don't have
diff --git a/doc/reference/h5structure.rst b/doc/reference/h5structure.rst
index 7ba99eef..76ff3a04 100644
--- a/doc/reference/h5structure.rst
+++ b/doc/reference/h5structure.rst
@@ -31,11 +31,11 @@ In order to be used with the DMFT routines, the following data needs to be provi
 
   * `n_shells`, numpy.int. Number of atomic shells for which post processing is possible. This is `not` the number of correlated orbitals! If you have two correlated atoms in the unit cell that are equivalent, `n_shells` is 2!
  
-  * `shells`, double list of numpy.int. First dimension: `n_shells`, second dimension: 4. Information about the atomic shells. For each shell, we give a list of 4 numbers: [index, sort, l, dim]. index is the atom index, sort defines the equivalency of the atoms. For instance, with two equivalent atoms in the unit cell, index runs from 0 to 1, but sort can take only one value 0. l is the angular quantum number, and dim the dimension of the atomic shell.
+  * `shells`, list of dictionaries {string:int}. List dimension: `n_shells`, dictionary dimension: 4. Information about the atomic shells. For each shell, we give a dict with keys ['atom', 'sort', 'l', 'dim']. atom is the atom index, sort defines the equivalency of the atoms. For instance, with two equivalent atoms in the unit cell, atom runs from 0 to 1, but sort can take only one value 0. l is the angular quantum number, and dim the dimension of the atomic shell.
 
   * `n_corr_shells`, numpy.int. Number of correlated atomic shells, for which correlations are included. This includes atoms which are equivalent by symmetry. If you have two correlated atoms in the unit cell that are equivalent, `n_corr_shells` is 2! 
 
-  * `corr_shells`, double list of numpy.int. First dimension: `n_corr_shells`, second dimension: 6. Information about the correlated orbitals. For each correlated shell, we give a list of 6 numbers: [index, sort, l, dim, SO, irep]. Similar as for `shells`, index is the atom index, sort defines the equivalency of the atoms. For instance, with two equivalent atoms in the unit cell, index runs from 0 to 1, but sort can take only one value 0. l is the angular quantum number, and dim the dimension of the atomic shell. If spin-orbit is included in the calculation, SO is 1, otherwise 0. irep is a dummy integer, set to 0.
+  * `corr_shells`, list of dictionaries {string:int}. List dimension: `n_corr_shells`, dictionary dimension: 6. Information about the correlated orbitals. For each correlated shell, we give a dict with keys ['atom', 'sort', 'l', 'dim', 'SO', 'irep']. As for `shells`, atom is the atom index, sort defines the equivalency of the atoms. For instance, with two equivalent atoms in the unit cell, atom runs from 0 to 1, but sort can take only one value 0. l is the angular quantum number, and dim the dimension of the atomic shell. If spin-orbit is included in the calculation, SO is 1, otherwise 0. irep is a dummy integer, set to 0.
 
   * `use_rotations`, numpy.int. If local and global coordinate systems are used, this falg is set to 1. Otherwise set to 0.
 
@@ -81,11 +81,11 @@ The above described converter of the Wien2k input is quite involved, since Wien2
 
   1                <- n_shells
 
-  1 1 2 3          <- shells, as above: iatom, isort, l, dim
+  1 1 2 3          <- shells, as above: atom, sort, l, dim
 
   1                <- n_corr_shells
 
-  1 1 2 3 0 0      <- corr_shells, as above: iatom, isort, l, dim, SO, dummy
+  1 1 2 3 0 0      <- corr_shells, as above: atom, sort, l, dim, SO, dummy
 
   2 2 3            <- n_reps, dim_reps (length 2, because eg/t2g splitting)
 
diff --git a/python/converters/converter_tools.py b/python/converters/converter_tools.py
index 6f64483d..ff5633bd 100644
--- a/python/converters/converter_tools.py
+++ b/python/converters/converter_tools.py
@@ -51,29 +51,32 @@ class ConverterTools:
             subprocess.call(["mv","-f","temphgfrt.h5","%s"%self.hdf_file])
             
 
-    def det_shell_equivalence(self,lst):
+    def det_shell_equivalence(self,corr_shells):
         """
-        The number of inequivalent shells is determined from lst, and a mapping is given as
-        corr_to_inequiv(i_corr_shells) = i_inequiv_corr_shells
-        inequiv_to_corr(i_inequiv_corr_shells) = i_corr_shells
+        The number of inequivalent shells is determined from corr_shells, and a mapping is given as
+        corr_to_inequiv(i_corr_shells) = i_inequiv_shells
+        inequiv_to_corr(i_inequiv_shells) = i_corr_shells
         in order to put the self energies to all equivalent shells, and for extracting Gloc
         """
 
-        corr_to_inequiv = [0 for i in range(len(lst))]
+        corr_to_inequiv = [0 for i in range(len(corr_shells))]
         inequiv_to_corr = [0]
         n_inequiv_shells = 1
-        tmp = [ lst[0][1:3] ]
-        if (len(lst)>1):
-            for i in range(len(lst)-1):
-                fnd = False
+
+        if len(corr_shells) > 1:
+            inequiv_sort = [ corr_shells[0]['sort'] ]
+            inequiv_l = [ corr_shells[0]['l'] ]
+            for i in range(len(corr_shells)-1):
+                is_equiv = False
                 for j in range(n_inequiv_shells):
-                    if (tmp[j]==lst[i+1][1:3]):
-                        fnd = True
+                    if (inequiv_sort[j]==corr_shells[i+1]['sort']) and (inequiv_l[j]==corr_shells[i+1]['l']):
+                        is_equiv = True
                         corr_to_inequiv[i+1] = j
-                if (fnd==False):
+                if is_equiv==False:
                     corr_to_inequiv[i+1] = n_inequiv_shells
                     n_inequiv_shells += 1
-                    tmp.append( lst[i+1][1:3] )
-                    inequiv_to_corr.append(i+1)
+                    inequiv_sort.append( corr_shells[i+1]['sort'] )
+                    inequiv_l.append( corr_shells[i+1]['l'] )
+                    inequiv_to_corr.append( i+1 )
 
-        return [n_inequiv_shells, corr_to_inequiv, inequiv_to_corr]
+        return n_inequiv_shells, corr_to_inequiv, inequiv_to_corr
diff --git a/python/converters/hk_converter.py b/python/converters/hk_converter.py
index 7c7ef3c5..14ad5f8e 100644
--- a/python/converters/hk_converter.py
+++ b/python/converters/hk_converter.py
@@ -74,36 +74,38 @@ class HkConverter(ConverterTools):
             # the information on the non-correlated shells is needed for defining dimension of matrices:
             n_shells = int(R.next())                      # number of shells considered in the Wanniers
                                                           # corresponds to index R in formulas
-            # now read the information about the shells:
-            shells = [ [ int(R.next()) for i in range(4) ] for icrsh in range(n_shells) ]    # reads iatom, sort, l, dim
+            # now read the information about the shells (atom, sort, l, dim):
+            shell_entries = ['atom', 'sort', 'l', 'dim']
+            shells = [ {name: int(val) for name, val in zip(shell_entries, R)} for ish in range(n_shells) ]
 
             n_corr_shells = int(R.next())                 # number of corr. shells (e.g. Fe d, Ce f) in the unit cell, 
                                                           # corresponds to index R in formulas
-            # now read the information about the shells:
-            corr_shells = [ [ int(R.next()) for i in range(6) ] for icrsh in range(n_corr_shells) ]    # reads iatom, sort, l, dim, SO flag, irep
+            # now read the information about the shells (atom, sort, l, dim, SO flag, irep):
+            corr_shell_entries = ['atom', 'sort', 'l', 'dim', 'SO', 'irep']
+            corr_shells = [ {name: int(val) for name, val in zip(corr_shell_entries, R)} for icrsh in range(n_corr_shells) ]
 
             # determine the number of inequivalent correlated shells and maps, needed for further reading
             [n_inequiv_shells, corr_to_inequiv, inequiv_to_corr] = ConverterTools.det_shell_equivalence(self,corr_shells)
 
             use_rotations = 0
-            rot_mat = [numpy.identity(corr_shells[icrsh][3],numpy.complex_) for icrsh in xrange(n_corr_shells)]
+            rot_mat = [numpy.identity(corr_shells[icrsh]['dim'],numpy.complex_) for icrsh in range(n_corr_shells)]
             rot_mat_time_inv = [0 for i in range(n_corr_shells)]
             
             # Representative representations are read from file
             n_reps = [1 for i in range(n_inequiv_shells)]
             dim_reps = [0 for i in range(n_inequiv_shells)]
             T = []
-            for icrsh in range(n_inequiv_shells):
-                n_reps[icrsh] = int(R.next())   # number of representatives ("subsets"), e.g. t2g and eg
-                dim_reps[icrsh] = [int(R.next()) for i in range(n_reps[icrsh])]   # dimensions of the subsets
+            for ish in range(n_inequiv_shells):
+                n_reps[ish] = int(R.next())   # number of representatives ("subsets"), e.g. t2g and eg
+                dim_reps[ish] = [int(R.next()) for i in range(n_reps[ish])]   # dimensions of the subsets
             
                 # The transformation matrix:
                 # is of dimension 2l+1, it is taken to be standard d (as in Wien2k)
-                ll = 2*corr_shells[inequiv_to_corr[icrsh]][2]+1
-                lmax = ll * (corr_shells[inequiv_to_corr[icrsh]][4] + 1)
+                ll = 2*corr_shells[inequiv_to_corr[ish]]['l']+1
+                lmax = ll * (corr_shells[inequiv_to_corr[ish]]['SO'] + 1)
                 T.append(numpy.zeros([lmax,lmax],numpy.complex_))
                 
-                T[icrsh] = numpy.array([[0.0, 0.0, 1.0, 0.0, 0.0],
+                T[ish] = numpy.array([[0.0, 0.0, 1.0, 0.0, 0.0],
                                        [1.0/sqrt(2.0), 0.0, 0.0, 0.0, 1.0/sqrt(2.0)],
                                        [-1.0/sqrt(2.0), 0.0, 0.0, 0.0, 1.0/sqrt(2.0)],
                                        [0.0, 1.0/sqrt(2.0), 0.0, -1.0/sqrt(2.0), 0.0],
@@ -113,28 +115,27 @@ class HkConverter(ConverterTools):
             n_spin_blocs = SP + 1 - SO   # number of spins to read for Norbs and Ham, NOT Projectors
         
             # define the number of n_orbitals for all k points: it is the number of total bands and independent of k!
-            n_orb = sum([ shells[ish][3] for ish in range(n_shells) ])
-            n_orbitals = numpy.ones([n_k,n_spin_blocs],numpy.int) * n_orb
+            n_orbitals = numpy.ones([n_k,n_spin_blocs],numpy.int) * sum([ sh['dim'] for sh in shells ])
 
             # Initialise the projectors:
-            proj_mat = numpy.zeros([n_k,n_spin_blocs,n_corr_shells,max(numpy.array(corr_shells)[:,3]),max(n_orbitals)],numpy.complex_)
+            proj_mat = numpy.zeros([n_k,n_spin_blocs,n_corr_shells,max([crsh['dim'] for crsh in corr_shells]),max(n_orbitals)],numpy.complex_)
 
             # Read the projectors from the file:
-            for ik in xrange(n_k):
+            for ik in range(n_k):
                 for icrsh in range(n_corr_shells):
                     for isp in range(n_spin_blocs):
 
                         # calculate the offset:
                         offset = 0
-                        no = 0
-                        for i in range(n_shells):
-                            if (no==0):
-                                if ((shells[i][0]==corr_shells[icrsh][0]) and (shells[i][1]==corr_shells[icrsh][1])):
-                                    no = corr_shells[icrsh][3]
+                        n_orb = 0
+                        for ish in range(n_shells):
+                            if (n_orb==0):
+                                if (shells[ish]['atom']==corr_shells[icrsh]['atom']) and (shells[ish]['sort']==corr_shells[icrsh]['sort']):
+                                    n_orb = corr_shells[icrsh]['dim']
                                 else:
-                                    offset += shells[i][3]
+                                    offset += shells[ish]['dim']
 
-                        proj_mat[ik,isp,icrsh,0:no,offset:offset+no] = numpy.identity(no)
+                        proj_mat[ik,isp,icrsh,0:n_orb,offset:offset+n_orb] = numpy.identity(n_orb)
                     
             # now define the arrays for weights and hopping ...
             bz_weights = numpy.ones([n_k],numpy.float_)/ float(n_k)  # w(k_index),  default normalisation 
@@ -142,7 +143,7 @@ class HkConverter(ConverterTools):
 
             if (weights_in_file):
                 # weights in the file
-                for ik in xrange(n_k) : bz_weights[ik] = R.next()
+                for ik in range(n_k) : bz_weights[ik] = R.next()
                 
             # if the sum over spins is in the weights, take it out again!!
             sm = sum(bz_weights)
@@ -150,36 +151,36 @@ class HkConverter(ConverterTools):
 
             # Grab the H
             for isp in range(n_spin_blocs):
-                for ik in xrange(n_k) :
-                    no = n_orbitals[ik,isp]
+                for ik in range(n_k) :
+                    n_orb = n_orbitals[ik,isp]
 
                     if (first_real_part_matrix): # first read all real components for given k, then read imaginary parts
                         
-                        for i in xrange(no):
+                        for i in range(n_orb):
                             if (only_upper_triangle):
                                 istart = i
                             else:
                                 istart = 0
-                            for j in xrange(istart,no):
+                            for j in range(istart,n_orb):
                                 hopping[ik,isp,i,j] = R.next()
                 
-                        for i in xrange(no):
+                        for i in range(n_orb):
                             if (only_upper_triangle):
                                 istart = i
                             else:
                                 istart = 0
-                            for j in xrange(istart,no):
+                            for j in range(istart,n_orb):
                                 hopping[ik,isp,i,j] += R.next() * 1j
                                 if ((only_upper_triangle)and(i!=j)): hopping[ik,isp,j,i] = hopping[ik,isp,i,j].conjugate()
                 
                     else: # read (real,im) tuple
                     
-                        for i in xrange(no):
+                        for i in range(n_orb):
                             if (only_upper_triangle):
                                 istart = i
                             else:
                                 istart = 0
-                            for j in xrange(istart,no):
+                            for j in range(istart,n_orb):
                                 hopping[ik,isp,i,j] = R.next()
                                 hopping[ik,isp,i,j] += R.next() * 1j
                             
diff --git a/python/converters/wien2k_converter.py b/python/converters/wien2k_converter.py
index 7eb852a6..6bd1d7bf 100644
--- a/python/converters/wien2k_converter.py
+++ b/python/converters/wien2k_converter.py
@@ -71,8 +71,8 @@ class Wien2kConverter(ConverterTools):
         # R is a generator : each R.Next() will return the next number in the file
         R = ConverterTools.read_fortran_file(self,self.dft_file,self.fortran_to_replace)
         try:
-            energy_unit = R.next()                         # read the energy convertion factor
-            n_k = int(R.next())                            # read the number of k points
+            energy_unit = R.next()                        # read the energy convertion factor
+            n_k = int(R.next())                           # read the number of k points
             k_dep_projection = 1                          
             SP = int(R.next())                            # flag for spin-polarised calculation
             SO = int(R.next())                            # flag for spin-orbit calculation
@@ -82,29 +82,32 @@ class Wien2kConverter(ConverterTools):
 
             # the information on the non-correlated shells is not important here, maybe skip:
             n_shells = int(R.next())                      # number of shells (e.g. Fe d, As p, O p) in the unit cell, 
-                                                               # corresponds to index R in formulas
-            shells = [ [ int(R.next()) for i in range(4) ] for icrsh in range(n_shells) ]    # reads iatom, sort, l, dim
+                                                          # corresponds to index R in formulas
+            # now read the information about the shells (atom, sort, l, dim):
+            shell_entries = ['atom', 'sort', 'l', 'dim']
+            shells = [ {name: int(val) for name, val in zip(shell_entries, R)} for ish in range(n_shells) ]
 
             n_corr_shells = int(R.next())                 # number of corr. shells (e.g. Fe d, Ce f) in the unit cell, 
                                                           # corresponds to index R in formulas
-            # now read the information about the shells:
-            corr_shells = [ [ int(R.next()) for i in range(6) ] for icrsh in range(n_corr_shells) ]    # reads iatom, sort, l, dim, SO flag, irep
+            # now read the information about the shells (atom, sort, l, dim, SO flag, irep):
+            corr_shell_entries = ['atom', 'sort', 'l', 'dim', 'SO', 'irep']
+            corr_shells = [ {name: int(val) for name, val in zip(corr_shell_entries, R)} for icrsh in range(n_corr_shells) ]
 
             # determine the number of inequivalent correlated shells and maps, needed for further reading
-            [n_inequiv_shells, corr_to_inequiv, inequiv_to_corr] = ConverterTools.det_shell_equivalence(self,corr_shells)
+            n_inequiv_shells, corr_to_inequiv, inequiv_to_corr = ConverterTools.det_shell_equivalence(self,corr_shells)
 
             use_rotations = 1
-            rot_mat = [numpy.identity(corr_shells[icrsh][3],numpy.complex_) for icrsh in xrange(n_corr_shells)]
+            rot_mat = [numpy.identity(corr_shells[icrsh]['dim'],numpy.complex_) for icrsh in range(n_corr_shells)]
            
             # read the matrices
             rot_mat_time_inv = [0 for i in range(n_corr_shells)]
 
-            for icrsh in xrange(n_corr_shells):
-                for i in xrange(corr_shells[icrsh][3]):    # read real part:
-                    for j in xrange(corr_shells[icrsh][3]):
+            for icrsh in range(n_corr_shells):
+                for i in range(corr_shells[icrsh]['dim']):    # read real part:
+                    for j in range(corr_shells[icrsh]['dim']):
                         rot_mat[icrsh][i,j] = R.next()
-                for i in xrange(corr_shells[icrsh][3]):    # read imaginary part:
-                    for j in xrange(corr_shells[icrsh][3]):
+                for i in range(corr_shells[icrsh]['dim']):    # read imaginary part:
+                    for j in range(corr_shells[icrsh]['dim']):
                         rot_mat[icrsh][i,j] += 1j * R.next()
 
                 if (SP==1):             # read time inversion flag:
@@ -114,23 +117,23 @@ class Wien2kConverter(ConverterTools):
             n_reps = [1 for i in range(n_inequiv_shells)]
             dim_reps = [0 for i in range(n_inequiv_shells)]
             T = []
-            for icrsh in range(n_inequiv_shells):
-                n_reps[icrsh] = int(R.next())   # number of representatives ("subsets"), e.g. t2g and eg
-                dim_reps[icrsh] = [int(R.next()) for i in range(n_reps[icrsh])]   # dimensions of the subsets
+            for ish in range(n_inequiv_shells):
+                n_reps[ish] = int(R.next())   # number of representatives ("subsets"), e.g. t2g and eg
+                dim_reps[ish] = [int(R.next()) for i in range(n_reps[ish])]   # dimensions of the subsets
             
                 # The transformation matrix:
                 # is of dimension 2l+1 without SO, and 2*(2l+1) with SO!
-                ll = 2*corr_shells[inequiv_to_corr[icrsh]][2]+1
-                lmax = ll * (corr_shells[inequiv_to_corr[icrsh]][4] + 1)
+                ll = 2*corr_shells[inequiv_to_corr[ish]]['l']+1
+                lmax = ll * (corr_shells[inequiv_to_corr[ish]]['SO'] + 1)
                 T.append(numpy.zeros([lmax,lmax],numpy.complex_))
                 
                 # now read it from file:
-                for i in xrange(lmax):
-                    for j in xrange(lmax):
-                        T[icrsh][i,j] = R.next()
-                for i in xrange(lmax):
-                    for j in xrange(lmax):
-                        T[icrsh][i,j] += 1j * R.next()
+                for i in range(lmax):
+                    for j in range(lmax):
+                        T[ish][i,j] = R.next()
+                for i in range(lmax):
+                    for j in range(lmax):
+                        T[ish][i,j] += 1j * R.next()
     
             # Spin blocks to be read:
             n_spin_blocs = SP + 1 - SO   
@@ -138,25 +141,25 @@ class Wien2kConverter(ConverterTools):
             # read the list of n_orbitals for all k points
             n_orbitals = numpy.zeros([n_k,n_spin_blocs],numpy.int)
             for isp in range(n_spin_blocs):
-                for ik in xrange(n_k):
+                for ik in range(n_k):
                     n_orbitals[ik,isp] = int(R.next())
             
             # Initialise the projectors:
-            proj_mat = numpy.zeros([n_k,n_spin_blocs,n_corr_shells,max(numpy.array(corr_shells)[:,3]),max(n_orbitals)],numpy.complex_)
+            proj_mat = numpy.zeros([n_k,n_spin_blocs,n_corr_shells,max([crsh['dim'] for crsh in corr_shells]),max(n_orbitals)],numpy.complex_)
 
             # Read the projectors from the file:
-            for ik in xrange(n_k):
+            for ik in range(n_k):
                 for icrsh in range(n_corr_shells):
-                    no = corr_shells[icrsh][3]
+                    n_orb = corr_shells[icrsh]['dim']
                     # first Real part for BOTH spins, due to conventions in dmftproj:
                     for isp in range(n_spin_blocs):
-                        for i in xrange(no):
-                            for j in xrange(n_orbitals[ik][isp]):
+                        for i in range(n_orb):
+                            for j in range(n_orbitals[ik][isp]):
                                 proj_mat[ik,isp,icrsh,i,j] = R.next()
                     # now Imag part:
                     for isp in range(n_spin_blocs):
-                        for i in xrange(no):
-                            for j in xrange(n_orbitals[ik][isp]):
+                        for i in range(n_orb):
+                            for j in range(n_orbitals[ik][isp]):
                                 proj_mat[ik,isp,icrsh,i,j] += 1j * R.next()
           
             # now define the arrays for weights and hopping ...
@@ -164,7 +167,7 @@ class Wien2kConverter(ConverterTools):
             hopping = numpy.zeros([n_k,n_spin_blocs,max(n_orbitals),max(n_orbitals)],numpy.complex_)
 
             # weights in the file
-            for ik in xrange(n_k) : bz_weights[ik] = R.next()         
+            for ik in range(n_k) : bz_weights[ik] = R.next()         
                 
             # if the sum over spins is in the weights, take it out again!!
             sm = sum(bz_weights)
@@ -173,9 +176,9 @@ class Wien2kConverter(ConverterTools):
             # Grab the H
             # we use now the convention of a DIAGONAL Hamiltonian -- convention for Wien2K.
             for isp in range(n_spin_blocs):
-                for ik in xrange(n_k) :
-                    no = n_orbitals[ik,isp]
-                    for i in xrange(no):
+                for ik in range(n_k) :
+                    n_orb = n_orbitals[ik,isp]
+                    for i in range(n_orb):
                         hopping[ik,isp,i,i] = R.next() * energy_unit
             
             # keep some things that we need for reading parproj:
@@ -211,7 +214,7 @@ class Wien2kConverter(ConverterTools):
         if not (mpi.is_master_node()): return
         mpi.report("Reading parproj input from %s..."%self.parproj_file)
 
-        dens_mat_below = [ [numpy.zeros([self.shells[ish][3],self.shells[ish][3]],numpy.complex_) for ish in range(self.n_shells)] 
+        dens_mat_below = [ [numpy.zeros([self.shells[ish]['dim'],self.shells[ish]['dim']],numpy.complex_) for ish in range(self.n_shells)] 
                            for isp in range(self.n_spin_blocs) ]
 
         R = ConverterTools.read_fortran_file(self,self.parproj_file,self.fortran_to_replace)
@@ -220,44 +223,44 @@ class Wien2kConverter(ConverterTools):
         n_parproj = numpy.array(n_parproj)
                 
         # Initialise P, here a double list of matrices:
-        proj_mat_pc = numpy.zeros([self.n_k,self.n_spin_blocs,self.n_shells,max(n_parproj),max(numpy.array(self.shells)[:,3]),max(self.n_orbitals)],numpy.complex_)
+        proj_mat_pc = numpy.zeros([self.n_k,self.n_spin_blocs,self.n_shells,max(n_parproj),max([sh['dim'] for sh in self.shells]),max(self.n_orbitals)],numpy.complex_)
         
-        rot_mat_all = [numpy.identity(self.shells[ish][3],numpy.complex_) for ish in xrange(self.n_shells)]
+        rot_mat_all = [numpy.identity(self.shells[ish]['dim'],numpy.complex_) for ish in range(self.n_shells)]
         rot_mat_all_time_inv = [0 for i in range(self.n_shells)]
 
         for ish in range(self.n_shells):
             # read first the projectors for this orbital:
-            for ik in xrange(self.n_k):
+            for ik in range(self.n_k):
                 for ir in range(n_parproj[ish]):
 
                     for isp in range(self.n_spin_blocs):
-                        for i in xrange(self.shells[ish][3]):    # read real part:
-                            for j in xrange(self.n_orbitals[ik][isp]):
+                        for i in range(self.shells[ish]['dim']):    # read real part:
+                            for j in range(self.n_orbitals[ik][isp]):
                                 proj_mat_pc[ik,isp,ish,ir,i,j] = R.next()
                             
                     for isp in range(self.n_spin_blocs):
-                        for i in xrange(self.shells[ish][3]):    # read imaginary part:
-                            for j in xrange(self.n_orbitals[ik][isp]):
+                        for i in range(self.shells[ish]['dim']):    # read imaginary part:
+                            for j in range(self.n_orbitals[ik][isp]):
                                 proj_mat_pc[ik,isp,ish,ir,i,j] += 1j * R.next()
                                         
                     
             # now read the Density Matrix for this orbital below the energy window:
             for isp in range(self.n_spin_blocs):
-                for i in xrange(self.shells[ish][3]):    # read real part:
-                    for j in xrange(self.shells[ish][3]):
+                for i in range(self.shells[ish]['dim']):    # read real part:
+                    for j in range(self.shells[ish]['dim']):
                         dens_mat_below[isp][ish][i,j] = R.next()
             for isp in range(self.n_spin_blocs):
-                for i in xrange(self.shells[ish][3]):    # read imaginary part:
-                    for j in xrange(self.shells[ish][3]):
+                for i in range(self.shells[ish]['dim']):    # read imaginary part:
+                    for j in range(self.shells[ish]['dim']):
                         dens_mat_below[isp][ish][i,j] += 1j * R.next()
                 if (self.SP==0): dens_mat_below[isp][ish] /= 2.0
 
             # Global -> local rotation matrix for this shell:
-            for i in xrange(self.shells[ish][3]):    # read real part:
-                for j in xrange(self.shells[ish][3]):
+            for i in range(self.shells[ish]['dim']):    # read real part:
+                for j in range(self.shells[ish]['dim']):
                     rot_mat_all[ish][i,j] = R.next()
-            for i in xrange(self.shells[ish][3]):    # read imaginary part:
-                for j in xrange(self.shells[ish][3]):
+            for i in range(self.shells[ish]['dim']):    # read imaginary part:
+                for j in range(self.shells[ish]['dim']):
                     rot_mat_all[ish][i,j] += 1j * R.next()
                     
             if (self.SP):
@@ -294,25 +297,25 @@ class Wien2kConverter(ConverterTools):
             # read the list of n_orbitals for all k points
             n_orbitals = numpy.zeros([n_k,self.n_spin_blocs],numpy.int)
             for isp in range(self.n_spin_blocs):
-                for ik in xrange(n_k):
+                for ik in range(n_k):
                     n_orbitals[ik,isp] = int(R.next())
 
             # Initialise the projectors:
-            proj_mat = numpy.zeros([n_k,self.n_spin_blocs,self.n_corr_shells,max(numpy.array(self.corr_shells)[:,3]),max(n_orbitals)],numpy.complex_)
+            proj_mat = numpy.zeros([n_k,self.n_spin_blocs,self.n_corr_shells,max([crsh['dim'] for crsh in corr_shells]),max(n_orbitals)],numpy.complex_)
 
             # Read the projectors from the file:
-            for ik in xrange(n_k):
+            for ik in range(n_k):
                 for icrsh in range(self.n_corr_shells):
-                    no = self.corr_shells[icrsh][3]
+                    n_orb = self.corr_shells[icrsh]['dim']
                     # first Real part for BOTH spins, due to conventions in dmftproj:
                     for isp in range(self.n_spin_blocs):
-                        for i in xrange(no):
-                            for j in xrange(n_orbitals[ik,isp]):
+                        for i in range(n_orb):
+                            for j in range(n_orbitals[ik,isp]):
                                 proj_mat[ik,isp,icrsh,i,j] = R.next()
                     # now Imag part:
                     for isp in range(self.n_spin_blocs):
-                        for i in xrange(no):
-                            for j in xrange(n_orbitals[ik,isp]):
+                        for i in range(n_orb):
+                            for j in range(n_orbitals[ik,isp]):
                                 proj_mat[ik,isp,icrsh,i,j] += 1j * R.next()
 
             hopping = numpy.zeros([n_k,self.n_spin_blocs,max(n_orbitals),max(n_orbitals)],numpy.complex_)
@@ -320,9 +323,9 @@ class Wien2kConverter(ConverterTools):
             # Grab the H
             # we use now the convention of a DIAGONAL Hamiltonian!!!!
             for isp in range(self.n_spin_blocs):
-                for ik in xrange(n_k) :
-                    no = n_orbitals[ik,isp]
-                    for i in xrange(no):
+                for ik in range(n_k) :
+                    n_orb = n_orbitals[ik,isp]
+                    for i in range(n_orb):
                         hopping[ik,isp,i,i] = R.next() * self.energy_unit
 
             # now read the partial projectors:
@@ -330,19 +333,19 @@ class Wien2kConverter(ConverterTools):
             n_parproj = numpy.array(n_parproj)
             
             # Initialise P, here a double list of matrices:
-            proj_mat_pc = numpy.zeros([n_k,self.n_spin_blocs,self.n_shells,max(n_parproj),max(numpy.array(self.shells)[:,3]),max(n_orbitals)],numpy.complex_)
+            proj_mat_pc = numpy.zeros([n_k,self.n_spin_blocs,self.n_shells,max(n_parproj),max([sh['dim'] for sh in self.shells]),max(n_orbitals)],numpy.complex_)
 
             for ish in range(self.n_shells):
-                for ik in xrange(n_k):
+                for ik in range(n_k):
                     for ir in range(n_parproj[ish]):
                         for isp in range(self.n_spin_blocs):
                                     
-                            for i in xrange(self.shells[ish][3]):    # read real part:
-                                for j in xrange(n_orbitals[ik,isp]):
+                            for i in range(self.shells[ish]['dim']):    # read real part:
+                                for j in range(n_orbitals[ik,isp]):
                                     proj_mat_pc[ik,isp,ish,ir,i,j] = R.next()
                             
-                            for i in xrange(self.shells[ish][3]):    # read imaginary part:
-                                for j in xrange(n_orbitals[ik,isp]):
+                            for i in range(self.shells[ish]['dim']):    # read imaginary part:
+                                for j in range(n_orbitals[ik,isp]):
                                     proj_mat_pc[ik,isp,ish,ir,i,j] += 1j * R.next()
 
         except StopIteration : # a more explicit error if the file is corrupted.
@@ -375,36 +378,36 @@ class Wien2kConverter(ConverterTools):
         try:
             n_symm = int(R.next())           # Number of symmetry operations
             n_atoms = int(R.next())       # number of atoms involved
-            perm = [ [int(R.next()) for i in xrange(n_atoms)] for j in xrange(n_symm) ]    # list of permutations of the atoms
+            perm = [ [int(R.next()) for i in range(n_atoms)] for j in range(n_symm) ]    # list of permutations of the atoms
             if SP: 
-                time_inv = [ int(R.next()) for j in xrange(n_symm) ]           # time inversion for SO coupling
+                time_inv = [ int(R.next()) for j in range(n_symm) ]           # time inversion for SO coupling
             else:
-                time_inv = [ 0 for j in xrange(n_symm) ]
+                time_inv = [ 0 for j in range(n_symm) ]
 
             # Now read matrices:
             mat = []  
-            for i_symm in xrange(n_symm):
+            for i_symm in range(n_symm):
                 
-                mat.append( [ numpy.zeros([orbits[orb][3], orbits[orb][3]],numpy.complex_) for orb in xrange(n_orbits) ] )
+                mat.append( [ numpy.zeros([orbits[orb]['dim'], orbits[orb]['dim']],numpy.complex_) for orb in range(n_orbits) ] )
                 for orb in range(n_orbits):
-                    for i in xrange(orbits[orb][3]):
-                        for j in xrange(orbits[orb][3]):
+                    for i in range(orbits[orb]['dim']):
+                        for j in range(orbits[orb]['dim']):
                             mat[i_symm][orb][i,j] = R.next()            # real part
-                    for i in xrange(orbits[orb][3]):
-                        for j in xrange(orbits[orb][3]):
+                    for i in range(orbits[orb]['dim']):
+                        for j in range(orbits[orb]['dim']):
                             mat[i_symm][orb][i,j] += 1j * R.next()      # imaginary part
 
-            mat_tinv = [numpy.identity(orbits[orb][3],numpy.complex_)
+            mat_tinv = [numpy.identity(orbits[orb]['dim'],numpy.complex_)
                         for orb in range(n_orbits)]
 
             if ((SO==0) and (SP==0)):
                 # here we need an additional time inversion operation, so read it:
                 for orb in range(n_orbits):
-                    for i in xrange(orbits[orb][3]):
-                        for j in xrange(orbits[orb][3]):
+                    for i in range(orbits[orb]['dim']):
+                        for j in range(orbits[orb]['dim']):
                             mat_tinv[orb][i,j] = R.next()            # real part
-                    for i in xrange(orbits[orb][3]):
-                        for j in xrange(orbits[orb][3]):
+                    for i in range(orbits[orb]['dim']):
+                        for j in range(orbits[orb]['dim']):
                             mat_tinv[orb][i,j] += 1j * R.next()      # imaginary part
                 
 
diff --git a/python/sumk_dft.py b/python/sumk_dft.py
index c7a7fc77..c7c8fd8b 100644
--- a/python/sumk_dft.py
+++ b/python/sumk_dft.py
@@ -73,7 +73,7 @@ class SumkDFT:
 
             # GF structure used for the local things in the k sums
             # Most general form allowing for all hybridisation, i.e. largest blocks possible
-            self.gf_struct_sumk = [ [ (sp, range( self.corr_shells[icrsh][3])) for sp in self.spin_block_names[self.corr_shells[icrsh][4]] ]
+            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) ]
 
             #-----
@@ -83,8 +83,8 @@ class SumkDFT:
                                                                               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]][3]) )
-                                                        for sp in self.spin_block_names[self.corr_shells[self.inequiv_to_corr[ish]][4]] ])
+                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
@@ -187,7 +187,7 @@ class SumkDFT:
 
         gf_downfolded = gf_inp.copy()
         isp = self.spin_names_to_ind[self.SO][bname]       # get spin index for proj. matrices
-        dim = self.corr_shells[icrsh][3]
+        dim = self.corr_shells[icrsh]['dim']
         n_orb = self.n_orbitals[ik,isp]
         projmat = self.proj_mat[ik,isp,icrsh,0:dim,0:n_orb]
 
@@ -201,7 +201,7 @@ class SumkDFT:
 
         gf_upfolded = gf_inp.copy()
         isp = self.spin_names_to_ind[self.SO][bname]       # get spin index for proj. matrices
-        dim = self.corr_shells[icrsh][3]
+        dim = self.corr_shells[icrsh]['dim']
         n_orb = self.n_orbitals[ik,isp]
         projmat = self.proj_mat[ik,isp,icrsh,0:dim,0:n_orb]
 
@@ -391,7 +391,7 @@ class SumkDFT:
         for ish in include_shells:
 
 
-            for sp in self.spin_block_names[self.corr_shells[self.inequiv_to_corr[ish]][4]]:
+            for sp in self.spin_block_names[self.corr_shells[self.inequiv_to_corr[ish]]['SO']]:
                 dmbool = (abs(dens_mat[ish][sp]) > threshold)  # gives an index list of entries larger that threshold
 
                 # Determine off-diagonal entries in upper triangular part of density matrix
@@ -474,8 +474,8 @@ class SumkDFT:
         """
         dens_mat = [ {} for icrsh in range(self.n_corr_shells)]
         for icrsh in range(self.n_corr_shells):
-            for sp in self.spin_block_names[self.corr_shells[icrsh][4]]:
-                dens_mat[icrsh][sp] = numpy.zeros([self.corr_shells[icrsh][3],self.corr_shells[icrsh][3]], numpy.complex_)
+            for sp in self.spin_block_names[self.corr_shells[icrsh]['SO']]:
+                dens_mat[icrsh][sp] = numpy.zeros([self.corr_shells[icrsh]['dim'],self.corr_shells[icrsh]['dim']], numpy.complex_)
 
         ikarray = numpy.array(range(self.n_k))
         for ik in mpi.slice_array(ikarray):
@@ -507,9 +507,9 @@ class SumkDFT:
                             MMat[isp][inu,inu] = 0.0
 
             for icrsh in range(self.n_corr_shells):
-                for isp, sp in enumerate(self.spin_block_names[self.corr_shells[icrsh][4]]):
-                    isp = self.spin_names_to_ind[self.corr_shells[icrsh][4]][sp]
-                    dim = self.corr_shells[icrsh][3]
+                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]
+                    dim = self.corr_shells[icrsh]['dim']
                     n_orb = self.n_orbitals[ik,isp]
                     projmat = self.proj_mat[ik,isp,icrsh,0:dim,0:n_orb]
                     if method == "using_gf":
@@ -545,8 +545,8 @@ class SumkDFT:
         # define matrices for inequivalent shells:
         eff_atlevels = [ {} for ish in range(self.n_inequiv_shells) ]
         for ish in range(self.n_inequiv_shells):
-            for sp in self.spin_block_names[self.corr_shells[self.inequiv_to_corr[ish]][4]]:
-                eff_atlevels[ish][sp] = numpy.identity(self.corr_shells[self.inequiv_to_corr[ish]][3], numpy.complex_)
+            for sp in self.spin_block_names[self.corr_shells[self.inequiv_to_corr[ish]]['SO']]:
+                eff_atlevels[ish][sp] = numpy.identity(self.corr_shells[self.inequiv_to_corr[ish]]['dim'], numpy.complex_)
 
         # Chemical Potential:
         for ish in range(self.n_inequiv_shells):
@@ -562,14 +562,14 @@ class SumkDFT:
             # calculate the sum over k. Does not depend on mu, so do it only once:
             self.Hsumk = [ {} for icrsh in range(self.n_corr_shells) ]
             for icrsh in range(self.n_corr_shells):
-                for sp in self.spin_block_names[self.corr_shells[icrsh][4]]:
-                    dim = self.corr_shells[icrsh][3]  #*(1+self.corr_shells[icrsh][4])
+                for sp in self.spin_block_names[self.corr_shells[icrsh]['SO']]:
+                    dim = self.corr_shells[icrsh]['dim']  #*(1+self.corr_shells[icrsh]['SO'])
                     self.Hsumk[icrsh][sp] = numpy.zeros([dim,dim],numpy.complex_)
 
             for icrsh in range(self.n_corr_shells):
-                dim = self.corr_shells[icrsh][3]
-                for isp, sp in enumerate(self.spin_block_names[self.corr_shells[icrsh][4]]):
-                    isp = self.spin_names_to_ind[self.corr_shells[icrsh][4]][sp]
+                dim = self.corr_shells[icrsh]['dim']
+                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]
                     for ik in range(self.n_k):
                         n_orb = self.n_orbitals[ik,isp]
                         MMat = numpy.identity(n_orb, numpy.complex_)
@@ -604,8 +604,8 @@ class SumkDFT:
         # construct the density matrix dm_imp and double counting arrays
         self.dc_imp = [ {} for icrsh in range(self.n_corr_shells)]
         for icrsh in range(self.n_corr_shells):
-            dim = self.corr_shells[icrsh][3]
-            spn = self.spin_block_names[self.corr_shells[icrsh][4]]
+            dim = self.corr_shells[icrsh]['dim']
+            spn = self.spin_block_names[self.corr_shells[icrsh]['SO']]
             for sp in spn: self.dc_imp[icrsh][sp] = numpy.zeros([dim,dim],numpy.float_)
         self.dc_energ = [0.0 for icrsh in range(self.n_corr_shells)]
 
@@ -624,9 +624,9 @@ class SumkDFT:
             if iorb != orb: continue # ignore this orbital
 
             Ncr = {}
-            dim = self.corr_shells[icrsh][3] #*(1+self.corr_shells[icrsh][4])
+            dim = self.corr_shells[icrsh]['dim'] #*(1+self.corr_shells[icrsh]['SO'])
 
-            spn = self.spin_block_names[self.corr_shells[icrsh][4]]
+            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
@@ -637,7 +637,7 @@ class SumkDFT:
 
             Ncrtot = 0.0
 
-            spn = self.spin_block_names[self.corr_shells[icrsh][4]]
+            spn = self.spin_block_names[self.corr_shells[icrsh]['SO']]
             for sp in spn: 
                 Ncrtot += Ncr[sp]
 
@@ -897,12 +897,12 @@ class SumkDFT:
     # Check that the density matrix from projectors (DM = P Pdagger) is correct (ie matches DFT)
     def check_projectors(self):
 
-        dens_mat = [numpy.zeros([self.corr_shells[icrsh][3],self.corr_shells[icrsh][3]],numpy.complex_)
+        dens_mat = [numpy.zeros([self.corr_shells[icrsh]['dim'],self.corr_shells[icrsh]['dim']],numpy.complex_)
                    for icrsh in range(self.n_corr_shells)]
 
         for ik in range(self.n_k):
             for icrsh in range(self.n_corr_shells):
-                dim = self.corr_shells[icrsh][3]
+                dim = self.corr_shells[icrsh]['dim']
                 n_orb = self.n_orbitals[ik,0]
                 projmat = self.proj_mat[ik,0,icrsh,0:dim,0:n_orb]
                 dens_mat[icrsh][:,:] += numpy.dot(projmat, projmat.transpose().conjugate()) * self.bz_weights[ik]
diff --git a/python/sumk_dft_tools.py b/python/sumk_dft_tools.py
index d9cfaafa..bb5dffb0 100644
--- a/python/sumk_dft_tools.py
+++ b/python/sumk_dft_tools.py
@@ -46,7 +46,7 @@ class SumkDFTTools(SumkDFT):
 
         gf_downfolded = gf_inp.copy()
         isp = self.spin_names_to_ind[self.SO][bname]       # get spin index for proj. matrices
-        dim = self.shells[ish][3]
+        dim = self.shells[ish]['dim']
         n_orb = self.n_orbitals[ik,isp]
         L=self.proj_mat_pc[ik,isp,ish,ir,0:dim,0:n_orb]
         R=self.proj_mat_pc[ik,isp,ish,ir,0:dim,0:n_orb].conjugate().transpose()
@@ -158,15 +158,15 @@ class SumkDFTTools(SumkDFT):
         DOSproj     = [ {} for ish in range(self.n_inequiv_shells) ]
         DOSproj_orb = [ {} for ish in range(self.n_inequiv_shells) ]
         for ish in range(self.n_inequiv_shells):
-            for bn in self.spin_block_names[self.corr_shells[self.inequiv_to_corr[ish]][4]]:
-                dim = self.corr_shells[self.inequiv_to_corr[ish]][3]
+            for bn in self.spin_block_names[self.corr_shells[self.inequiv_to_corr[ish]]['SO']]:
+                dim = self.corr_shells[self.inequiv_to_corr[ish]]['dim']
                 DOSproj[ish][bn] = numpy.zeros([n_om],numpy.float_)
                 DOSproj_orb[ish][bn] = numpy.zeros([n_om,dim,dim],numpy.float_)
 
         # init:
         Gloc = []
         for icrsh in range(self.n_corr_shells):
-            spn = self.spin_block_names[self.corr_shells[icrsh][4]]
+            spn = self.spin_block_names[self.corr_shells[icrsh]['SO']]
             glist = lambda : [ GfReFreq(indices = inner, window = (om_min,om_max), n_points = n_om) for block,inner in self.gf_struct_sumk[icrsh]]
             Gloc.append(BlockGf(name_list = spn, block_list = glist(),make_copies=False))
         for icrsh in range(self.n_corr_shells): Gloc[icrsh].zero()                        # initialize to zero
@@ -214,8 +214,8 @@ class SumkDFTTools(SumkDFT):
                     for i in range(n_om): f.write("%s    %s\n"%(om_mesh[i],DOSproj[ish][bn][i]))
                     f.close()
 
-                    for i in range(self.corr_shells[self.inequiv_to_corr[ish]][3]):
-                        for j in range(i,self.corr_shells[self.inequiv_to_corr[ish]][3]):
+                    for i in range(self.corr_shells[self.inequiv_to_corr[ish]]['dim']):
+                        for j in range(i,self.corr_shells[self.inequiv_to_corr[ish]]['dim']):
                             Fname = 'DOS'+bn+'_proj'+str(ish)+'_'+str(i)+'_'+str(j)+'.dat'
                             f=open(Fname,'w')
                             for iom in range(n_om): f.write("%s    %s\n"%(om_mesh[iom],DOSproj_orb[ish][bn][iom,i,j]))
@@ -246,7 +246,7 @@ class SumkDFTTools(SumkDFT):
 
         mu = self.chemical_potential
 
-        gf_struct_proj = [ [ (sp, range(self.shells[i][3])) for sp in self.spin_block_names[self.SO] ]  for i in range(self.n_shells) ]
+        gf_struct_proj = [ [ (sp, range(self.shells[i]['dim'])) for sp in self.spin_block_names[self.SO] ]  for i in range(self.n_shells) ]
         Gproj = [BlockGf(name_block_generator = [ (block,GfReFreq(indices = inner, mesh = self.Sigma_imp[0].mesh)) for block,inner in gf_struct_proj[ish] ], make_copies = False )
                  for ish in range(self.n_shells)]
         for ish in range(self.n_shells): Gproj[ish].zero()
@@ -262,7 +262,7 @@ class SumkDFTTools(SumkDFT):
         DOSproj_orb = [ {} for ish in range(self.n_shells) ]
         for ish in range(self.n_shells):
             for bn in self.spin_block_names[self.SO]:
-                dim = self.shells[ish][3]
+                dim = self.shells[ish]['dim']
                 DOSproj[ish][bn] = numpy.zeros([n_om],numpy.float_)
                 DOSproj_orb[ish][bn] = numpy.zeros([n_om,dim,dim],numpy.float_)
 
@@ -317,8 +317,8 @@ class SumkDFTTools(SumkDFT):
                     for i in range(n_om): f.write("%s    %s\n"%(Msh[i],DOSproj[ish][bn][i]))
                     f.close()
 
-                    for i in range(self.shells[ish][3]):
-                        for j in range(i,self.shells[ish][3]):
+                    for i in range(self.shells[ish]['dim']):
+                        for j in range(i,self.shells[ish]['dim']):
                             Fname = './DOScorr'+bn+'_proj'+str(ish)+'_'+str(i)+'_'+str(j)+'.dat'
                             f=open(Fname,'w')
                             for iom in range(n_om): f.write("%s    %s\n"%(Msh[iom],DOSproj_orb[ish][bn][iom,i,j]))
@@ -384,7 +384,7 @@ class SumkDFTTools(SumkDFT):
             for sp in spn: Akw[sp] = numpy.zeros([self.n_k, n_om ],numpy.float_)
         else:
             Akw = {}
-            for sp in spn: Akw[sp] = numpy.zeros([self.shells[ishell][3],self.n_k, n_om ],numpy.float_)
+            for sp in spn: Akw[sp] = numpy.zeros([self.shells[ishell]['dim'],self.n_k, n_om ],numpy.float_)
 
         if fermi_surface:
             om_minplot = -2.0*broadening
@@ -393,7 +393,7 @@ class SumkDFTTools(SumkDFT):
             for sp in spn: Akw[sp] = numpy.zeros([self.n_k,1],numpy.float_)
 
         if not ishell is None:
-            GFStruct_proj =  [ (sp, range(self.shells[ishell][3])) for sp in spn ]
+            GFStruct_proj =  [ (sp, range(self.shells[ishell]['dim'])) for sp in spn ]
             Gproj = BlockGf(name_block_generator = [ (block,GfReFreq(indices = inner, mesh = self.Sigma_imp[0].mesh)) for block,inner in GFStruct_proj ], make_copies = False)
             Gproj.zero()
 
@@ -427,7 +427,7 @@ class SumkDFTTools(SumkDFT):
 
                 for iom in range(n_om):
                     if (M[iom] > om_minplot) and (M[iom] < om_maxplot):
-                        for ish in range(self.shells[ishell][3]):
+                        for ish in range(self.shells[ishell]['dim']):
                             for ibn in spn:
                                 Akw[ibn][ish,ik,iom] = Gproj[ibn].data[iom,ish,ish].imag/(-3.1415926535)
 
@@ -471,7 +471,7 @@ class SumkDFTTools(SumkDFT):
 
             else:
                 for ibn in spn:
-                    for ish in range(self.shells[ishell][3]):
+                    for ish in range(self.shells[ishell]['dim']):
 
                         if invert_Akw:
                             maxAkw=Akw[ibn][ish,:,:].max()
@@ -505,11 +505,11 @@ class SumkDFTTools(SumkDFT):
         # Density matrix in the window
         spn = self.spin_block_names[self.SO]
         ntoi = self.spin_names_to_ind[self.SO]
-        self.dens_mat_window = [ [numpy.zeros([self.shells[ish][3],self.shells[ish][3]],numpy.complex_) for ish in range(self.n_shells)]
+        self.dens_mat_window = [ [numpy.zeros([self.shells[ish]['dim'],self.shells[ish]['dim']],numpy.complex_) for ish in range(self.n_shells)]
                                  for isp in range(len(spn)) ]    # init the density matrix
 
         mu = self.chemical_potential
-        GFStruct_proj = [ [ (sp, range(self.shells[i][3])) for sp in spn ]  for i in range(self.n_shells) ]
+        GFStruct_proj = [ [ (sp, range(self.shells[i]['dim'])) for sp in spn ]  for i in range(self.n_shells) ]
         if hasattr(self,"Sigma_imp"):
             Gproj = [BlockGf(name_block_generator = [ (block,GfImFreq(indices = inner, mesh = self.Sigma_imp[0].mesh)) for block,inner in GFStruct_proj[ish] ], make_copies = False)
                      for ish in range(self.n_shells)]
diff --git a/python/symmetry.py b/python/symmetry.py
index 0d3dec26..03842fc3 100644
--- a/python/symmetry.py
+++ b/python/symmetry.py
@@ -67,7 +67,7 @@ class Symmetry:
         for i_symm in range(self.n_symm):
             for iorb in range(self.n_orbits):
                 srch = copy.deepcopy(self.orbits[iorb])
-                srch[0] = self.perm[i_symm][self.orbits[iorb][0]-1]
+                srch['atom'] = self.perm[i_symm][self.orbits[iorb]['atom']-1]
                 self.orb_map[i_symm][iorb] = self.orbits.index(srch)
 
 
@@ -90,8 +90,8 @@ class Symmetry:
 
         for i_symm in range(self.n_symm):
             for iorb in range(self.n_orbits):
-                l = self.orbits[iorb][2]         # s, p, d, or f
-                dim = self.orbits[iorb][3]
+                l = self.orbits[iorb]['l']         # s, p, d, or f
+                dim = self.orbits[iorb]['dim']
                 jorb = self.orb_map[i_symm][iorb]
 
                 if isinstance(obj[0],BlockGf):
diff --git a/python/trans_basis.py b/python/trans_basis.py
index f81f4085..e483ea11 100644
--- a/python/trans_basis.py
+++ b/python/trans_basis.py
@@ -28,7 +28,7 @@ class TransBasis:
             self.SK = SK
 
         self.T = copy.deepcopy(self.SK.T[0])
-        self.w = numpy.identity(SK.corr_shells[0][3])
+        self.w = numpy.identity(SK.corr_shells[0]['dim'])
 
 
     def __call__(self, prop_to_be_diagonal = 'eal'):
@@ -88,7 +88,7 @@ class TransBasis:
 
         f = open(filename,'w')
         Tnew = self.T.conjugate()
-        dim = self.SK.corr_shells[0][3]
+        dim = self.SK.corr_shells[0]['dim']
 
         if self.SK.SO == 0:
 
diff --git a/python/update_archive.py b/python/update_archive.py
index 4c068aa5..936a3783 100644
--- a/python/update_archive.py
+++ b/python/update_archive.py
@@ -15,24 +15,37 @@ Please keep a copy of your old archive as this script is
 If you encounter any problem please report it on github!
 """
 
-def det_shell_equivalence(lst):
-    corr_to_inequiv = [0 for i in range(len(lst))]
+def convert_shells(shells):
+    shell_entries = ['atom', 'sort', 'l', 'dim']
+    return [ {name: int(val) for name, val in zip(shell_entries, shells[ish])} for ish in range(len(shells)) ]    
+
+def convert_corr_shells(corr_shells):
+    corr_shell_entries = ['atom', 'sort', 'l', 'dim', 'SO', 'irep']
+    return [ {name: int(val) for name, val in zip(corr_shell_entries, corr_shells[icrsh])} for icrsh in range(len(corr_shells)) ]    
+
+def det_shell_equivalence(corr_shells):
+    corr_to_inequiv = [0 for i in range(len(corr_shells))]
     inequiv_to_corr = [0]
     n_inequiv_shells = 1
-    tmp = [ lst[0][1:3] ]
-    if (len(lst)>1):
-        for i in range(len(lst)-1):
-            fnd = False
+
+    if len(corr_shells) > 1:
+        inequiv_sort = [ corr_shells[0]['sort'] ]
+        inequiv_l = [ corr_shells[0]['l'] ]
+        for i in range(len(corr_shells)-1):
+            is_equiv = False
             for j in range(n_inequiv_shells):
-                if (tmp[j]==lst[i+1][1:3]):
-                    fnd = True
+                if (inequiv_sort[j]==corr_shells[i+1]['sort']) and (inequiv_l[j]==corr_shells[i+1]['l']):
+                    is_equiv = True
                     corr_to_inequiv[i+1] = j
-            if (fnd==False):
+            if is_equiv==False:
                 corr_to_inequiv[i+1] = n_inequiv_shells
                 n_inequiv_shells += 1
-                tmp.append( lst[i+1][1:3] )
-                inequiv_to_corr.append(i+1)
-    return [n_inequiv_shells, corr_to_inequiv, inequiv_to_corr]
+                inequiv_sort.append( corr_shells[i+1]['sort'] )
+                inequiv_l.append( corr_shells[i+1]['l'] )
+                inequiv_to_corr.append( i+1 )
+
+    return n_inequiv_shells, corr_to_inequiv, inequiv_to_corr
+
 
 ### Main ###
 
@@ -60,13 +73,24 @@ for obj in move_to_output:
        A.copy('dft_input/'+obj,'dft_output/'+obj)
        del(A['dft_input'][obj])
 
+# Update shells and corr_shells to list of dicts
+shells_old = HDFArchive(filename,'r')['dft_input']['shells'] 
+corr_shells_old = HDFArchive(filename,'r')['dft_input']['corr_shells']
+shells = convert_shells(shells_old)
+corr_shells = convert_corr_shells(corr_shells_old)
+del(A['dft_input']['shells'])
+del(A['dft_input']['corr_shells'])
+A.close()
+# Need to use HDFArchive for the following
+HDFArchive(filename,'a')['dft_input']['shells'] = shells
+HDFArchive(filename,'a')['dft_input']['corr_shells'] = corr_shells
+A = h5py.File(filename)
+
 # Add shell equivalency quantities
-B = A['dft_input']
-corr_shells = HDFArchive(filename,'r')['dft_input']['corr_shells']
 equiv_shell_info = det_shell_equivalence(corr_shells)
-B['n_inequiv_shells'] = equiv_shell_info[0]
-B['corr_to_inequiv'] = equiv_shell_info[1]
-B['inequiv_to_corr'] = equiv_shell_info[2]
+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']
@@ -75,7 +99,13 @@ for group in groups:
     if group not in A.keys(): continue
     print "Changing n_s to n_symm ..."
     A[group].move('n_s','n_symm')
-
+    # Convert orbits to list of dicts
+    orbits_old = HDFArchive(filename,'r')[group]['orbits']
+    orbits = convert_corr_shells(orbits_old)
+    del(A[group]['orbits'])
+    A.close()
+    HDFArchive(filename,'a')[group]['orbits'] = orbits
+    A = h5py.File(filename)
 A.close()
 
 # Repack to reclaim disk space
diff --git a/test/SrVO3.h5 b/test/SrVO3.h5
index c141e14f..b8d81a0f 100644
Binary files a/test/SrVO3.h5 and b/test/SrVO3.h5 differ
diff --git a/test/hk_convert.output.h5 b/test/hk_convert.output.h5
index dd5bd4e3..6f288a7f 100644
Binary files a/test/hk_convert.output.h5 and b/test/hk_convert.output.h5 differ
diff --git a/test/srvo3_Gloc.py b/test/srvo3_Gloc.py
index a257d570..3f83b9e1 100644
--- a/test/srvo3_Gloc.py
+++ b/test/srvo3_Gloc.py
@@ -31,10 +31,10 @@ beta = 40
 # Init the SumK class
 SK=SumkDFT(hdf_file='SrVO3.h5',use_dft_blocks=True)
 
-num_orbitals = SK.corr_shells[0][3]
-l = SK.corr_shells[0][2]
-spin_names = ["up","down"]
-orb_names = ["%s"%i for i in range(num_orbitals)]
+num_orbitals = SK.corr_shells[0]['dim']
+l = SK.corr_shells[0]['l']
+spin_names = ['up','down']
+orb_names = ['%s'%i for i in range(num_orbitals)]
 orb_hybridized = False
 
 S = Solver(beta=beta, gf_struct=set_operator_structure(spin_names,orb_names,orb_hybridized))
diff --git a/test/wien2k_convert.output.h5 b/test/wien2k_convert.output.h5
index 7d24e52a..5f0fdbb4 100644
Binary files a/test/wien2k_convert.output.h5 and b/test/wien2k_convert.output.h5 differ