Tidy up of symmetry

*changed map -> orb_map in symmetry to avoid using python keyword

python/TODOFIX

@@ -8,6 +8,7 @@ Substitutions:
 * read_symmetry_input -> convert_symmetry_input
 * Symm_corr -> symmcorr
 * gf_struct_corr -> gf_struct_sumk
+* n_s -> n_symm
 
 internal substitutions:
 Symm_par --> symmpar

python/converters/wien2k_converter.py

@@ -374,26 +374,26 @@ class Wien2kConverter(ConverterTools):
         R = ConverterTools.read_fortran_file(self,symm_file,self.fortran_to_replace)
 
         try:
-            n_s = int(R.next())           # Number of symmetry operations
+            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_s) ]    # list of permutations of the atoms
+            perm = [ [int(R.next()) for i in xrange(n_atoms)] for j in xrange(n_symm) ]    # list of permutations of the atoms
             if SP: 
-                time_inv = [ int(R.next()) for j in xrange(n_s) ]           # timeinversion for SO xoupling
+                time_inv = [ int(R.next()) for j in xrange(n_symm) ]           # time inversion for SO coupling
             else:
-                time_inv = [ 0 for j in xrange(n_s) ] 
+                time_inv = [ 0 for j in xrange(n_symm) ]
 
             # Now read matrices:
             mat = []  
-            for in_s in xrange(n_s):
+            for i_symm in xrange(n_symm):
                 
                 mat.append( [ numpy.zeros([orbits[orb][3], orbits[orb][3]],numpy.complex_) for orb in xrange(n_orbits) ] )
                 for orb in range(n_orbits):
                     for i in xrange(orbits[orb][3]):
                         for j in xrange(orbits[orb][3]):
-                            mat[in_s][orb][i,j] = R.next()            # real part
+                            mat[i_symm][orb][i,j] = R.next()            # real part
                     for i in xrange(orbits[orb][3]):
                         for j in xrange(orbits[orb][3]):
-                            mat[in_s][orb][i,j] += 1j * R.next()      # imaginary part
+                            mat[i_symm][orb][i,j] += 1j * R.next()      # imaginary part
 
             mat_tinv = [numpy.identity(orbits[orb][3],numpy.complex_)
                         for orb in range(n_orbits)]
@@ -419,6 +419,6 @@ class Wien2kConverter(ConverterTools):
         # Save it to the HDF:
         ar=HDFArchive(self.hdf_file,'a')
         if not (symm_subgrp in ar): ar.create_group(symm_subgrp)
-        things_to_save = ['n_s','n_atoms','perm','orbits','SO','SP','time_inv','mat','mat_tinv']
+        things_to_save = ['n_symm','n_atoms','perm','orbits','SO','SP','time_inv','mat','mat_tinv']
         for it in things_to_save: ar[symm_subgrp][it] = locals()[it]
         del ar

python/symmetry.py

@@ -20,15 +20,12 @@
 #
 ################################################################################
 
-
 import copy,numpy
-import string
 from types import *
 from pytriqs.gf.local import *
 from pytriqs.archive import *
 import pytriqs.utility.mpi as mpi
 
-
 class Symmetry:
     """This class provides the routines for applying symmetry operations for the k sums.
        It contains the permutations of the atoms in the unti cell, and the corresponding
@@ -38,18 +35,19 @@ class Symmetry:
         """Initialises the class.
            Reads the permutations and rotation matrizes from the file, and constructs the mapping for
            the given orbitals. For each orbit a matrix is read!!!
-           SO: Flag for SO coupled calculations.
+           SO: Flag for spin-orbit coupling.
-           SP: Spin polarisation yes/no
+           SP: Flag for spin polarisation.
            """
 
-        assert type(hdf_file)==StringType,"hdf_file must be a filename"; self.hdf_file = hdf_file
+        assert type(hdf_file) == StringType, "hdf_file must be a filename"
-        things_to_read = ['n_s','n_atoms','perm','orbits','SO','SP','time_inv','mat','mat_tinv']
+        self.hdf_file = hdf_file
+        things_to_read = ['n_symm','n_atoms','perm','orbits','SO','SP','time_inv','mat','mat_tinv']
         for it in things_to_read: setattr(self,it,0)
 
-        if (mpi.is_master_node()):
+        if mpi.is_master_node():
             #Read the stuff on master:
             ar = HDFArchive(hdf_file,'a')
-            if (subgroup is None):
+            if subgroup is None:
                 ar2 = ar
             else:
                 ar2 = ar[subgroup]
@@ -62,81 +60,68 @@ class Symmetry:
         for it in things_to_read: setattr(self,it,mpi.bcast(getattr(self,it)))
 
         # now define the mapping of orbitals:
-        # self.map[iorb]=jorb gives the permutation of the orbitals as given in the list, when the
+        # self.orb_map[iorb] = jorb gives the permutation of the orbitals as given in the list, when the
         # permutation of the atoms is done:
         self.n_orbits = len(self.orbits)
-
+        self.orb_map = [ [0 for iorb in range(self.n_orbits)] for i_symm in range(self.n_symm) ]
-        self.map = [ [0 for iorb in range(self.n_orbits)] for in_s in range(self.n_s) ]
+        for i_symm in range(self.n_symm):
-        for in_s in range(self.n_s):
             for iorb in range(self.n_orbits):
-
                 srch = copy.deepcopy(self.orbits[iorb])
-                srch[0] = self.perm[in_s][self.orbits[iorb][0]-1]
+                srch[0] = self.perm[i_symm][self.orbits[iorb][0]-1]
-                self.map[in_s][iorb] = self.orbits.index(srch)
+                self.orb_map[i_symm][iorb] = self.orbits.index(srch)
-
 
 
     def symmetrize(self,obj):
 
-        assert isinstance(obj,list),"obj has to be a list of objects!"
+        assert isinstance(obj,list), "symmetry: obj has to be a list of objects."
-        assert len(obj)==self.n_orbits,"obj has to be a list of the same length as defined in the init"
+        assert len(obj) == self.n_orbits, "symmetry: obj has to be a list of the same length as defined in the init."
 
-        if (isinstance(obj[0],BlockGf)):
+        if isinstance(obj[0],BlockGf):
             symm_obj = [ obj[i].copy() for i in range(len(obj)) ]        # here the result is stored, it is a BlockGf!
             for iorb in range(self.n_orbits): symm_obj[iorb].zero()      # set to zero
         else:
             # if not a BlockGf, we assume it is a matrix (density matrix), has to be complex since self.mat is complex!
-            #symm_obj = [ numpy.zeros([self.orbits[iorb][3],self.orbits[iorb][3]],numpy.complex_) for iorb in range(self.n_orbits) ]
             symm_obj = [ copy.deepcopy(obj[i]) for i in range(len(obj)) ]
-
             for iorb in range(self.n_orbits):
-                if (type(symm_obj[iorb])==DictType):
+                if type(symm_obj[iorb]) == DictType:
                     for ii in symm_obj[iorb]: symm_obj[iorb][ii] *= 0.0
                 else:
                     symm_obj[iorb] *= 0.0
 
-
+        for i_symm in range(self.n_symm):
-        for in_s in range(self.n_s):
-
             for iorb in range(self.n_orbits):
-
                 l = self.orbits[iorb][2]         # s, p, d, or f
                 dim = self.orbits[iorb][3]
-                jorb = self.map[in_s][iorb]
+                jorb = self.orb_map[i_symm][iorb]
 
-
+                if isinstance(obj[0],BlockGf):
-                if (isinstance(obj[0],BlockGf)):
 
                     tmp = obj[iorb].copy()
-                    if (self.time_inv[in_s]): tmp << tmp.transpose()
+                    if self.time_inv[i_symm]: tmp << tmp.transpose()
-                    for bname,gf in tmp: tmp[bname].from_L_G_R(self.mat[in_s][iorb],tmp[bname],self.mat[in_s][iorb].conjugate().transpose())
+                    for bname,gf in tmp: tmp[bname].from_L_G_R(self.mat[i_symm][iorb],tmp[bname],self.mat[i_symm][iorb].conjugate().transpose())
-                    tmp *= 1.0/self.n_s
+                    tmp *= 1.0/self.n_symm
                     symm_obj[jorb] += tmp
 
                 else:
 
-                    if (type(obj[iorb])==DictType):
+                    if type(obj[iorb]) == DictType:
-
                         for ii in obj[iorb]:
-                            if (self.time_inv[in_s]==0):
+                            if self.time_inv[i_symm] == 0:
-                                symm_obj[jorb][ii] += numpy.dot(numpy.dot(self.mat[in_s][iorb],obj[iorb][ii]),
+                                symm_obj[jorb][ii] += numpy.dot(numpy.dot(self.mat[i_symm][iorb],obj[iorb][ii]),
-                                                                self.mat[in_s][iorb].conjugate().transpose()) / self.n_s
+                                                                self.mat[i_symm][iorb].conjugate().transpose()) / self.n_symm
                             else:
-                                symm_obj[jorb][ii] += numpy.dot(numpy.dot(self.mat[in_s][iorb],obj[iorb][ii].conjugate()),
+                                symm_obj[jorb][ii] += numpy.dot(numpy.dot(self.mat[i_symm][iorb],obj[iorb][ii].conjugate()),
-                                                                self.mat[in_s][iorb].conjugate().transpose()) / self.n_s
+                                                                self.mat[i_symm][iorb].conjugate().transpose()) / self.n_symm
-
-
-
                     else:
-                        if (self.time_inv[in_s]==0):
+                        if self.time_inv[i_symm] == 0:
-                            symm_obj[jorb] += numpy.dot(numpy.dot(self.mat[in_s][iorb],obj[iorb]),self.mat[in_s][iorb].conjugate().transpose()) / self.n_s
+                            symm_obj[jorb] += numpy.dot(numpy.dot(self.mat[i_symm][iorb],obj[iorb]),
+                                                        self.mat[i_symm][iorb].conjugate().transpose()) / self.n_symm
                         else:
-                            symm_obj[jorb] += numpy.dot(numpy.dot(self.mat[in_s][iorb],obj[iorb].conjugate()),
+                            symm_obj[jorb] += numpy.dot(numpy.dot(self.mat[i_symm][iorb],obj[iorb].conjugate()),
-                                                        self.mat[in_s][iorb].conjugate().transpose()) / self.n_s
+                                                        self.mat[i_symm][iorb].conjugate().transpose()) / self.n_symm
-
 
 # Markus: This does not what it is supposed to do, check how this should work (keep for now)
-#        if ((self.SO==0) and (self.SP==0)):
+#        if (self.SO == 0) and (self.SP == 0):
 #            # add time inv:
             #mpi.report("Add time inversion")
 #            for iorb in range(self.n_orbits):
@@ -148,7 +133,7 @@ class Symmetry:
 #                    symm_obj[iorb] /= 2.0
 #
 #                else:
-#                    if (type(symm_obj[iorb])==DictType):
+#                    if type(symm_obj[iorb]) == DictType:
 #                        for ii in symm_obj[iorb]:
 #                            symm_obj[iorb][ii] += numpy.dot(numpy.dot(self.mat_tinv[iorb],symm_obj[iorb][ii].conjugate()),
 #                                                            self.mat_tinv[iorb].transpose().conjugate())
@@ -158,9 +143,4 @@ class Symmetry:
 #                                                    self.mat_tinv[iorb].transpose().conjugate())
 #                        symm_obj[iorb] /= 2.0
 
-
         return symm_obj
-
-
-
-

python/update_archive.py

@@ -17,6 +17,7 @@ If you encounter any problem please report it on github!
 filename = sys.argv[1]
 A = h5py.File(filename)
 
+# Rename groups
 old_to_new = {'SumK_LDA':'lda_input', 'SumK_LDA_ParProj':'lda_parproj_input', 
  'SymmCorr':'lda_symmcorr_input', 'SymmPar':'lda_symmpar_input', 'SumK_LDA_Bands':'lda_bands_input'}
 
@@ -26,6 +27,7 @@ for old, new in old_to_new.iteritems():
     A.copy(old,new)
     del(A[old])
 
+# Move output items from lda_input to lda_output
 move_to_output = ['gf_struct_solver','map_inv','map',
                   'chemical_potential','dc_imp','dc_energ','deg_shells',
                   'h_field']
@@ -36,6 +38,14 @@ for obj in move_to_output:
        A.copy('lda_input/'+obj,'lda_output/'+obj)
        del(A['lda_input'][obj])
 
+# Rename variables
+groups = ['lda_symmcorr_input','lda_symmpar_input']
+
+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')
+
 A.close()
 
 # Repack to reclaim disk space