Tidying up

* Replace <<= with <<
* Reordering in put_Sigma
* Remove all instances of wien2triqs in doc

doc/Ce-gamma.py

@@ -55,7 +55,7 @@ if (previous_present):
     mpi.report("Using stored data for initialisation")
     if (mpi.is_master_node()):
         ar = HDFArchive(HDFfilename,'a')
-        S.Sigma <<= ar['SigmaImFreq']
+        S.Sigma << ar['SigmaImFreq']
         del ar
     S.Sigma = mpi.bcast(S.Sigma)
     SK.load()
@@ -75,7 +75,7 @@ for Iteration_Number in range(1,Loops+1):
             mpi.report("No adjustment of chemical potential\nTotal density  = %.3f"%SK.total_density(mu=Chemical_potential))
 
         # Density:
-        S.G <<= SK.extract_G_loc()[0]
+        S.G << SK.extract_G_loc()[0]
         mpi.report("Total charge of Gloc : %.6f"%S.G.total_density())
         dm = S.G.density()
 
@@ -104,9 +104,9 @@ for Iteration_Number in range(1,Loops+1):
             if (mpi.is_master_node()and (Mix<1.0)):
                 mpi.report("Mixing Sigma and G with factor %s"%Mix)
                 if ('SigmaImFreq' in ar):
-                    S.Sigma <<= Mix * S.Sigma + (1.0-Mix) * ar['SigmaImFreq']
+                    S.Sigma << Mix * S.Sigma + (1.0-Mix) * ar['SigmaImFreq']
                 if ('GF' in ar):
-                    S.G <<= Mix * S.G + (1.0-Mix) * ar['GF']
+                    S.G << Mix * S.G + (1.0-Mix) * ar['GF']
 
             S.G = mpi.bcast(S.G)
             S.Sigma = mpi.bcast(S.Sigma)

doc/LDADMFTmain.rst

@@ -93,8 +93,8 @@ set up the loop over DMFT iterations and the self-consistency condition::
 
           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.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.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,

doc/about.rst

@@ -1,8 +1,8 @@
 License
 =======
 
-Wien2TRIQS  is published under the `GNU General Public License, version 3
+The dft_tools package is published under the 
-<http://www.gnu.org/licenses/gpl.html>`_.
+`GNU General Public License, version 3 <http://www.gnu.org/licenses/gpl.html>`_.
 
 Authors & Quotation
 =======================

doc/advanced.rst

@@ -76,7 +76,7 @@ of the last iteration::
   if (previous_present):
     if (mpi.is_master_node()):
         ar = HDFArchive(lda_filename+'.h5','a')
-        S.Sigma <<= ar['SigmaImFreq']
+        S.Sigma << ar['SigmaImFreq']
         del ar
     S.Sigma = mpi.bcast(S.Sigma)
     
@@ -90,16 +90,16 @@ previous section, with some additional refinement::
         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
-        S.G <<= SK.extract_G_loc()[0]                           # calculation of the local Green function
+        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)
-            S.Sigma <<= SK.dc_imp[0]['up'][0,0]
+            S.Sigma << SK.dc_imp[0]['up'][0,0]
   
-        S.G0 <<= inverse(S.Sigma + inverse(S.G))
+        S.G0 << inverse(S.Sigma + inverse(S.G))
   
         # Solve the impurity problem:
         S.solve(U_interact=U,J_hund=J,use_spinflip=use_spinflip,use_matrix=use_matrix,
@@ -109,7 +109,7 @@ previous section, with some additional refinement::
         # solution done, do the post-processing:
         mpi.report("Total charge of impurity problem : %.6f"%S.G.total_density())
   
-        S.Sigma <<=(inverse(S.G0)-inverse(S.G))
+        S.Sigma <<(inverse(S.G0)-inverse(S.G))
         # Solve the impurity problem:
         S.solve(U_interact=U,J_hund=J,use_spinflip=use_spinflip,use_matrix=use_matrix,
                      l=l,T=SK.T[0], dim_reps=SK.dim_reps[0], irep=2, n_cycles=qmc_cycles,
@@ -118,15 +118,15 @@ previous section, with some additional refinement::
         # solution done, do the post-processing:
         mpi.report("Total charge of impurity problem : %.6f"%S.G.total_density())
   
-        S.Sigma <<=(inverse(S.G0)-inverse(S.G))
+        S.Sigma <<(inverse(S.G0)-inverse(S.G))
   
         # Now mix Sigma and G with factor Mix, if wanted:
         if ((iteration_number>1) or (previous_present)):
             if (mpi.is_master_node()):
                 ar = HDFArchive(lda_filename+'.h5','a')
                 mpi.report("Mixing Sigma and G with factor %s"%mix)
-                S.Sigma <<= mix * S.Sigma + (1.0-mix) * ar['Sigma']
+                S.Sigma << mix * S.Sigma + (1.0-mix) * ar['Sigma']
-                S.G <<= mix * S.G + (1.0-mix) * ar['GF']
+                S.G << mix * S.G + (1.0-mix) * ar['GF']
                 del ar
             S.G = mpi.bcast(S.G)
             S.Sigma = mpi.bcast(S.Sigma)

doc/conf.py.in

@@ -20,12 +20,12 @@ templates_path = ['@CMAKE_SOURCE_DIR@/doc/_templates']
 html_theme = 'triqs'
 html_theme_path = ['@TRIQS_THEMES_PATH@']
 html_show_sphinx = False
-html_context = {'header_title': 'Wien2TRIQS',
+html_context = {'header_title': 'dft_tools',
                 'header_subtitle': 'Connecting <a class="triqs" style="font-size: 12px" href="http://ipht.cea.fr/triqs">TRIQS</a>  to the Wien2k package',
                 'header_links': [['Install', 'install'],
                                  ['Documentation', 'documentation'],
                                  ['Issues', 'issues'],
-                                 ['About Wien2TRIQS', 'about']]}
+                                 ['About dft_tools', 'about']]}
 html_static_path = ['@CMAKE_SOURCE_DIR@/doc/_static']
 html_sidebars = {'index': ['sideb.html', 'searchbox.html']}
 

doc/index.rst

@@ -4,13 +4,13 @@
 
 .. _wien2k:
 
-Wien2TRIQS
+dft_tools
 ========================================================
 
 This application is aimed at DMFT calculations with 
 realistic band structure calculations. 
 A priori TRIQS can be connected to various realistic band structure codes.
-In this release, we provide the Wien2TRIQS extension module which contains an 
+In this release, we provide the dft_tools extension module which contains an 
 interface to the `Wien2k package <http://www.wien2k.at>`_.
 
 

doc/selfcons.rst

@@ -26,7 +26,7 @@ and store it in a format such that :program:`Wien2k` can read it. Therefore, aft
 previous section, we symmetrise the self energy, and recalculate the impurity Green function::
 
   SK.symm_deg_gf(S.Sigma,orb=0)
-  S.G <<= inverse(S.G0) - S.Sigma
+  S.G << inverse(S.G0) - S.Sigma
   S.G.invert()
 
 These steps are not necessary, but can help to reduce fluctuations in the total energy. 

lda_dmft_cthyb.py

@@ -24,7 +24,7 @@ use_blocks = True                # use bloc structure from LDA input
 prec_mu = 0.0001
 
 # Solver parameters
-p = SolverCore.solve_parameters()
+p = {}
 p["max_time"] = -1
 p["random_name"] = ""
 p["random_seed"] = 123 * mpi.rank + 567
@@ -58,16 +58,18 @@ spin_names = ["up","down"]
 orb_names = ["%s"%i for i in range(num_orbitals)]
 orb_hybridized = False
 
+# Construct U matrix for density-density calculations
+gf_struct = set_operator_structure(spin_names,orb_names,orb_hybridized)
 # Construct U matrix for density-density calculations
 Umat, Upmat = U_matrix_kanamori(n_orb=n_orb, U_int=U, J_hund=J)
 # Construct Hamiltonian and solver
-L = LocalProblem(spin_names, orb_names, orb_hybridized, h_loc_type="density", U=Umat, Uprime=Upmat, H_dump="H.txt")
+H = h_loc_density(spin_names, orb_names, orb_hybridized, U=Umat, Uprime=Upmat, H_dump="H.txt")
-S = Solver(beta=beta, gf_struct=L.gf_struct)
+S = Solver(beta=beta, gf_struct=gf_struct)
 
 if (previous_present):
   if (mpi.is_master_node()):
       ar = HDFArchive(lda_filename+'.h5','a')
-      S.Sigma_iw <<= ar['Sigma_iw']
+      S.Sigma_iw << ar['Sigma_iw']
       del ar
   S.Sigma_iw = mpi.bcast(S.Sigma_iw)
 
@@ -76,33 +78,33 @@ 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
-      S.G_iw <<= SK.extract_G_loc()[0]                           # extract the local Green function
+      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)
-          S.Sigma_iw <<= SK.dc_imp[0]['up'][0,0]
+          S.Sigma_iw << SK.dc_imp[0]['up'][0,0]
 
       # now calculate new G0_iw to input into the solver:
       if (mpi.is_master_node()):
           # We can do a mixing of Delta in order to stabilize the DMFT iterations:
-          S.G0_iw <<= S.Sigma_iw + inverse(S.G_iw)
+          S.G0_iw << S.Sigma_iw + inverse(S.G_iw)
           ar = HDFArchive(lda_filename+'.h5','a')
           if ((iteration_number>1) or (previous_present)):
               mpi.report("Mixing input Delta with factor %s"%delta_mix)
               Delta = (delta_mix * S.G0_iw.delta()) + (1.0-delta_mix) * ar['Delta_iw']
-              S.G0_iw <<= S.G0_iw + S.G0_iw.delta() - Delta
+              S.G0_iw << S.G0_iw + S.G0_iw.delta() - Delta
 
           ar['Delta_iw'] = S.G0_iw.delta()
-          S.G0_iw <<= inverse(S.G0_iw)
+          S.G0_iw << inverse(S.G0_iw)
           del ar
 
       S.G0_iw = mpi.bcast(S.G0_iw)
 
       # Solve the impurity problem:
-      S.solve(h_loc=L.h_loc, params=p)
+      S.solve(h_loc=h_loc, **p)
 
       # solution done, do the post-processing:
       mpi.report("Total charge of impurity problem : %.6f"%S.G_iw.total_density())
@@ -112,8 +114,8 @@ for iteration_number in range(1,loops+1):
           if (mpi.is_master_node()):
               ar = HDFArchive(lda_filename+'.h5','a')
               mpi.report("Mixing Sigma and G with factor %s"%sigma_mix)
-              S.Sigma_iw <<= sigma_mix * S.Sigma_iw + (1.0-sigma_mix) * ar['Sigma_iw']
+              S.Sigma_iw << sigma_mix * S.Sigma_iw + (1.0-sigma_mix) * ar['Sigma_iw']
-              S.G_iw <<= sigma_mix * S.G_iw + (1.0-sigma_mix) * ar['G_iw']
+              S.G_iw << sigma_mix * S.G_iw + (1.0-sigma_mix) * ar['G_iw']
               del ar
           S.G_iw = mpi.bcast(S.G_iw)
           S.Sigma_iw = mpi.bcast(S.Sigma_iw)

python/TODOFIX

@@ -0,0 +1,2 @@
+* realfreq --> w?
+* matsubara --> iw?

python/sumk_lda.py

@@ -29,9 +29,8 @@ from pytriqs.archive import *
 import pytriqs.utility.mpi as mpi
 from math import cos,sin
 
-# FIXME PS: These two aren't used it seems
+# FIXME PS: This isn't used it seems
-from pytriqs.operators.operators2 import *
+# from pytriqs.operators.operators2 import *
-import string, pickle
 
 class SumkLDA:
     """This class provides a general SumK method for combining ab-initio code and pytriqs."""
@@ -54,7 +53,7 @@ class SumkLDA:
             self.symm_corr_data = symm_corr_data
             self.block_names = [ ['up','down'], ['ud'] ]
             self.n_spin_blocks_gf = [2,1]
-            self.Gupf = None
+            self.G_upfold = None
             self.h_field = h_field
 
             # read input from HDF:
@@ -63,11 +62,11 @@ class SumkLDA:
                               'rot_mat_time_inv','n_reps','dim_reps','T','n_orbitals','proj_mat','bz_weights','hopping']
             optional_things = ['gf_struct_solver','map_inv','map','chemical_potential','dc_imp','dc_energ','deg_shells']
 
-            self.retval = self.read_input_from_hdf(subgrp=self.lda_data,things_to_read=things_to_read,optional_things=optional_things)
+            self.read_value = self.read_input_from_hdf(subgrp=self.lda_data,things_to_read=things_to_read,optional_things=optional_things)
 
             if (self.SO) and (abs(self.h_field)>0.000001):
                 self.h_field=0.0
-                mpi.report("For SO, the external magnetic field is not implemented, setting it to 0!!")
+                mpi.report("For SO, the external magnetic field is not implemented, setting it to 0!")
 
 
             # determine the number of inequivalent correlated shells (self.n_inequiv_corr_shells)
@@ -81,10 +80,11 @@ class SumkLDA:
                     self.names_to_ind[ibl][self.block_names[ibl][inm]] = inm * self.SP #(self.Nspinblocs-1)
 
             # 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_corr = [ [ (al, range( self.corr_shells[i][3])) for al in self.block_names[self.corr_shells[i][4]] ]
                                    for i in xrange(self.n_corr_shells) ]
 
-            if not (self.retval['gf_struct_solver']):
+            if not (self.read_value['gf_struct_solver']):
                 # No gf_struct was stored in HDF, so first set a standard one:
                 self.gf_struct_solver = [ dict([ (al, range(self.corr_shells[self.invshellmap[i]][3]) )
                                                         for al in self.block_names[self.corr_shells[self.invshellmap[i]][4]] ])
@@ -97,14 +97,14 @@ class SumkLDA:
                         self.map[i][al] = [al for j in range( self.corr_shells[self.invshellmap[i]][3] ) ]
                         self.map_inv[i][al] = al
 
-            if not (self.retval['dc_imp']):
+            if not (self.read_value['dc_imp']):
                 # init the double counting:
                 self.__init_dc()
 
-            if not (self.retval['chemical_potential']):
+            if not (self.read_value['chemical_potential']):
                 self.chemical_potential = mu
 
-            if not (self.retval['deg_shells']):
+            if not (self.read_value['deg_shells']):
                 self.deg_shells = [ [] for i in range(self.n_inequiv_corr_shells)]
 
             if self.symm_op:
@@ -132,7 +132,7 @@ class SumkLDA:
         Reads data from the HDF file
         """
 
-        retval = True
+        read_value = True
         # init variables on all nodes:
         for it in things_to_read: exec "self.%s = 0"%it
         for it in optional_things: exec "self.%s = 0"%it
@@ -146,20 +146,20 @@ class SumkLDA:
                         exec "self.%s = ar['%s']['%s']"%(it,subgrp,it)
                     else:
                         mpi.report("Loading %s failed!"%it)
-                        retval = False
+                        read_value = False
 
-                if ((retval) and (len(optional_things)>0)):
+                if ((read_value) and (len(optional_things)>0)):
                     # if necessary things worked, now read optional things:
-                    retval = {}
+                    read_value = {}
                     for it in optional_things:
                         if (it in ar[subgrp]):
                             exec "self.%s = ar['%s']['%s']"%(it,subgrp,it)
-                            retval['%s'%it] = True
+                            read_value['%s'%it] = True
                         else:
-                            retval['%s'%it] = False
+                            read_value['%s'%it] = False
             else:
                 mpi.report("Loading failed: No %s subgroup in HDF5!"%subgrp)
-                retval = False
+                read_value = False
 
             del ar
 
@@ -168,9 +168,9 @@ class SumkLDA:
         for it in optional_things: exec "self.%s = mpi.bcast(self.%s)"%(it,it)
 
 
-        retval = mpi.bcast(retval)
+        read_value = mpi.bcast(read_value)
 
-        return retval
+        return read_value
 
 
 
@@ -191,8 +191,8 @@ class SumkLDA:
 
         things_to_read=['chemical_potential','dc_imp','dc_energ']
 
-        retval = self.read_input_from_hdf(subgrp=self.lda_data,things_to_read=things_to_read)
+        read_value = self.read_input_from_hdf(subgrp=self.lda_data,things_to_read=things_to_read)
-        return retval
+        return read_value
 
 
     def downfold(self,ik,icrsh,sig,gf_to_downfold,gf_inp):
@@ -232,7 +232,7 @@ class SumkLDA:
         if (direction=='toGlobal'):
 
             if ((self.rot_mat_time_inv[icrsh]==1) and (self.SO)):
-                gf_rotated <<= gf_rotated.transpose()
+                gf_rotated << gf_rotated.transpose()
                 gf_rotated.from_L_G_R(self.rot_mat[icrsh].conjugate(),gf_rotated,self.rot_mat[icrsh].transpose())
             else:
                 gf_rotated.from_L_G_R(self.rot_mat[icrsh],gf_rotated,self.rot_mat[icrsh].conjugate().transpose())
@@ -240,7 +240,7 @@ class SumkLDA:
         elif (direction=='toLocal'):
 
             if ((self.rot_mat_time_inv[icrsh]==1)and(self.SO)):
-                gf_rotated <<= gf_rotated.transpose()
+                gf_rotated << gf_rotated.transpose()
                 gf_rotated.from_L_G_R(self.rot_mat[icrsh].transpose(),gf_rotated,self.rot_mat[icrsh].conjugate())
             else:
                 gf_rotated.from_L_G_R(self.rot_mat[icrsh].conjugate().transpose(),gf_rotated,self.rot_mat[icrsh])
@@ -262,8 +262,8 @@ class SumkLDA:
             beta = self.Sigma_imp[0].mesh.beta        #override beta if Sigma is present
 
 #
-#        if (self.Gupf is None):
+#        if (self.G_upfold is None):
-#            # first setting up of Gupf
+#            # first setting up of G_upfold
 #            BS = [ range(self.n_orbitals[ik,ntoi[ib]]) for ib in bln ]
 #            gf_struct = [ (bln[ib], BS[ib]) for ib in range(self.n_spin_blocks_gf[self.SO]) ]
 #            a_list = [a for a,al in gf_struct]
@@ -271,16 +271,16 @@ class SumkLDA:
 #                glist = lambda : [ GfImFreq(indices = al, mesh = self.Sigma_imp[0].mesh) for a,al in gf_struct]
 #            else:
 #                glist = lambda : [ GfImFreq(indices = al, beta = beta) for a,al in gf_struct]
-#            self.Gupf = BlockGf(name_list = a_list, block_list = glist(),make_copies=False)
+#            self.G_upfold = BlockGf(name_list = a_list, block_list = glist(),make_copies=False)
-#            self.Gupf.zero()
+#            self.G_upfold.zero()
-#            self.Gupf_id = self.Gupf.copy()
+#            self.G_upfold_id = self.G_upfold.copy()
-#            self.Gupf_id <<= iOmega_n
+#            self.G_upfold_id << iOmega_n
 #
-#        GFsize = [ gf.N1 for sig,gf in self.Gupf]
+#        GFsize = [ gf.N1 for sig,gf in self.G_upfold]
 #        unchangedsize = all( [ self.n_orbitals[ik,ntoi[bln[ib]]]==GFsize[ib]
 #                               for ib in range(self.n_spin_blocks_gf[self.SO]) ] )
 #
-#        if ((not unchangedsize)or(self.Gupf.mesh.beta!=beta)):
+#        if ((not unchangedsize)or(self.G_upfold.mesh.beta!=beta)):
 #            BS = [ range(self.n_orbitals[ik,ntoi[ib]]) for ib in bln ]
 #            gf_struct = [ (bln[ib], BS[ib]) for ib in range(self.n_spin_blocks_gf[self.SO]) ]
 #            a_list = [a for a,al in gf_struct]
@@ -288,25 +288,25 @@ class SumkLDA:
 #                glist = lambda : [ GfImFreq(indices = al, mesh = self.Sigma_imp[0].mesh) for a,al in gf_struct]
 #            else:
 #                glist = lambda : [ GfImFreq(indices = al, beta = beta) for a,al in gf_struct]
-#            self.Gupf = BlockGf(name_list = a_list, block_list = glist(),make_copies=False)
+#            self.G_upfold = BlockGf(name_list = a_list, block_list = glist(),make_copies=False)
-#            self.Gupf.zero()
+#            self.G_upfold.zero()
-#            self.Gupf_id = self.Gupf.copy()
+#            self.G_upfold_id = self.G_upfold.copy()
-#            self.Gupf_id <<= iOmega_n
+#            self.G_upfold_id << iOmega_n
 #
 ###
 # FIXME PS Remove commented out code above if this works
-        # Do we need to set up Gupf?
+        # Do we need to set up G_upfold?
-        set_up_Gupf = False
+        set_up_G_upfold = False
-        if self.Gupf == None:
+        if self.G_upfold == None:
-            set_up_Gupf = True
+            set_up_G_upfold = True
         else:
-            GFsize = [ gf.N1 for sig,gf in self.Gupf]
+            GFsize = [ gf.N1 for sig,gf in self.G_upfold]
             unchangedsize = all( [ self.n_orbitals[ik,ntoi[bln[ib]]]==GFsize[ib]
                                    for ib in range(self.n_spin_blocks_gf[self.SO]) ] )
-            if ((not unchangedsize)or(self.Gupf.mesh.beta!=beta)): set_up_Gupf = True
+            if ((not unchangedsize)or(self.G_upfold.mesh.beta!=beta)): set_up_G_upfold = True
 
-        # Set up Gupf
+        # Set up G_upfold
-        if set_up_Gupf:
+        if set_up_G_upfold:
             BS = [ range(self.n_orbitals[ik,ntoi[ib]]) for ib in bln ]
             gf_struct = [ (bln[ib], BS[ib]) for ib in range(self.n_spin_blocks_gf[self.SO]) ]
             a_list = [a for a,al in gf_struct]
@@ -314,55 +314,31 @@ class SumkLDA:
                 glist = lambda : [ GfImFreq(indices = al, mesh = self.Sigma_imp[0].mesh) for a,al in gf_struct]
             else:
                 glist = lambda : [ GfImFreq(indices = al, beta = beta) for a,al in gf_struct]
-            self.Gupf = BlockGf(name_list = a_list, block_list = glist(),make_copies=False)
+            self.G_upfold = BlockGf(name_list = a_list, block_list = glist(),make_copies=False)
-            self.Gupf.zero()
+            self.G_upfold.zero()
-            self.Gupf_id = self.Gupf.copy()
+            self.G_upfold_id = self.G_upfold.copy()
-            self.Gupf_id <<= iOmega_n
+            self.G_upfold_id << iOmega_n
 ###
 
         idmat = [numpy.identity(self.n_orbitals[ik,ntoi[bl]],numpy.complex_) for bl in bln]
 
-        self.Gupf <<= self.Gupf_id
+        self.G_upfold << self.G_upfold_id
         M = copy.deepcopy(idmat)
         for ibl in range(self.n_spin_blocks_gf[self.SO]):
             ind = ntoi[bln[ibl]]
             n_orb = self.n_orbitals[ik,ind]
             M[ibl] = self.hopping[ik,ind,0:n_orb,0:n_orb] - (idmat[ibl]*mu) - (idmat[ibl] * self.h_field * (1-2*ibl))
-        self.Gupf -= M
+        self.G_upfold -= M
 
         if (with_Sigma):
             for icrsh in xrange(self.n_corr_shells):
-                for sig,gf in self.Gupf: gf -= self.upfold(ik,icrsh,sig,stmp[icrsh][sig],gf)
+                for sig,gf in self.G_upfold: gf -= self.upfold(ik,icrsh,sig,stmp[icrsh][sig],gf)
 
-        self.Gupf.invert()
+        self.G_upfold.invert()
 
-	return self.Gupf
+	return self.G_upfold
 
 
-    def check_projectors(self):
-
-        dens_mat = [numpy.zeros([self.corr_shells[ish][3],self.corr_shells[ish][3]],numpy.complex_)
-                   for ish in range(self.n_corr_shells)]
-
-        for ik in range(self.n_k):
-
-            for ish in range(self.n_corr_shells):
-                dim = self.corr_shells[ish][3]
-                n_orb = self.n_orbitals[ik,0]
-                dens_mat[ish][:,:] += numpy.dot(self.proj_mat[ik,0,ish,0:dim,0:n_orb],self.proj_mat[ik,0,ish,0:dim,0:n_orb].transpose().conjugate()) * self.bz_weights[ik]
-
-        if (self.symm_op!=0): dens_mat = self.Symm_corr.symmetrize(dens_mat)
-
-        # Rotate to local coordinate system:
-        if (self.use_rotations):
-            for icrsh in xrange(self.n_corr_shells):
-                if (self.rot_mat_time_inv[icrsh]==1): dens_mat[icrsh] = dens_mat[icrsh].conjugate()
-                dens_mat[icrsh] = numpy.dot( numpy.dot(self.rot_mat[icrsh].conjugate().transpose(),dens_mat[icrsh]) ,
-                                            self.rot_mat[icrsh] )
-
-
-        return dens_mat
-
 
 
     def simple_point_dens_mat(self):
@@ -391,7 +367,7 @@ class SumkLDA:
             for ibl,bl in enumerate(bln):
                 ind = ntoi[bl]
                 for inu in range(self.n_orbitals[ik,ind]):
-                    if ( (self.hopping[ik,ind,inu,inu]-self.h_field*(1-2*ibl)) < 0.0):
+                    if ( (self.hopping[ik,ind,inu,inu]-self.h_field*(1-2*ibl)) < 0.0): # ONLY WORKS FOR DIAGONAL HOPPING MATRIX (TRUE IN WIEN2K)
                         MMat[ibl][inu,inu] = 1.0
                     else:
                         MMat[ibl][inu,inu] = 0.0
@@ -427,7 +403,7 @@ class SumkLDA:
 
         return dens_mat
 
-
+    # calculate upfolded gf, then density matrix -- no assumptions on structure (ie diagonal or not)
     def density_gf(self,beta):
         """Calculates the density without setting up Gloc. It is useful for Hubbard I, and very fast."""
 
@@ -440,9 +416,9 @@ class SumkLDA:
 
         for ik in mpi.slice_array(ikarray):
 
-            Gupf = self.lattice_gf_matsubara(ik=ik, beta=beta, mu=self.chemical_potential)
+            G_upfold = self.lattice_gf_matsubara(ik=ik, beta=beta, mu=self.chemical_potential)
-            Gupf *= self.bz_weights[ik]
+            G_upfold *= self.bz_weights[ik]
-            dm = Gupf.density()
+            dm = G_upfold.density()
             MMat = [dm[bl] for bl in self.block_names[self.SO]]
 
             for icrsh in range(self.n_corr_shells):
@@ -579,9 +555,9 @@ class SumkLDA:
             ss.zero()
             Ndeg = len(degsh)
             for bl in degsh: ss += gf_to_symm[bl] / (1.0*Ndeg)
-            for bl in degsh: gf_to_symm[bl] <<= ss
+            for bl in degsh: gf_to_symm[bl] << ss
-
 
+# for simply dft input, get crystal field splittings. 
     def eff_atomic_levels(self):
         """Calculates the effective atomic levels needed as input for the Hubbard I Solver."""
 
@@ -618,7 +594,7 @@ class SumkLDA:
                         n_orb = self.n_orbitals[ik,isp]
                         MMat = numpy.identity(n_orb, numpy.complex_)
                         MMat = self.hopping[ik,isp,0:n_orb,0:n_orb] - (1-2*ibn) * self.h_field * MMat
-                        self.Hsumk[icrsh][bn] += self.bz_weights[ik] * numpy.dot( numpy.dot(self.proj_mat[ik,isp,icrsh,0:dim,0:n_orb],MMat), #self.hopping[ik][isp]) ,
+                        self.Hsumk[icrsh][bn] += self.bz_weights[ik] * numpy.dot( numpy.dot(self.proj_mat[ik,isp,icrsh,0:dim,0:n_orb],MMat),
                                                                                   self.proj_mat[ik,isp,icrsh,0:dim,0:n_orb].conjugate().transpose() )
 
             # symmetrisation:
@@ -703,19 +679,24 @@ class SumkLDA:
                 if (use_val is None):
 
                     if (use_dc_formula==0): # FLL
+
                         self.dc_energ[icrsh] = U_interact / 2.0 * Ncrtot * (Ncrtot-1.0)
                         for bl in a_list:
                             Uav = U_interact*(Ncrtot-0.5) - J_hund*(Ncr[bl] - 0.5)
                             self.dc_imp[icrsh][bl] *= Uav
                             self.dc_energ[icrsh]  -= J_hund / 2.0 * (Ncr[bl]) * (Ncr[bl]-1.0)
                             mpi.report("DC for shell %(icrsh)i and block %(bl)s = %(Uav)f"%locals())
+
                     elif (use_dc_formula==1): # Held's formula, with U_interact the interorbital onsite interaction
+
                         self.dc_energ[icrsh] = (U_interact + (M-1)*(U_interact-2.0*J_hund) + (M-1)*(U_interact-3.0*J_hund))/(2*M-1) / 2.0 * Ncrtot * (Ncrtot-1.0)
                         for bl in a_list:
                             Uav =(U_interact + (M-1)*(U_interact-2.0*J_hund) + (M-1)*(U_interact-3.0*J_hund))/(2*M-1) * (Ncrtot-0.5)
                             self.dc_imp[icrsh][bl] *= Uav
                             mpi.report("DC for shell %(icrsh)i and block %(bl)s = %(Uav)f"%locals())
+
                     elif (use_dc_formula==2): # AMF
+
                         self.dc_energ[icrsh] = 0.5 * U_interact * Ncrtot * Ncrtot
                         for bl in a_list:
                             Uav = U_interact*(Ncrtot - Ncr[bl]/M) - J_hund * (Ncr[bl] - Ncr[bl]/M)
@@ -740,57 +721,23 @@ class SumkLDA:
 
 
 
-
-
-    def find_dc(self,orb,guess,dens_mat,dens_req=None,precision=0.01):
-        """Searches for DC in order to fulfill charge neutrality.
-           If dens_req is given, then DC is set such that the LOCAL charge of orbital
-           orb coincides with dens_req."""
-
-        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)
-            if (dens_req is None):
-                return self.total_density(mu=mu)
-            else:
-                return self.extract_G_loc()[orb].total_density()
-
-
-        if (dens_req is None):
-            Dens_rel = self.density_required - self.charge_below
-        else:
-            Dens_rel = dens_req
-
-        dcnew = dichotomy.dichotomy(function = F,
-                                    x_init = guess, y_value = Dens_rel,
-                                    precision_on_y = precision, delta_x=0.5,
-                                    max_loops = 100, x_name="Double-Counting", y_name= "Total Density",
-                                    verbosity = 3)[0]
-
-        return dcnew
-
-
-
-
-
     def put_Sigma(self, Sigma_imp):
         """Puts the impurity self energies for inequivalent atoms into the class, respects the multiplicity of the atoms."""
 
         assert isinstance(Sigma_imp,list), "Sigma_imp has to be a list of Sigmas for the correlated shells, even if it is of length 1!"
         assert len(Sigma_imp)==self.n_inequiv_corr_shells, "give exactly one Sigma for each inequivalent corr. shell!"
 
-
         # init self.Sigma_imp:
-        if Sigma_imp[0].note == 'ReFreq':
+        if all(type(g) == GfImFreq for name,g in Sigma_imp[0]):
-            # Real frequency Sigma:
-            self.Sigma_imp = [ BlockGf( name_block_generator = [ (a,GfReFreq(indices = al, mesh = Sigma_imp[0].mesh)) for a,al in self.gf_struct_corr[i] ],
-                                  make_copies = False) for i in xrange(self.n_corr_shells) ]
-            for i in xrange(self.n_corr_shells): self.Sigma_imp[i].note='ReFreq'
-        else:
             # Imaginary frequency Sigma:
             self.Sigma_imp = [ BlockGf( name_block_generator = [ (a,GfImFreq(indices = al, mesh = Sigma_imp[0].mesh)) for a,al in self.gf_struct_corr[i] ],
                                   make_copies = False) for i in xrange(self.n_corr_shells) ]
+        elif all(type(g) == GfReFreq for name,g in Sigma_imp[0]):
+            # Real frequency Sigma:
+            self.Sigma_imp = [ BlockGf( name_block_generator = [ (a,GfReFreq(indices = al, mesh = Sigma_imp[0].mesh)) for a,al in self.gf_struct_corr[i] ],
+                                  make_copies = False) for i in xrange(self.n_corr_shells) ]
+        else:
+            raise ValueError, "This type of Sigma is not handled."
 
         # transform the CTQMC blocks to the full matrix:
         for icrsh in xrange(self.n_corr_shells):
@@ -816,12 +763,12 @@ class SumkLDA:
                         ind2 = orblist[j]	
                         ind1_imp = map_ind[bl][ind1]
                         ind2_imp = map_ind[bl][ind2]
-                        self.Sigma_imp[icrsh][self.map_inv[s][bl]][ind1_imp,ind2_imp] <<= Sigma_imp[s][bl][ind1,ind2]
+                        self.Sigma_imp[icrsh][self.map_inv[s][bl]][ind1_imp,ind2_imp] << Sigma_imp[s][bl][ind1,ind2]
 
         # rotation from local to global coordinate system:
         if (self.use_rotations):
             for icrsh in xrange(self.n_corr_shells):
-                for sig,gf in self.Sigma_imp[icrsh]: self.Sigma_imp[icrsh][sig] <<= self.rotloc(icrsh,gf,direction='toGlobal')
+                for sig,gf in self.Sigma_imp[icrsh]: self.Sigma_imp[icrsh][sig] << self.rotloc(icrsh,gf,direction='toGlobal')
 
 
 
@@ -832,10 +779,11 @@ class SumkLDA:
         sres = [s.copy() for s in self.Sigma_imp]
         for icrsh in xrange(self.n_corr_shells):
             for bl,gf in sres[icrsh]:
+                # Transform dc_imp to global coordinate system
                 dccont = numpy.dot(self.rot_mat[icrsh],numpy.dot(self.dc_imp[icrsh][bl],self.rot_mat[icrsh].conjugate().transpose()))
                 sres[icrsh][bl] -= dccont
 
-        return sres
+        return sres # list of self energies corrected by DC
 
 
 
@@ -846,34 +794,6 @@ class SumkLDA:
 
 
 
-    def sorts_of_atoms(self,lst):
-        """
-        This routine should determine the number of sorts in the double list lst
-        """
-        sortlst = [ lst[i][1] for i in xrange(len(lst)) ]
-        sortlst.sort()
-        sorts = 1
-        for i in xrange(len(sortlst)-1):
-            if sortlst[i+1]>sortlst[i]: sorts += 1
-
-        return sorts
-
-
-
-    def number_of_atoms(self,lst):
-        """
-        This routine should determine the number of atoms in the double list lst
-        """
-        atomlst = [ lst[i][0] for i in xrange(len(lst)) ]
-        atomlst.sort()
-        atoms = 1
-        for i in xrange(len(atomlst)-1):
-            if atomlst[i+1]>atomlst[i]: atoms += 1
-
-        return atoms
-
-
-
     def inequiv_shells(self,lst):
         """
         The number of inequivalent shells is calculated from lst, and a mapping is given as
@@ -948,7 +868,7 @@ class SumkLDA:
         return self.chemical_potential
 
 
-
+#FIXME NOT USED!
     def find_mu_nonint(self, dens_req, orb = None, beta = 40, precision = 0.01):
 
         def F(mu):
@@ -998,12 +918,12 @@ class SumkLDA:
 
             for icrsh in xrange(self.n_corr_shells):
                 tmp = Gloc[icrsh].copy()                  # init temporary storage
-                for sig,gf in tmp: tmp[sig] <<= self.downfold(ik,icrsh,sig,S[sig],gf)
+                for sig,gf in tmp: tmp[sig] << self.downfold(ik,icrsh,sig,S[sig],gf)
                 Gloc[icrsh] += tmp
 
         #collect data from mpi:
         for icrsh in xrange(self.n_corr_shells):
-            Gloc[icrsh] <<= mpi.all_reduce(mpi.world,Gloc[icrsh],lambda x,y : x+y)
+            Gloc[icrsh] << mpi.all_reduce(mpi.world,Gloc[icrsh],lambda x,y : x+y)
         mpi.barrier()
 
 
@@ -1014,7 +934,7 @@ class SumkLDA:
         # Gloc is rotated to the local coordinate system:
         if (self.use_rotations):
             for icrsh in xrange(self.n_corr_shells):
-                for sig,gf in Gloc[icrsh]: Gloc[icrsh][sig] <<= self.rotloc(icrsh,gf,direction='toLocal')
+                for sig,gf in Gloc[icrsh]: Gloc[icrsh][sig] << self.rotloc(icrsh,gf,direction='toLocal')
 
         # transform to CTQMC blocks:
         Glocret = [ BlockGf( name_block_generator = [ (a,GfImFreq(indices = al, mesh = Gloc[0].mesh)) for a,al in self.gf_struct_solver[i].iteritems() ],
@@ -1041,7 +961,7 @@ class SumkLDA:
                         ind2 = orblist[j]	
                         ind1_imp = map_ind[bl][ind1]
                         ind2_imp = map_ind[bl][ind2]
-                        Glocret[ish][bl][ind1,ind2] <<= Gloc[self.invshellmap[ish]][self.map_inv[ish][bl]][ind1_imp,ind2_imp]
+                        Glocret[ish][bl][ind1,ind2] << Gloc[self.invshellmap[ish]][self.map_inv[ish][bl]][ind1_imp,ind2_imp]
 
 
         # return only the inequivalent shells:
@@ -1146,3 +1066,90 @@ class SumkLDA:
 
 
         return deltaN, dens
+
+
+
+################
+# FIXME LEAVE UNDOCUMENTED
+################
+
+    def find_dc(self,orb,guess,dens_mat,dens_req=None,precision=0.01):
+        """Searches for DC in order to fulfill charge neutrality.
+           If dens_req is given, then DC is set such that the LOCAL charge of orbital
+           orb coincides with dens_req."""
+
+        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)
+            if (dens_req is None):
+                return self.total_density(mu=mu)
+            else:
+                return self.extract_G_loc()[orb].total_density()
+
+
+        if (dens_req is None):
+            Dens_rel = self.density_required - self.charge_below
+        else:
+            Dens_rel = dens_req
+
+        dcnew = dichotomy.dichotomy(function = F,
+                                    x_init = guess, y_value = Dens_rel,
+                                    precision_on_y = precision, delta_x=0.5,
+                                    max_loops = 100, x_name="Double-Counting", y_name= "Total Density",
+                                    verbosity = 3)[0]
+
+        return dcnew
+
+
+# FIXME Check that dens matrix from projectors (DM=PPdagger) is correct (ie matches DFT)
+    def check_projectors(self):
+
+        dens_mat = [numpy.zeros([self.corr_shells[ish][3],self.corr_shells[ish][3]],numpy.complex_)
+                   for ish in range(self.n_corr_shells)]
+
+        for ik in range(self.n_k):
+
+            for ish in range(self.n_corr_shells):
+                dim = self.corr_shells[ish][3]
+                n_orb = self.n_orbitals[ik,0]
+                dens_mat[ish][:,:] += numpy.dot(self.proj_mat[ik,0,ish,0:dim,0:n_orb],self.proj_mat[ik,0,ish,0:dim,0:n_orb].transpose().conjugate()) * self.bz_weights[ik]
+
+        if (self.symm_op!=0): dens_mat = self.Symm_corr.symmetrize(dens_mat)
+
+        # Rotate to local coordinate system:
+        if (self.use_rotations):
+            for icrsh in xrange(self.n_corr_shells):
+                if (self.rot_mat_time_inv[icrsh]==1): dens_mat[icrsh] = dens_mat[icrsh].conjugate()
+                dens_mat[icrsh] = numpy.dot( numpy.dot(self.rot_mat[icrsh].conjugate().transpose(),dens_mat[icrsh]) ,
+                                            self.rot_mat[icrsh] )
+
+
+        return dens_mat
+
+
+# FIXME DETERMINE EQUIVALENCY OF SHELLS
+    def sorts_of_atoms(self,lst):
+        """
+        This routine should determine the number of sorts in the double list lst
+        """
+        sortlst = [ lst[i][1] for i in xrange(len(lst)) ]
+        sortlst.sort()
+        sorts = 1
+        for i in xrange(len(sortlst)-1):
+            if sortlst[i+1]>sortlst[i]: sorts += 1
+
+        return sorts
+
+
+    def number_of_atoms(self,lst):
+        """
+        This routine should determine the number of atoms in the double list lst
+        """
+        atomlst = [ lst[i][0] for i in xrange(len(lst)) ]
+        atomlst.sort()
+        atoms = 1
+        for i in xrange(len(atomlst)-1):
+            if atomlst[i+1]>atomlst[i]: atoms += 1
+
+        return atoms

python/sumk_lda_tools.py

@@ -47,7 +47,7 @@ class SumkLDATools(SumkLDA):
     def __init__(self, hdf_file, mu = 0.0, h_field = 0.0, use_lda_blocks = False, lda_data = 'SumK_LDA', symm_corr_data = 'SymmCorr',
                  par_proj_data = 'SumK_LDA_ParProj', symm_par_data = 'SymmPar', bands_data = 'SumK_LDA_Bands'):
 
-        self.Gupf_refreq = None
+        self.G_upfold_refreq = None
         SumkLDA.__init__(self,hdf_file=hdf_file,mu=mu,h_field=h_field,use_lda_blocks=use_lda_blocks,lda_data=lda_data,
                           symm_corr_data=symm_corr_data,par_proj_data=par_proj_data,symm_par_data=symm_par_data,
                           bands_data=bands_data)
@@ -77,14 +77,14 @@ class SumkLDATools(SumkLDA):
         gf_rotated = gf_to_rotate.copy()
         if (direction=='toGlobal'):
             if ((self.rot_mat_all_time_inv[ish]==1) and (self.SO)):
-                gf_rotated <<= gf_rotated.transpose()
+                gf_rotated << gf_rotated.transpose()
                 gf_rotated.from_L_G_R(self.rot_mat_all[ish].conjugate(),gf_rotated,self.rot_mat_all[ish].transpose())
             else:
                 gf_rotated.from_L_G_R(self.rot_mat_all[ish],gf_rotated,self.rot_mat_all[ish].conjugate().transpose())
 
         elif (direction=='toLocal'):
             if ((self.rot_mat_all_time_inv[ish]==1)and(self.SO)):
-                gf_rotated <<= gf_rotated.transpose()
+                gf_rotated << gf_rotated.transpose()
                 gf_rotated.from_L_G_R(self.rot_mat_all[ish].transpose(),gf_rotated,self.rot_mat_all[ish].conjugate())
             else:
                 gf_rotated.from_L_G_R(self.rot_mat_all[ish].conjugate().transpose(),gf_rotated,self.rot_mat_all[ish])
@@ -102,13 +102,13 @@ class SumkLDATools(SumkLDA):
 
         if (not hasattr(self,"Sigma_imp")): with_Sigma=False
         if (with_Sigma):
-            assert self.Sigma_imp[0].note == 'ReFreq', "Real frequency Sigma needed for lattice_gf_realfreq!"
+            assert all(type(g) == GfReFreq for name,g in self.Sigma_imp[0]), "Real frequency Sigma needed for lattice_gf_realfreq!"
             stmp = self.add_dc()
         else:
             assert (not (mesh is None)),"Without Sigma, give the mesh=(om_min,om_max,n_points) for lattice_gf_realfreq!"
 
-        if (self.Gupf_refreq is None):
+        if (self.G_upfold_refreq is None):
-            # first setting up of Gupf_refreq
+            # first setting up of G_upfold_refreq
             BS = [ range(self.n_orbitals[ik,ntoi[ib]]) for ib in bln ]
             gf_struct = [ (bln[ib], BS[ib]) for ib in range(self.n_spin_blocks_gf[self.SO]) ]
             a_list = [a for a,al in gf_struct]
@@ -116,10 +116,10 @@ class SumkLDATools(SumkLDA):
                 glist = lambda : [ GfReFreq(indices = al, mesh =self.Sigma_imp[0].mesh) for a,al in gf_struct]
             else:
                 glist = lambda : [ GfReFreq(indices = al, window=(mesh[0],mesh[1]),n_points=mesh[2]) for a,al in gf_struct]
-            self.Gupf_refreq = BlockGf(name_list = a_list, block_list = glist(),make_copies=False)
+            self.G_upfold_refreq = BlockGf(name_list = a_list, block_list = glist(),make_copies=False)
-            self.Gupf_refreq.zero()
+            self.G_upfold_refreq.zero()
 
-        GFsize = [ gf.N1 for sig,gf in self.Gupf_refreq]
+        GFsize = [ gf.N1 for sig,gf in self.G_upfold_refreq]
         unchangedsize = all( [ self.n_orbitals[ik,ntoi[bln[ib]]]==GFsize[ib]
                                for ib in range(self.n_spin_blocks_gf[self.SO]) ] )
 
@@ -131,28 +131,28 @@ class SumkLDATools(SumkLDA):
                 glist = lambda : [ GfReFreq(indices = al, mesh =self.Sigma_imp[0].mesh) for a,al in gf_struct]
             else:
                 glist = lambda : [ GfReFreq(indices = al, window=(mesh[0],mesh[1]),n_points=mesh[2]) for a,al in gf_struct]
-            self.Gupf_refreq = BlockGf(name_list = a_list, block_list = glist(),make_copies=False)
+            self.G_upfold_refreq = BlockGf(name_list = a_list, block_list = glist(),make_copies=False)
-            self.Gupf_refreq.zero()
+            self.G_upfold_refreq.zero()
 
         idmat = [numpy.identity(self.n_orbitals[ik,ntoi[bl]],numpy.complex_) for bl in bln]
 
-        self.Gupf_refreq <<= Omega + 1j*broadening
+        self.G_upfold_refreq << Omega + 1j*broadening
         M = copy.deepcopy(idmat)
         for ibl in range(self.n_spin_blocks_gf[self.SO]):
             ind = ntoi[bln[ibl]]
             n_orb = self.n_orbitals[ik,ind]
             M[ibl] = self.hopping[ik,ind,0:n_orb,0:n_orb] - (idmat[ibl]*mu) - (idmat[ibl] * self.h_field * (1-2*ibl))
-        self.Gupf_refreq -= M
+        self.G_upfold_refreq -= M
 
         if (with_Sigma):
-            tmp = self.Gupf_refreq.copy()    # init temporary storage
+            tmp = self.G_upfold_refreq.copy()    # init temporary storage
             for icrsh in xrange(self.n_corr_shells):
-                for sig,gf in tmp: tmp[sig] <<= self.upfold(ik,icrsh,sig,stmp[icrsh][sig],gf)
+                for sig,gf in tmp: tmp[sig] << self.upfold(ik,icrsh,sig,stmp[icrsh][sig],gf)
-                self.Gupf_refreq -= tmp      # adding to the upfolded GF
+                self.G_upfold_refreq -= tmp      # adding to the upfolded GF
 
-        self.Gupf_refreq.invert()
+        self.G_upfold_refreq.invert()
 
-        return self.Gupf_refreq
+        return self.G_upfold_refreq
 
 
 
@@ -186,18 +186,18 @@ class SumkLDATools(SumkLDA):
 
         for ik in xrange(self.n_k):
 
-            Gupf=self.lattice_gf_realfreq(ik=ik,mu=self.chemical_potential,broadening=broadening,mesh=(om_min,om_max,n_om),with_Sigma=False)
+            G_upfold=self.lattice_gf_realfreq(ik=ik,mu=self.chemical_potential,broadening=broadening,mesh=(om_min,om_max,n_om),with_Sigma=False)
-            Gupf *= self.bz_weights[ik]
+            G_upfold *= self.bz_weights[ik]
 
             # non-projected DOS
             for iom in range(n_om):
-                for sig,gf in Gupf:
+                for sig,gf in G_upfold:
                     asd = gf.data[iom,:,:].imag.trace()/(-3.1415926535)
                     DOS[sig][iom] += asd
 
             for icrsh in xrange(self.n_corr_shells):
                 tmp = Gloc[icrsh].copy()
-                for sig,gf in tmp: tmp[sig] <<= self.downfold(ik,icrsh,sig,Gupf[sig],gf) # downfolding G
+                for sig,gf in tmp: tmp[sig] << self.downfold(ik,icrsh,sig,G_upfold[sig],gf) # downfolding G
                 Gloc[icrsh] += tmp
 
 
@@ -206,7 +206,7 @@ class SumkLDATools(SumkLDA):
 
         if (self.use_rotations):
             for icrsh in xrange(self.n_corr_shells):
-                for sig,gf in Gloc[icrsh]: Gloc[icrsh][sig] <<= self.rotloc(icrsh,gf,direction='toLocal')
+                for sig,gf in Gloc[icrsh]: Gloc[icrsh][sig] << self.rotloc(icrsh,gf,direction='toLocal')
 
         # Gloc can now also be used to look at orbitally resolved quantities
         for ish in range(self.n_inequiv_corr_shells):
@@ -243,8 +243,8 @@ class SumkLDATools(SumkLDA):
         """
 
         thingstoread = ['dens_mat_below','n_parproj','proj_mat_pc','rot_mat_all','rot_mat_all_time_inv']
-        retval = self.read_input_from_hdf(subgrp=self.par_proj_data,things_to_read = thingstoread)
+        read_value = self.read_input_from_hdf(subgrp=self.par_proj_data,things_to_read = thingstoread)
-        return retval
+        return read_value
 
 
 
@@ -254,9 +254,9 @@ class SumkLDATools(SumkLDA):
         assert hasattr(self,"Sigma_imp"), "Set Sigma First!!"
 
         #thingstoread = ['Dens_Mat_below','N_parproj','Proj_Mat_pc','rotmat_all']
-        #retval = self.read_input_from_HDF(SubGrp=self.par_proj_data, thingstoread=thingstoread)
+        #read_value = self.read_input_from_HDF(SubGrp=self.par_proj_data, thingstoread=thingstoread)
-        retval = self.read_par_proj_input_from_hdf()
+        read_value = self.read_par_proj_input_from_hdf()
-        if not retval: return retval
+        if not read_value: return read_value
         if self.symm_op: self.Symm_par = Symmetry(self.hdf_file,subgroup=self.symm_par_data)
 
         mu = self.chemical_potential
@@ -296,14 +296,14 @@ class SumkLDATools(SumkLDA):
             for ish in xrange(self.n_shells):
                 tmp = Gproj[ish].copy()
                 for ir in xrange(self.n_parproj[ish]):
-                    for sig,gf in tmp: tmp[sig] <<= self.downfold_pc(ik,ir,ish,sig,S[sig],gf)
+                    for sig,gf in tmp: tmp[sig] << self.downfold_pc(ik,ir,ish,sig,S[sig],gf)
                     Gproj[ish] += tmp
 
         # collect data from mpi:
         for sig in DOS:
             DOS[sig] = mpi.all_reduce(mpi.world,DOS[sig],lambda x,y : x+y)
         for ish in xrange(self.n_shells):
-            Gproj[ish] <<= mpi.all_reduce(mpi.world,Gproj[ish],lambda x,y : x+y)
+            Gproj[ish] << mpi.all_reduce(mpi.world,Gproj[ish],lambda x,y : x+y)
         mpi.barrier()
 
         if (self.symm_op!=0): Gproj = self.Symm_par.symmetrize(Gproj)
@@ -311,7 +311,7 @@ class SumkLDATools(SumkLDA):
         # rotation to local coord. system:
         if (self.use_rotations):
             for ish in xrange(self.n_shells):
-                for sig,gf in Gproj[ish]: Gproj[ish][sig] <<= self.rotloc_all(ish,gf,direction='toLocal')
+                for sig,gf in Gproj[ish]: Gproj[ish][sig] << self.rotloc_all(ish,gf,direction='toLocal')
 
         for ish in range(self.n_shells):
             for sig,gf in Gproj[ish]:
@@ -348,8 +348,8 @@ class SumkLDATools(SumkLDA):
 
         assert hasattr(self,"Sigma_imp"), "Set Sigma First!!"
         thingstoread = ['n_k','n_orbitals','proj_mat','hopping','n_parproj','proj_mat_pc']
-        retval = self.read_input_from_hdf(subgrp=self.bands_data,things_to_read=thingstoread)
+        read_value = self.read_input_from_hdf(subgrp=self.bands_data,things_to_read=thingstoread)
-        if not retval: return retval
+        if not read_value: return read_value
 
         if fermi_surface: ishell=None
 
@@ -430,13 +430,13 @@ class SumkLDATools(SumkLDA):
                 Gproj.zero()
                 tmp = Gproj.copy()
                 for ir in xrange(self.n_parproj[ishell]):
-                    for sig,gf in tmp: tmp[sig] <<= self.downfold_pc(ik,ir,ishell,sig,S[sig],gf)
+                    for sig,gf in tmp: tmp[sig] << self.downfold_pc(ik,ir,ishell,sig,S[sig],gf)
                     Gproj += tmp
 
                 # TO BE FIXED:
                 # rotate to local frame
                 #if (self.use_rotations):
-                #    for sig,gf in Gproj: Gproj[sig] <<= self.rotloc(0,gf,direction='toLocal')
+                #    for sig,gf in Gproj: Gproj[sig] << self.rotloc(0,gf,direction='toLocal')
 
                 for iom in range(n_om):
                     if (M[iom]>om_minplot) and (M[iom]<om_maxplot):
@@ -601,9 +601,9 @@ class SumkLDATools(SumkLDA):
 
 
         #thingstoread = ['Dens_Mat_below','N_parproj','Proj_Mat_pc','rotmat_all']
-        #retval = self.read_input_from_HDF(SubGrp=self.par_proj_data,thingstoread=thingstoread)
+        #read_value = self.read_input_from_HDF(SubGrp=self.par_proj_data,thingstoread=thingstoread)
-        retval = self.read_par_proj_input_from_hdf()
+        read_value = self.read_par_proj_input_from_hdf()
-        if not retval: return retval
+        if not read_value: return read_value
         if self.symm_op: self.Symm_par = Symmetry(self.hdf_file,subgroup=self.symm_par_data)
 
         # Density matrix in the window
@@ -636,13 +636,13 @@ class SumkLDATools(SumkLDA):
             for ish in xrange(self.n_shells):
                 tmp = Gproj[ish].copy()
                 for ir in xrange(self.n_parproj[ish]):
-                    for sig,gf in tmp: tmp[sig] <<= self.downfold_pc(ik,ir,ish,sig,S[sig],gf)
+                    for sig,gf in tmp: tmp[sig] << self.downfold_pc(ik,ir,ish,sig,S[sig],gf)
                     Gproj[ish] += tmp
 
         #print "K-Sum done on node",mpi.rank," at ",datetime.now()
         #collect data from mpi:
         for ish in xrange(self.n_shells):
-            Gproj[ish] <<= mpi.all_reduce(mpi.world,Gproj[ish],lambda x,y : x+y)
+            Gproj[ish] << mpi.all_reduce(mpi.world,Gproj[ish],lambda x,y : x+y)
         mpi.barrier()
 
         #print "Data collected on node",mpi.rank," at ",datetime.now()
@@ -655,7 +655,7 @@ class SumkLDATools(SumkLDA):
 
             # Rotation to local:
             if (self.use_rotations):
-                for sig,gf in Gproj[ish]: Gproj[ish][sig] <<= self.rotloc_all(ish,gf,direction='toLocal')
+                for sig,gf in Gproj[ish]: Gproj[ish][sig] << self.rotloc_all(ish,gf,direction='toLocal')
 
             isp = 0
             for sig,gf in Gproj[ish]: #dmg.append(Gproj[ish].density()[sig])

python/symmetry.py

@@ -108,7 +108,7 @@ class Symmetry:
                 if (isinstance(obj[0],BlockGf)):
 
                     tmp = obj[iorb].copy()
-                    if (self.time_inv[in_s]): tmp <<= tmp.transpose()
+                    if (self.time_inv[in_s]): tmp << tmp.transpose()
                     for sig,gf in tmp: tmp[sig].from_L_G_R(self.mat[in_s][iorb],tmp[sig],self.mat[in_s][iorb].conjugate().transpose())
                     tmp *= 1.0/self.n_s
                     symm_obj[jorb] += tmp
@@ -142,7 +142,7 @@ class Symmetry:
 #            for iorb in range(self.n_orbits):
 #                if (isinstance(symm_obj[0],BlockGf)):
 #                    tmp = symm_obj[iorb].copy()
-#                    tmp <<= tmp.transpose()
+#                    tmp << tmp.transpose()
 #                    for sig,gf in tmp: tmp[sig].from_L_G_R(self.mat_tinv[iorb],tmp[sig],self.mat_tinv[iorb].transpose().conjugate())
 #                    symm_obj[iorb] += tmp
 #                    symm_obj[iorb] /= 2.0

python/trans_basis.py

@@ -83,7 +83,7 @@ class TransBasis:
                 for j in range(len(orblist)):
                     ind1 = orblist[i]
                     ind2 = orblist[j]
-                    gfrotated[self.SK.map_inv[s][bl]][ind1,ind2] <<= gf_to_rot[bl][ind1,ind2]
+                    gfrotated[self.SK.map_inv[s][bl]][ind1,ind2] << gf_to_rot[bl][ind1,ind2]
 
         # Rotate using the matrix w
         for sig,bn in gfrotated:
@@ -96,7 +96,7 @@ class TransBasis:
                 for j in range(len(orblist)):
                     ind1 = orblist[i]
                     ind2 = orblist[j]
-                    gfreturn[bl][ind1,ind2] <<= gfrotated[self.SK.map_inv[0][bl]][ind1,ind2]
+                    gfreturn[bl][ind1,ind2] << gfrotated[self.SK.map_inv[0][bl]][ind1,ind2]
 
         return gfreturn