From 259fd648244d2481a09b17276ebc201519511877 Mon Sep 17 00:00:00 2001 From: Priyanka Seth Date: Sun, 7 Dec 2014 00:29:39 +0100 Subject: [PATCH] Merged lattice_gf_matsubara and lattice_gf_realfreq into single function --- python/lattice_gf.py | 67 ---------- python/sumk_dft.py | 276 +++++++++++++++++++-------------------- python/sumk_dft_tools.py | 212 +++++++++++------------------- python/symmetry.py | 6 +- 4 files changed, 214 insertions(+), 347 deletions(-) delete mode 100644 python/lattice_gf.py diff --git a/python/lattice_gf.py b/python/lattice_gf.py deleted file mode 100644 index 9060fd7f..00000000 --- a/python/lattice_gf.py +++ /dev/null @@ -1,67 +0,0 @@ - def lattice_gf(self, ik, mu, iw_or_w="iw", beta=40, broadening, mesh=None, with_Sigma=True, with_dc=True): - """Calculates the lattice Green function from the DFT hopping and the self energy at k-point number ik - and chemical potential mu.""" - - ntoi = self.spin_names_to_ind[self.SO] - spn = self.spin_block_names[self.SO] - if (iw_or_w != "iw") and (iw_or_w != "w"): - raise ValueError, "lattice_gf: implemented only for Re/Im frequency functions." - - # Are we including Sigma? - if with_Sigma: - if with_dc: sigma_inc_dc = self.add_dc() - Sigma_imp = getattr(self,"Sigma_imp_"+iw_or_w) - if iw_or_w == "iw": beta = Sigma_imp[0].mesh.beta # override beta if Sigma_iw is present - else: - if (iw_or_w == "w") and (mesh is None): - raise ValueError, "Give the mesh=(om_min,om_max,n_points) for the lattice GfReFreq." - - # Check if G_upfold is present - set_up_G_upfold = False # Assume not - if not hasattr(self,"G_upfold_"+iw_or_w)): - set_up_G_upfold = True # Need to create G_upfold_iw - setattr(self,"G_upfold_"+iw_or_w,0) # Set temporarily to zero - else: # Check that existing GF is consistent - GFsize = [ gf.N1 for bname,gf in getattr(self,"G_upfold_"+iw_or_w)] - unchangedsize = all( [ self.n_orbitals[ik,ntoi[spn[isp]]]==GFsize[isp] - for isp in range(self.n_spin_blocks[self.SO]) ] ) - if not unchangedsize: set_up_G_upfold = True - if (iw_or_w != "iw") and (self.G_upfold_iw.mesh.beta != beta): # additional check for ImFreq - set_up_G_upfold = True - G_upfold = getattr(self,"G_upfold_"+iw_or_w) - - # Set up G_upfold - if set_up_G_upfold: - block_structure = [ range(self.n_orbitals[ik,ntoi[sp]]) for sp in spn ] - gf_struct = [ (spn[isp], block_structure[isp]) for isp in range(self.n_spin_blocks[self.SO]) ] - block_ind_list = [block for block,inner in gf_struct] - if with_Sigma: - glist = lambda : [ GfImFreq(indices=inner,mesh=Sigma_imp[0].mesh) for block,inner in gf_struct] - else: - if iw_or_w == "iw": - glist = lambda : [ GfImFreq(indices=inner,beta=beta) for block,inner in gf_struct] - elif iw_or_w == "w": - glist = lambda : [ GfImFreq(indices=inner,window=(mesh[0],mesh[1]),n_points=mesh[2]) for block,inner in gf_struct] - G_upfold = BlockGf(name_list = block_ind_list, block_list = glist(), make_copies = False) - G_upfold.zero() - - if iw_or_w == "iw": - G_upfold << iOmega_n - elif iw_or_w == "w": - G_upfold << Omega + 1j*broadening - - idmat = [numpy.identity(self.n_orbitals[ik,ntoi[sp]],numpy.complex_) for sp in spn] - M = copy.deepcopy(idmat) - for ibl in range(self.n_spin_blocks[self.SO]): - ind = ntoi[spn[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)) - G_upfold -= M - - if with_Sigma: - for icrsh in range(self.n_corr_shells): - for bname,gf in G_upfold: gf -= self.upfold(ik,icrsh,bname,sigma_inc_dc[icrsh][bname],gf) - - G_upfold.invert() - - return G_upfold diff --git a/python/sumk_dft.py b/python/sumk_dft.py index c3f2dc7c..87964943 100644 --- a/python/sumk_dft.py +++ b/python/sumk_dft.py @@ -49,7 +49,7 @@ class SumkDFT: self.symmpar_data = symmpar_data self.bands_data = bands_data self.transp_data = transp_data - self.G_upfold = None + #self.G_latt_iw = None # DEBUG -- remove later self.h_field = h_field # Read input from HDF: @@ -128,7 +128,6 @@ class SumkDFT: if (len(things_to_read) != 0): mpi.report("Loading failed: No %s subgroup in HDF5!"%subgrp) subgroup_present = False value_read = False - del ar # now do the broadcasting: for it in things_to_read: setattr(self,it,mpi.bcast(getattr(self,it))) @@ -163,7 +162,7 @@ class SumkDFT: try: list_to_return.append(ar[subgrp][it]) except: - raise ValueError, "%s not found, and so not loaded."%it + raise ValueError, "load: %s not found, and so not loaded."%it del ar return list_to_return @@ -203,7 +202,7 @@ class SumkDFT: """Local <-> Global rotation of a GF block. direction: 'toLocal' / 'toGlobal' """ - assert ((direction == 'toLocal') or (direction == 'toGlobal')),"Give direction 'toLocal' or 'toGlobal' in rotloc!" + assert ((direction == 'toLocal') or (direction == 'toGlobal')),"rotloc: Give direction 'toLocal' or 'toGlobal'." gf_rotated = gf_to_rotate.copy() @@ -226,76 +225,99 @@ class SumkDFT: return gf_rotated - def lattice_gf_matsubara(self,ik,mu,beta=40,with_Sigma=True): + def lattice_gf(self, ik, mu, iw_or_w="iw", beta=40, broadening=None, mesh=None, with_Sigma=True, with_dc=True): """Calculates the lattice Green function from the DFT hopping and the self energy at k-point number ik and chemical potential mu.""" ntoi = self.spin_names_to_ind[self.SO] spn = self.spin_block_names[self.SO] + if (iw_or_w != "iw") and (iw_or_w != "w"): raise ValueError, "lattice_gf: Implemented only for Re/Im frequency functions." + if not hasattr(self,"Sigma_imp_"+iw_or_w): with_Sigma = False + if broadening is None: + if mesh is None: + broadening = 0.01 + else: # broadening = 2 * \Delta omega, where \Delta omega is the spacing of omega points + broadening = 2.0 * ( (mesh[1]-mesh[0])/(mesh[2]-1) ) - if not hasattr(self,"Sigma_imp"): with_Sigma = False - + # Are we including Sigma? if with_Sigma: - stmp = self.add_dc() - beta = self.Sigma_imp[0].mesh.beta # override beta if Sigma is present + if with_dc: sigma_minus_dc = self.add_dc(iw_or_w) + Sigma_imp = getattr(self,"Sigma_imp_"+iw_or_w) + if iw_or_w == "iw": beta = Sigma_imp[0].mesh.beta # override beta if Sigma_iw is present + else: + if (iw_or_w == "w") and (mesh is None): + raise ValueError, "lattice_gf: Give the mesh=(om_min,om_max,n_points) for the lattice GfReFreq." - # Do we need to set up G_upfold? - set_up_G_upfold = False # assume not - if self.G_upfold is None: # yes if not G_upfold provided - set_up_G_upfold = True - else: # yes if inconsistencies present in existing G_upfold - GFsize = [ gf.N1 for bname,gf in self.G_upfold] - unchangedsize = all( [ self.n_orbitals[ik,ntoi[spn[isp]]] == GFsize[isp] - for isp in range(self.n_spin_blocks[self.SO]) ] ) - if (not unchangedsize) or (self.G_upfold.mesh.beta != beta): set_up_G_upfold = True + # Check if G_latt is present + set_up_G_latt = False # Assume not + if not hasattr(self,"G_latt_"+iw_or_w): + set_up_G_latt = True # Need to create G_latt_(i)w + else: # Check that existing GF is consistent + G_latt = getattr(self,"G_latt_"+iw_or_w) + GFsize = [ gf.N1 for bname,gf in G_latt] + unchangedsize = all( [ self.n_orbitals[ik,ntoi[spn[isp]]]==GFsize[isp] for isp in range(self.n_spin_blocks[self.SO]) ] ) + if not unchangedsize: set_up_G_latt = True + if (iw_or_w == "iw") and (self.G_latt_iw.mesh.beta != beta): set_up_G_latt = True # additional check for ImFreq - # Set up G_upfold - if set_up_G_upfold: + # Set up G_latt + if set_up_G_latt: block_structure = [ range(self.n_orbitals[ik,ntoi[sp]]) for sp in spn ] gf_struct = [ (spn[isp], block_structure[isp]) for isp in range(self.n_spin_blocks[self.SO]) ] block_ind_list = [block for block,inner in gf_struct] - if with_Sigma: - glist = lambda : [ GfImFreq(indices = inner, mesh = self.Sigma_imp[0].mesh) for block,inner in gf_struct] - else: - glist = lambda : [ GfImFreq(indices = inner, beta = beta) for block,inner in gf_struct] - self.G_upfold = BlockGf(name_list = block_ind_list, block_list = glist(), make_copies = False) - self.G_upfold.zero() + if iw_or_w == "iw": + glist = lambda : [ GfImFreq(indices=inner,beta=beta) for block,inner in gf_struct] + elif iw_or_w == "w": + if with_Sigma: + glist = lambda : [ GfReFreq(indices=inner,mesh=Sigma_imp[0].mesh) for block,inner in gf_struct] + else: + glist = lambda : [ GfReFreq(indices=inner,window=(mesh[0],mesh[1]),n_points=mesh[2]) for block,inner in gf_struct] + G_latt = BlockGf(name_list = block_ind_list, block_list = glist(), make_copies = False) + G_latt.zero() + + if iw_or_w == "iw": + G_latt << iOmega_n + elif iw_or_w == "w": + G_latt << Omega + 1j*broadening - self.G_upfold << iOmega_n idmat = [numpy.identity(self.n_orbitals[ik,ntoi[sp]],numpy.complex_) for sp in spn] M = copy.deepcopy(idmat) - for isp in range(self.n_spin_blocks[self.SO]): - ind = ntoi[spn[isp]] + for ibl in range(self.n_spin_blocks[self.SO]): + ind = ntoi[spn[ibl]] n_orb = self.n_orbitals[ik,ind] - M[isp] = self.hopping[ik,ind,0:n_orb,0:n_orb] - (idmat[isp]*mu) - (idmat[isp] * self.h_field * (1-2*isp)) - self.G_upfold -= M + M[ibl] = self.hopping[ik,ind,0:n_orb,0:n_orb] - (idmat[ibl]*mu) - (idmat[ibl] * self.h_field * (1-2*ibl)) + G_latt -= M if with_Sigma: for icrsh in range(self.n_corr_shells): - for bname,gf in self.G_upfold: gf -= self.upfold(ik,icrsh,bname,stmp[icrsh][bname],gf) + for bname,gf in G_latt: gf -= self.upfold(ik,icrsh,bname,sigma_minus_dc[icrsh][bname],gf) - self.G_upfold.invert() + G_latt.invert() + setattr(self,"G_latt_"+iw_or_w,G_latt) - return self.G_upfold + return G_latt def put_Sigma(self, Sigma_imp): - """Puts the impurity self energies for inequivalent atoms into the class, respects the multiplicity of the atoms.""" + """Sets 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_shells, "give exactly one Sigma for each inequivalent corr. shell!" + assert isinstance(Sigma_imp,list), "put_Sigma: 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_shells, "put_Sigma: give exactly one Sigma for each inequivalent corr. shell!" - # init self.Sigma_imp: + # init self.Sigma_imp_(i)w: if all(type(gf) == GfImFreq for bname,gf in Sigma_imp[0]): # Imaginary frequency Sigma: - self.Sigma_imp = [ BlockGf( name_block_generator = [ (block,GfImFreq(indices = inner, mesh = Sigma_imp[0].mesh)) for block,inner in self.gf_struct_sumk[icrsh] ], - make_copies = False) for icrsh in range(self.n_corr_shells) ] + self.Sigma_imp_iw = [ BlockGf( name_block_generator = [ (block,GfImFreq(indices = inner, mesh = Sigma_imp[0].mesh)) + for block,inner in self.gf_struct_sumk[icrsh] ], make_copies = False) + for icrsh in range(self.n_corr_shells) ] + SK_Sigma_imp = self.Sigma_imp_iw elif all(type(gf) == GfReFreq for bname,gf in Sigma_imp[0]): # Real frequency Sigma: - self.Sigma_imp = [ BlockGf( name_block_generator = [ (block,GfReFreq(indices = inner, mesh = Sigma_imp[0].mesh)) for block,inner in self.gf_struct_sumk[icrsh] ], - make_copies = False) for icrsh in range(self.n_corr_shells) ] + self.Sigma_imp_w = [ BlockGf( name_block_generator = [ (block,GfReFreq(indices = inner, mesh = Sigma_imp[0].mesh)) + for block,inner in self.gf_struct_sumk[icrsh] ], make_copies = False) + for icrsh in range(self.n_corr_shells) ] + SK_Sigma_imp = self.Sigma_imp_w else: - raise ValueError, "This type of Sigma is not handled." + raise ValueError, "put_Sigma: This type of Sigma is not handled." # transform the CTQMC blocks to the full matrix: for icrsh in range(self.n_corr_shells): @@ -305,12 +327,12 @@ class SumkDFT: for ind2 in inner: block_sumk,ind1_sumk = self.solver_to_sumk[ish][(block,ind1)] block_sumk,ind2_sumk = self.solver_to_sumk[ish][(block,ind2)] - self.Sigma_imp[icrsh][block_sumk][ind1_sumk,ind2_sumk] << Sigma_imp[ish][block][ind1,ind2] + SK_Sigma_imp[icrsh][block_sumk][ind1_sumk,ind2_sumk] << Sigma_imp[ish][block][ind1,ind2] # rotation from local to global coordinate system: if self.use_rotations: for icrsh in range(self.n_corr_shells): - for bname,gf in self.Sigma_imp[icrsh]: self.Sigma_imp[icrsh][bname] << self.rotloc(icrsh, gf, direction = 'toGlobal') + for bname,gf in SK_Sigma_imp[icrsh]: gf << self.rotloc(icrsh, gf, direction = 'toGlobal') def extract_G_loc(self, mu = None, with_Sigma = True): @@ -321,19 +343,19 @@ class SumkDFT: """ if mu is None: mu = self.chemical_potential - Gloc = [ self.Sigma_imp[icrsh].copy() for icrsh in range(self.n_corr_shells) ] # this list will be returned - for icrsh in range(self.n_corr_shells): Gloc[icrsh].zero() # initialize to zero + Gloc = [ self.Sigma_imp_iw[icrsh].copy() for icrsh in range(self.n_corr_shells) ] # this list will be returned + for icrsh in range(self.n_corr_shells): Gloc[icrsh].zero() # initialize to zero beta = Gloc[0].mesh.beta ikarray = numpy.array(range(self.n_k)) for ik in mpi.slice_array(ikarray): - S = self.lattice_gf_matsubara(ik = ik, mu = mu, with_Sigma = with_Sigma, beta = beta) - S *= self.bz_weights[ik] + G_latt_iw = self.lattice_gf(ik = ik, mu = mu, iw_or_w = "iw", with_Sigma = with_Sigma, beta = beta) + G_latt_iw *= self.bz_weights[ik] for icrsh in range(self.n_corr_shells): tmp = Gloc[icrsh].copy() # init temporary storage - for bname,gf in tmp: tmp[bname] << self.downfold(ik,icrsh,bname,S[bname],gf) + for bname,gf in tmp: tmp[bname] << self.downfold(ik,icrsh,bname,G_latt_iw[bname],gf) Gloc[icrsh] += tmp # collect data from mpi @@ -350,19 +372,18 @@ class SumkDFT: for bname,gf in Gloc[icrsh]: Gloc[icrsh][bname] << self.rotloc(icrsh,gf,direction = 'toLocal') # transform to CTQMC blocks: - Glocret = [ BlockGf( name_block_generator = [ (block,GfImFreq(indices = inner, mesh = Gloc[0].mesh)) for block,inner in self.gf_struct_solver[ish].iteritems() ], + Gloc_inequiv = [ BlockGf( name_block_generator = [ (block,GfImFreq(indices = inner, mesh = Gloc[0].mesh)) for block,inner in self.gf_struct_solver[ish].iteritems() ], make_copies = False) for ish in range(self.n_inequiv_shells) ] for ish in range(self.n_inequiv_shells): - for block,inner in self.gf_struct_solver[ish].iteritems(): for ind1 in inner: for ind2 in inner: block_sumk,ind1_sumk = self.solver_to_sumk[ish][(block,ind1)] block_sumk,ind2_sumk = self.solver_to_sumk[ish][(block,ind2)] - Glocret[ish][block][ind1,ind2] << Gloc[self.inequiv_to_corr[ish]][block_sumk][ind1_sumk,ind2_sumk] + Gloc_inequiv[ish][block][ind1,ind2] << Gloc[self.inequiv_to_corr[ish]][block_sumk][ind1_sumk,ind2_sumk] # return only the inequivalent shells: - return Glocret + return Gloc_inequiv def analyse_block_structure(self, threshold = 0.00001, include_shells = None, dm = None): @@ -451,7 +472,7 @@ class SumkDFT: def density_matrix(self, method = 'using_gf', beta = 40.0): """Calculate density matrices in one of two ways: - if 'using_gf': First get upfolded gf (g_loc is not set up), then density matrix. + if 'using_gf': First get lattice gf (g_loc is not set up), then density matrix. It is useful for Hubbard I, and very quick. No assumption on the hopping structure is made (ie diagonal or not). if 'using_point_integration': Only works for diagonal hopping matrix (true in wien2k). @@ -466,9 +487,9 @@ class SumkDFT: if method == "using_gf": - G_upfold = self.lattice_gf_matsubara(ik = ik, beta = beta, mu = self.chemical_potential) - G_upfold *= self.bz_weights[ik] - dm = G_upfold.density() + G_latt_iw = self.lattice_gf(ik = ik, mu = self.chemical_potential, iw_or_w = "iw", beta = beta) + G_latt_iw *= self.bz_weights[ik] + dm = G_latt_iw.density() MMat = [dm[sp] for sp in self.spin_block_names[self.SO]] elif method == "using_point_integration": @@ -494,10 +515,10 @@ class SumkDFT: for icrsh in range(self.n_corr_shells): for isp, sp in enumerate(self.spin_block_names[self.corr_shells[icrsh]['SO']]): - isp = self.spin_names_to_ind[self.corr_shells[icrsh]['SO']][sp] + ind = 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] + n_orb = self.n_orbitals[ik,ind] + projmat = self.proj_mat[ik,ind,icrsh,0:dim,0:n_orb] if method == "using_gf": dens_mat[icrsh][sp] += numpy.dot( numpy.dot(projmat,MMat[isp]), projmat.transpose().conjugate() ) @@ -507,8 +528,8 @@ class SumkDFT: # get data from nodes: for icrsh in range(self.n_corr_shells): - for bname in dens_mat[icrsh]: - dens_mat[icrsh][bname] = mpi.all_reduce(mpi.world, dens_mat[icrsh][bname], lambda x,y : x+y) + for sp in dens_mat[icrsh]: + dens_mat[icrsh][sp] = mpi.all_reduce(mpi.world, dens_mat[icrsh][sp], lambda x,y : x+y) mpi.barrier() if self.symm_op != 0: dens_mat = self.symmcorr.symmetrize(dens_mat) @@ -516,9 +537,9 @@ class SumkDFT: # Rotate to local coordinate system: if self.use_rotations: for icrsh in range(self.n_corr_shells): - for bl in dens_mat[icrsh]: - if self.rot_mat_time_inv[icrsh] == 1: dens_mat[icrsh][bl] = dens_mat[icrsh][bl].conjugate() - dens_mat[icrsh][bl] = numpy.dot( numpy.dot(self.rot_mat[icrsh].conjugate().transpose(),dens_mat[icrsh][bl]), + for sp in dens_mat[icrsh]: + if self.rot_mat_time_inv[icrsh] == 1: dens_mat[icrsh][sp] = dens_mat[icrsh][sp].conjugate() + dens_mat[icrsh][sp] = numpy.dot( numpy.dot(self.rot_mat[icrsh].conjugate().transpose(),dens_mat[icrsh][sp]), self.rot_mat[icrsh] ) return dens_mat @@ -533,54 +554,41 @@ class SumkDFT: for ish in range(self.n_inequiv_shells): 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): - for ii in eff_atlevels[ish]: - eff_atlevels[ish][ii] *= -self.chemical_potential - - # double counting term: - for ish in range(self.n_inequiv_shells): - for ii in eff_atlevels[ish]: - eff_atlevels[ish][ii] -= self.dc_imp[self.inequiv_to_corr[ish]][ii] + eff_atlevels[ish][sp] *= -self.chemical_potential + eff_atlevels[ish][sp] -= self.dc_imp[self.inequiv_to_corr[ish]][sp] # sum over k: if not hasattr(self,"Hsumk"): # 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]['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]['dim'] + for sp in self.spin_block_names[self.corr_shells[icrsh]['SO']]: + self.Hsumk[icrsh][sp] = numpy.zeros([dim,dim],numpy.complex_) 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] + ind = 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] + n_orb = self.n_orbitals[ik,ind] MMat = numpy.identity(n_orb, numpy.complex_) - MMat = self.hopping[ik,isp,0:n_orb,0:n_orb] - (1-2*isp) * self.h_field * MMat - projmat = self.proj_mat[ik,isp,icrsh,0:dim,0:n_orb] + MMat = self.hopping[ik,ind,0:n_orb,0:n_orb] - (1-2*isp) * self.h_field * MMat + projmat = self.proj_mat[ik,ind,icrsh,0:dim,0:n_orb] self.Hsumk[icrsh][sp] += self.bz_weights[ik] * numpy.dot( numpy.dot(projmat,MMat), projmat.conjugate().transpose() ) - # symmetrisation: if self.symm_op != 0: self.Hsumk = self.symmcorr.symmetrize(self.Hsumk) # Rotate to local coordinate system: if self.use_rotations: for icrsh in range(self.n_corr_shells): - for bl in self.Hsumk[icrsh]: - - if self.rot_mat_time_inv[icrsh] == 1: self.Hsumk[icrsh][bl] = self.Hsumk[icrsh][bl].conjugate() - self.Hsumk[icrsh][bl] = numpy.dot( numpy.dot(self.rot_mat[icrsh].conjugate().transpose(),self.Hsumk[icrsh][bl]) , + for sp in self.Hsumk[icrsh]: + if self.rot_mat_time_inv[icrsh] == 1: self.Hsumk[icrsh][sp] = self.Hsumk[icrsh][sp].conjugate() + self.Hsumk[icrsh][sp] = numpy.dot( numpy.dot(self.rot_mat[icrsh].conjugate().transpose(),self.Hsumk[icrsh][sp]) , self.rot_mat[icrsh] ) # add to matrix: for ish in range(self.n_inequiv_shells): - for bl in eff_atlevels[ish]: - eff_atlevels[ish][bl] += self.Hsumk[self.inequiv_to_corr[ish]][bl] + for sp in eff_atlevels[ish]: + eff_atlevels[ish][sp] += self.Hsumk[self.inequiv_to_corr[ish]][sp] return eff_atlevels @@ -588,6 +596,7 @@ class SumkDFT: def init_dc(self): """ Initialise the double counting terms to have the correct shape.""" + self.dc_imp = [ {} for icrsh in range(self.n_corr_shells)] for icrsh in range(self.n_corr_shells): dim = self.corr_shells[icrsh]['dim'] @@ -672,25 +681,24 @@ class SumkDFT: mpi.report("DC energy = %s"%self.dc_energ[icrsh]) - def add_dc(self): + def add_dc(self,iw_or_w="iw"): """Substracts the double counting term from the impurity self energy.""" # Be careful: Sigma_imp is already in the global coordinate system!! - sres = [s.copy() for s in self.Sigma_imp] + sigma_minus_dc = [s.copy() for s in getattr(self,"Sigma_imp_"+iw_or_w)] for icrsh in range(self.n_corr_shells): - for bname,gf in sres[icrsh]: + for bname,gf in sigma_minus_dc[icrsh]: # Transform dc_imp to global coordinate system dccont = numpy.dot(self.rot_mat[icrsh],numpy.dot(self.dc_imp[icrsh][bname],self.rot_mat[icrsh].conjugate().transpose())) - sres[icrsh][bname] -= dccont + sigma_minus_dc[icrsh][bname] -= dccont - return sres # list of self energies corrected by DC + return sigma_minus_dc def symm_deg_gf(self,gf_to_symm,orb): """Symmetrises a GF for the given degenerate shells self.deg_shells""" for degsh in self.deg_shells[orb]: - #loop over degenerate shells: ss = gf_to_symm[degsh[0]].copy() ss.zero() n_deg = len(degsh) @@ -710,10 +718,8 @@ class SumkDFT: dens = 0.0 ikarray = numpy.array(range(self.n_k)) for ik in mpi.slice_array(ikarray): - - S = self.lattice_gf_matsubara(ik = ik, mu = mu) - dens += self.bz_weights[ik] * S.total_density() - + G_latt_iw = self.lattice_gf(ik = ik, mu = mu, iw_or_w = "iw") + dens += self.bz_weights[ik] * G_latt_iw.total_density() # collect data from mpi: dens = mpi.all_reduce(mpi.world, dens, lambda x,y : x+y) mpi.barrier() @@ -748,7 +754,7 @@ class SumkDFT: def calc_density_correction(self,filename = 'dens_mat.dat'): """ Calculates the density correction in order to feed it back to the DFT calculations.""" - assert type(filename) == StringType, "filename has to be a string!" + assert type(filename) == StringType, "calc_density_correction: filename has to be a string!" ntoi = self.spin_names_to_ind[self.SO] spn = self.spin_block_names[self.SO] @@ -761,12 +767,12 @@ class SumkDFT: ikarray = numpy.array(range(self.n_k)) for ik in mpi.slice_array(ikarray): - S = self.lattice_gf_matsubara(ik = ik, mu = self.chemical_potential) - for bname,gf in S: - deltaN[bname][ik] = S[bname].density() - dens[bname] += self.bz_weights[ik] * S[bname].total_density() + G_latt_iw = self.lattice_gf(ik = ik, mu = self.chemical_potential, iw_or_w = "iw") + for bname,gf in G_latt_iw: + deltaN[bname][ik] = G_latt_iw[bname].density() + dens[bname] += self.bz_weights[ik] * G_latt_iw[bname].total_density() - #put mpi Barrier: + # mpi reduce: for bname in deltaN: for ik in range(self.n_k): deltaN[bname][ik] = mpi.all_reduce(mpi.world, deltaN[bname][ik], lambda x,y : x+y) @@ -783,9 +789,9 @@ class SumkDFT: # write chemical potential (in Rydberg): f.write("%.14f\n"%(self.chemical_potential/self.energy_unit)) if self.SP != 0: f1.write("%.14f\n"%(self.chemical_potential/self.energy_unit)) - # write beta in ryderg-1 - f.write("%.14f\n"%(S.mesh.beta*self.energy_unit)) - if self.SP != 0: f1.write("%.14f\n"%(S.mesh.beta*self.energy_unit)) + # write beta in rydberg-1 + f.write("%.14f\n"%(G_latt_iw.mesh.beta*self.energy_unit)) + if self.SP != 0: f1.write("%.14f\n"%(G_latt_iw.mesh.beta*self.energy_unit)) if self.SP == 0: # no spin-polarization @@ -847,12 +853,10 @@ class SumkDFT: return self.chemical_potential - def find_dc(self,orb,guess,dens_mat,dens_req=None,precision=0.01): + def calc_dc_for_density(self,orb,dc_init,dens_mat,density=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 + If density is given, then DC is set such that the LOCAL charge of orbital + orb coincides with the given density.""" def F(dc): self.calc_dc(dens_mat = dens_mat, U_interact = 0, J_hund = 0, orb = orb, use_dc_value = dc) @@ -861,24 +865,20 @@ class SumkDFT: else: return self.extract_G_loc()[orb].total_density() + if density is None: density = self.density_required - self.charge_below - if dens_req is None: - density = self.density_required - self.charge_below - else: - density = dens_req - - dcnew = dichotomy.dichotomy(function = F, - x_init = guess, y_value = density, + dc = dichotomy.dichotomy(function = F, + x_init = dc_init, y_value = density, precision_on_y = precision, delta_x = 0.5, max_loops = 100, x_name = "Double Counting", y_name= "Total Density", verbosity = 3)[0] - return dcnew + return dc - # Check that the density matrix from projectors (DM = P Pdagger) is correct (ie matches DFT) def check_projectors(self): - + """Calculated the density matrix from projectors (DM = P Pdagger) to check that it is correct and + specifically that it matches DFT.""" dens_mat = [numpy.zeros([self.corr_shells[icrsh]['dim'],self.corr_shells[icrsh]['dim']],numpy.complex_) for icrsh in range(self.n_corr_shells)] @@ -895,27 +895,25 @@ class SumkDFT: if self.use_rotations: for icrsh in range(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]) , + dens_mat[icrsh] = numpy.dot( numpy.dot(self.rot_mat[icrsh].conjugate().transpose(),dens_mat[icrsh]), self.rot_mat[icrsh] ) return dens_mat - # Determine the number of equivalent shells - def sorts_of_atoms(self,lst): + def sorts_of_atoms(self,shells): """ - This routine should determine the number of sorts in the double list lst + Determine the number of inequivalent sorts. """ - sortlst = [ lst[i][1] for i in range(len(lst)) ] - sorts = len(set(sortlst)) - return sorts + sortlst = [ shells[i]['sort'] for i in range(len(shells)) ] + n_sorts = len(set(sortlst)) + return n_sorts - # Determine the number of atoms - def number_of_atoms(self,lst): + def number_of_atoms(self,shells): """ - This routine should determine the number of atoms in the double list lst + Determine the number of inequivalent atoms. """ - atomlst = [ lst[i][0] for i in range(len(lst)) ] - atoms = len(set(atomlst)) - return atoms + atomlst = [ shells[i]['atom'] for i in range(len(shells)) ] + n_atoms = len(set(atomlst)) + return n_atoms diff --git a/python/sumk_dft_tools.py b/python/sumk_dft_tools.py index 2d809ab1..8901e9b2 100644 --- a/python/sumk_dft_tools.py +++ b/python/sumk_dft_tools.py @@ -22,7 +22,6 @@ from types import * import numpy -import pytriqs.utility.dichotomy as dichotomy from pytriqs.gf.local import * import pytriqs.utility.mpi as mpi from symmetry import * @@ -36,7 +35,7 @@ class SumkDFTTools(SumkDFT): parproj_data = 'dft_parproj_input', symmpar_data = 'dft_symmpar_input', bands_data = 'dft_bands_input', transp_data = 'dft_transp_input'): - self.G_upfold_refreq = None + #self.G_latt_w = None # DEBUG -- remove later SumkDFT.__init__(self, hdf_file=hdf_file, h_field=h_field, use_dft_blocks=use_dft_blocks, dft_data=dft_data, symmcorr_data=symmcorr_data, parproj_data=parproj_data, symmpar_data=symmpar_data, bands_data=bands_data, transp_data=transp_data) @@ -60,7 +59,7 @@ class SumkDFTTools(SumkDFT): """Local <-> Global rotation of a GF block. direction: 'toLocal' / 'toGlobal' """ - assert (direction == 'toLocal' or direction == 'toGlobal'),"Give direction 'toLocal' or 'toGlobal' in rotloc!" + assert (direction == 'toLocal' or direction == 'toGlobal'),"rotloc_all: Give direction 'toLocal' or 'toGlobal'." gf_rotated = gf_to_rotate.copy() @@ -82,69 +81,6 @@ class SumkDFTTools(SumkDFT): return gf_rotated - def lattice_gf_realfreq(self, ik, mu, broadening, mesh=None, with_Sigma=True): - """Calculates the lattice Green function on the real frequency axis. If self energy is - present and with_Sigma=True, the mesh is taken from Sigma. Otherwise, the mesh has to be given.""" - - ntoi = self.spin_names_to_ind[self.SO] - spn = self.spin_block_names[self.SO] - - if not hasattr(self,"Sigma_imp"): with_Sigma=False - if with_Sigma: - assert all(type(gf) == GfReFreq for bname,gf 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.G_upfold_refreq is None: - # first setting up of G_upfold_refreq - block_structure = [ range(self.n_orbitals[ik,ntoi[sp]]) for sp in spn ] - gf_struct = [ (spn[isp], block_structure[isp]) for isp in range(self.n_spin_blocks[self.SO]) ] - block_ind_list = [block for block,inner in gf_struct] - if with_Sigma: - glist = lambda : [ GfReFreq(indices = inner, mesh=self.Sigma_imp[0].mesh) for block,inner in gf_struct] - else: - glist = lambda : [ GfReFreq(indices = inner, window=(mesh[0],mesh[1]), n_points=mesh[2]) for block,inner in gf_struct] - self.G_upfold_refreq = BlockGf(name_list = block_ind_list, block_list = glist(),make_copies=False) - self.G_upfold_refreq.zero() - - GFsize = [ gf.N1 for bname,gf in self.G_upfold_refreq] - unchangedsize = all( [ self.n_orbitals[ik,ntoi[spn[isp]]] == GFsize[isp] - for isp in range(self.n_spin_blocks[self.SO]) ] ) - - if not unchangedsize: - block_structure = [ range(self.n_orbitals[ik,ntoi[sp]]) for sp in spn ] - gf_struct = [ (spn[isp], block_structure[isp]) for isp in range(self.n_spin_blocks[self.SO]) ] - block_ind_list = [block for block,inner in gf_struct] - if with_Sigma: - glist = lambda : [ GfReFreq(indices = inner, mesh =self.Sigma_imp[0].mesh) for block,inner in gf_struct] - else: - glist = lambda : [ GfReFreq(indices = inner, window=(mesh[0],mesh[1]), n_points=mesh[2]) for block,inner in gf_struct] - self.G_upfold_refreq = BlockGf(name_list = block_ind_list, block_list = glist(),make_copies=False) - self.G_upfold_refreq.zero() - - idmat = [numpy.identity(self.n_orbitals[ik,ntoi[sp]],numpy.complex_) for sp in spn] - - self.G_upfold_refreq << Omega + 1j*broadening - M = copy.deepcopy(idmat) - for isp in range(self.n_spin_blocks[self.SO]): - ind = ntoi[spn[isp]] - n_orb = self.n_orbitals[ik,ind] - M[isp] = self.hopping[ik,ind,0:n_orb,0:n_orb] - (idmat[isp]*mu) - (idmat[isp] * self.h_field * (1-2*isp)) - self.G_upfold_refreq -= M - - if with_Sigma: - tmp = self.G_upfold_refreq.copy() # init temporary storage - for icrsh in range(self.n_corr_shells): - for bname,gf in tmp: tmp[bname] << self.upfold(ik,icrsh,bname,stmp[icrsh][bname],gf) - self.G_upfold_refreq -= tmp # adding to the upfolded GF - - self.G_upfold_refreq.invert() - - return self.G_upfold_refreq - - - def check_input_dos(self, om_min, om_max, n_om, beta=10, broadening=0.01): @@ -153,16 +89,16 @@ class SumkDFTTools(SumkDFT): for i in range(n_om): om_mesh[i] = om_min + delta_om * i DOS = {} - for bn in self.spin_block_names[self.SO]: - DOS[bn] = numpy.zeros([n_om],numpy.float_) + for sp in self.spin_block_names[self.SO]: + DOS[sp] = numpy.zeros([n_om],numpy.float_) 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]]['SO']]: + for sp 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_) + DOSproj[ish][sp] = numpy.zeros([n_om],numpy.float_) + DOSproj_orb[ish][sp] = numpy.zeros([n_om,dim,dim],numpy.float_) # init: Gloc = [] @@ -174,18 +110,18 @@ class SumkDFTTools(SumkDFT): for ik in range(self.n_k): - G_upfold=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.bz_weights[ik] + G_latt_w=self.lattice_gf(ik=ik,mu=self.chemical_potential,iw_or_w="w",broadening=broadening,mesh=(om_min,om_max,n_om),with_Sigma=False) + G_latt_w *= self.bz_weights[ik] # non-projected DOS for iom in range(n_om): - for bname,gf in G_upfold: + for bname,gf in G_latt_w: asd = gf.data[iom,:,:].imag.trace()/(-3.1415926535) DOS[bname][iom] += asd for icrsh in range(self.n_corr_shells): tmp = Gloc[icrsh].copy() - for bname,gf in tmp: tmp[bname] << self.downfold(ik,icrsh,bname,G_upfold[bname],gf) # downfolding G + for bname,gf in tmp: tmp[bname] << self.downfold(ik,icrsh,bname,G_latt_w[bname],gf) # downfolding G Gloc[icrsh] += tmp @@ -205,21 +141,21 @@ class SumkDFTTools(SumkDFT): # output: if mpi.is_master_node(): - for bn in self.spin_block_names[self.SO]: - f=open('DOS%s.dat'%bn, 'w') - for i in range(n_om): f.write("%s %s\n"%(om_mesh[i],DOS[bn][i])) + for sp in self.spin_block_names[self.SO]: + f=open('DOS%s.dat'%sp, 'w') + for i in range(n_om): f.write("%s %s\n"%(om_mesh[i],DOS[sp][i])) f.close() for ish in range(self.n_inequiv_shells): - f=open('DOS%s_proj%s.dat'%(bn,ish),'w') - for i in range(n_om): f.write("%s %s\n"%(om_mesh[i],DOSproj[ish][bn][i])) + f=open('DOS%s_proj%s.dat'%(sp,ish),'w') + for i in range(n_om): f.write("%s %s\n"%(om_mesh[i],DOSproj[ish][sp][i])) f.close() 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' + Fname = 'DOS'+sp+'_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])) + for iom in range(n_om): f.write("%s %s\n"%(om_mesh[iom],DOSproj_orb[ish][sp][iom,i,j])) f.close() @@ -239,7 +175,7 @@ class SumkDFTTools(SumkDFT): def dos_partial(self,broadening=0.01): """calculates the orbitally-resolved DOS""" - assert hasattr(self,"Sigma_imp"), "Set Sigma First!!" + assert hasattr(self,"Sigma_imp_w"), "dos_partial: Set Sigma_imp_w first." value_read = self.read_parproj_input_from_hdf() if not value_read: return value_read @@ -248,41 +184,41 @@ class SumkDFTTools(SumkDFT): mu = self.chemical_potential 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)] + Gproj = [BlockGf(name_block_generator = [ (block,GfReFreq(indices = inner, mesh = self.Sigma_imp_w[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() - Msh = [x.real for x in self.Sigma_imp[0].mesh] - n_om = len(Msh) + mesh = [x.real for x in self.Sigma_imp_w[0].mesh] + n_om = len(mesh) DOS = {} - for bn in self.spin_block_names[self.SO]: - DOS[bn] = numpy.zeros([n_om],numpy.float_) + for sp in self.spin_block_names[self.SO]: + DOS[sp] = numpy.zeros([n_om],numpy.float_) DOSproj = [ {} for ish in range(self.n_shells) ] 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]: + for sp in self.spin_block_names[self.SO]: 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_) + DOSproj[ish][sp] = numpy.zeros([n_om],numpy.float_) + DOSproj_orb[ish][sp] = numpy.zeros([n_om,dim,dim],numpy.float_) ikarray=numpy.array(range(self.n_k)) for ik in mpi.slice_array(ikarray): - S = self.lattice_gf_realfreq(ik=ik,mu=mu,broadening=broadening) - S *= self.bz_weights[ik] + G_latt_w = self.lattice_gf(ik=ik,mu=mu,iw_or_w="w",broadening=broadening) + G_latt_w *= self.bz_weights[ik] # non-projected DOS for iom in range(n_om): - for bname,gf in S: DOS[bname][iom] += gf.data[iom,:,:].imag.trace()/(-3.1415926535) + for bname,gf in G_latt_w: DOS[bname][iom] += gf.data[iom,:,:].imag.trace()/(-3.1415926535) #projected DOS: for ish in range(self.n_shells): tmp = Gproj[ish].copy() for ir in range(self.n_parproj[ish]): - for bname,gf in tmp: tmp[bname] << self.downfold_pc(ik,ir,ish,bname,S[bname],gf) + for bname,gf in tmp: tmp[bname] << self.downfold_pc(ik,ir,ish,bname,G_latt_w[bname],gf) Gproj[ish] += tmp # collect data from mpi: @@ -307,22 +243,22 @@ class SumkDFTTools(SumkDFT): if mpi.is_master_node(): # output to files - for bn in self.spin_block_names[self.SO]: - f=open('./DOScorr%s.dat'%bn, 'w') - for i in range(n_om): f.write("%s %s\n"%(Msh[i],DOS[bn][i])) + for sp in self.spin_block_names[self.SO]: + f=open('./DOScorr%s.dat'%sp, 'w') + for i in range(n_om): f.write("%s %s\n"%(mesh[i],DOS[sp][i])) f.close() # partial for ish in range(self.n_shells): - f=open('DOScorr%s_proj%s.dat'%(bn,ish),'w') - for i in range(n_om): f.write("%s %s\n"%(Msh[i],DOSproj[ish][bn][i])) + f=open('DOScorr%s_proj%s.dat'%(sp,ish),'w') + for i in range(n_om): f.write("%s %s\n"%(mesh[i],DOSproj[ish][sp][i])) f.close() 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' + Fname = './DOScorr'+sp+'_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])) + for iom in range(n_om): f.write("%s %s\n"%(mesh[iom],DOSproj_orb[ish][sp][iom,i,j])) f.close() @@ -332,7 +268,7 @@ class SumkDFTTools(SumkDFT): """ Calculates the correlated band structure with a real-frequency self energy. ATTENTION: Many things from the original input file are overwritten!!!""" - assert hasattr(self,"Sigma_imp"), "Set Sigma First!!" + assert hasattr(self,"Sigma_imp_w"), "spaghettis: Set Sigma_imp_w first." things_to_read = ['n_k','n_orbitals','proj_mat','hopping','n_parproj','proj_mat_pc'] value_read = self.read_input_from_hdf(subgrp=self.bands_data,things_to_read=things_to_read) if not value_read: return value_read @@ -370,12 +306,12 @@ class SumkDFTTools(SumkDFT): spn = self.spin_block_names[self.SO] # init DOS: - M = [x.real for x in self.Sigma_imp[0].mesh] - n_om = len(M) + mesh = [x.real for x in self.Sigma_imp_w[0].mesh] + n_om = len(mesh) if plot_range is None: - om_minplot = M[0]-0.001 - om_maxplot = M[n_om-1] + 0.001 + om_minplot = mesh[0]-0.001 + om_maxplot = mesh[n_om-1] + 0.001 else: om_minplot = plot_range[0] om_maxplot = plot_range[1] @@ -395,20 +331,20 @@ class SumkDFTTools(SumkDFT): if not ishell is None: 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 = BlockGf(name_block_generator = [ (block,GfReFreq(indices = inner, mesh = self.Sigma_imp_w[0].mesh)) for block,inner in GFStruct_proj ], make_copies = False) Gproj.zero() for ik in range(self.n_k): - S = self.lattice_gf_realfreq(ik=ik,mu=mu,broadening=broadening) + G_latt_w = self.lattice_gf(ik=ik,mu=mu,iw_or_w="w",broadening=broadening) if ishell is None: # non-projected A(k,w) for iom in range(n_om): - if (M[iom] > om_minplot) and (M[iom] < om_maxplot): + if (mesh[iom] > om_minplot) and (mesh[iom] < om_maxplot): if fermi_surface: - for bname,gf in S: Akw[bname][ik,0] += gf.data[iom,:,:].imag.trace()/(-3.1415926535) * (M[1]-M[0]) + for bname,gf in G_latt_w: Akw[bname][ik,0] += gf.data[iom,:,:].imag.trace()/(-3.1415926535) * (mesh[1]-mesh[0]) else: - for bname,gf in S: Akw[bname][ik,iom] += gf.data[iom,:,:].imag.trace()/(-3.1415926535) + for bname,gf in G_latt_w: Akw[bname][ik,iom] += gf.data[iom,:,:].imag.trace()/(-3.1415926535) Akw[bname][ik,iom] += ik*shift # shift Akw for plotting in xmgrace -- REMOVE @@ -417,7 +353,7 @@ class SumkDFTTools(SumkDFT): Gproj.zero() tmp = Gproj.copy() for ir in range(self.n_parproj[ishell]): - for bname,gf in tmp: tmp[bname] << self.downfold_pc(ik,ir,ishell,bname,S[bname],gf) + for bname,gf in tmp: tmp[bname] << self.downfold_pc(ik,ir,ishell,bname,G_latt_w[bname],gf) Gproj += tmp # FIXME NEED TO READ IN ROTMAT_ALL FROM PARPROJ SUBGROUP, REPLACE ROTLOC WITH ROTLOC_ALL @@ -427,7 +363,7 @@ class SumkDFTTools(SumkDFT): # for bname,gf in Gproj: Gproj[bname] << self.rotloc(0,gf,direction='toLocal') for iom in range(n_om): - if (M[iom] > om_minplot) and (M[iom] < om_maxplot): + if (mesh[iom] > om_minplot) and (mesh[iom] < om_maxplot): 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) @@ -458,13 +394,13 @@ class SumkDFTTools(SumkDFT): f.write('%s %s\n'%(ik,Akw[ibn][ik,0])) else: for iom in range(n_om): - if (M[iom] > om_minplot) and (M[iom] < om_maxplot): + if (mesh[iom] > om_minplot) and (mesh[iom] < om_maxplot): if invert_Akw: Akw[ibn][ik,iom] = 1.0/(minAkw-maxAkw)*(Akw[ibn][ik,iom] - maxAkw) if shift > 0.0001: - f.write('%s %s\n'%(M[iom],Akw[ibn][ik,iom])) + f.write('%s %s\n'%(mesh[iom],Akw[ibn][ik,iom])) else: - f.write('%s %s %s\n'%(ik,M[iom],Akw[ibn][ik,iom])) + f.write('%s %s %s\n'%(ik,mesh[iom],Akw[ibn][ik,iom])) f.write('\n') @@ -482,13 +418,13 @@ class SumkDFTTools(SumkDFT): for ik in range(self.n_k): for iom in range(n_om): - if (M[iom] > om_minplot) and (M[iom] < om_maxplot): + if (mesh[iom] > om_minplot) and (mesh[iom] < om_maxplot): if invert_Akw: Akw[ibn][ish,ik,iom] = 1.0/(minAkw-maxAkw)*(Akw[ibn][ish,ik,iom] - maxAkw) if shift > 0.0001: - f.write('%s %s\n'%(M[iom],Akw[ibn][ish,ik,iom])) + f.write('%s %s\n'%(mesh[iom],Akw[ibn][ish,ik,iom])) else: - f.write('%s %s %s\n'%(ik,M[iom],Akw[ibn][ish,ik,iom])) + f.write('%s %s %s\n'%(ik,mesh[iom],Akw[ibn][ish,ik,iom])) f.write('\n') @@ -511,10 +447,10 @@ class SumkDFTTools(SumkDFT): mu = self.chemical_potential 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) + if hasattr(self,"Sigma_imp_iw"): + Gproj = [BlockGf(name_block_generator = [ (block,GfImFreq(indices = inner, mesh = self.Sigma_imp_iw[0].mesh)) for block,inner in GFStruct_proj[ish] ], make_copies = False) for ish in range(self.n_shells)] - beta = self.Sigma_imp[0].mesh.beta + beta = self.Sigma_imp_iw[0].mesh.beta else: Gproj = [BlockGf(name_block_generator = [ (block,GfImFreq(indices = inner, beta = beta)) for block,inner in GFStruct_proj[ish] ], make_copies = False) for ish in range(self.n_shells)] @@ -524,13 +460,13 @@ class SumkDFTTools(SumkDFT): ikarray=numpy.array(range(self.n_k)) for ik in mpi.slice_array(ikarray): - S = self.lattice_gf_matsubara(ik=ik,mu=mu,beta=beta) - S *= self.bz_weights[ik] + G_latt_iw = self.lattice_gf(ik=ik,mu=mu,iw_or_w="iw",beta=beta) + G_latt_iw *= self.bz_weights[ik] for ish in range(self.n_shells): tmp = Gproj[ish].copy() for ir in range(self.n_parproj[ish]): - for bname,gf in tmp: tmp[bname] << self.downfold_pc(ik,ir,ish,bname,S[bname],gf) + for bname,gf in tmp: tmp[bname] << self.downfold_pc(ik,ir,ish,bname,G_latt_iw[bname],gf) Gproj[ish] += tmp #collect data from mpi: @@ -619,7 +555,7 @@ class SumkDFTTools(SumkDFT): # Define mesh for Greens function and the used energy range if (with_Sigma == False): - self.omega = numpy.array([round(x.real,12) for x in self.Sigma_imp[0].mesh]) + self.omega = numpy.array([round(x.real,12) for x in self.Sigma_imp_w[0].mesh]) mu = self.chemical_potential n_om = len(self.omega) print "Using omega mesh provided by Sigma." @@ -632,15 +568,15 @@ class SumkDFTTools(SumkDFT): # Truncate Sigma to given omega window for icrsh in range(self.n_corr_shells): - Sigma_save = self.Sigma_imp[icrsh].copy() + Sigma_save = self.Sigma_imp_w[icrsh].copy() spn = self.spin_block_names[self.corr_shells[icrsh]['SO']] glist = lambda : [ GfReFreq(indices = inner, window=(self.omega[0], self.omega[-1]),n_points=n_om) for block, inner in self.gf_struct_sumk[icrsh]] - self.Sigma_imp[icrsh] = BlockGf(name_list = spn, block_list = glist(),make_copies=False) - for i,g in self.Sigma_imp[icrsh]: + self.Sigma_imp_w[icrsh] = BlockGf(name_list = spn, block_list = glist(),make_copies=False) + for i,g in self.Sigma_imp_w[icrsh]: for iL in g.indices: for iR in g.indices: for iom in xrange(n_om): - g.data[iom,iL,iR]= Sigma_save[i].data[ioffset+iom,iL,iR] + g.data[iom,iL,iR] = Sigma_save[i].data[ioffset+iom,iL,iR] # FIXME IS THIS CLEAN? MANIPULATING data DIRECTLY? else: assert n_om is not None, "Number of omega points (n_om) needed!" @@ -684,8 +620,8 @@ class SumkDFTTools(SumkDFT): ikarray = numpy.array(range(self.n_k)) for ik in mpi.slice_array(ikarray): - unchangesize = all([ self.n_orbitals[ik][isp] == mupat[isp].shape[0] for isp in range(n_inequiv_spin_blocks)]) - if (not unchangesize): + unchangedsize = all([ self.n_orbitals[ik][isp] == mupat[isp].shape[0] for isp in range(n_inequiv_spin_blocks)]) + if not unchangedsize: # recontruct green functions. S = BlockGf(name_block_generator=[(bln[isp], GfReFreq(indices=range(self.n_orbitals[ik][isp]), window = (self.omega[0], self.omega[-1]), n_points = n_om)) for isp in range(n_inequiv_spin_blocks) ], make_copies=False) @@ -707,10 +643,10 @@ class SumkDFTTools(SumkDFT): if (with_Sigma == False): tmp = S.copy() # init temporary storage # form self energy from impurity self energy and double counting term. - stmp = self.add_dc() - ## substract self energy + sigma_minus_dc = self.add_dc(iw_or_w="w") + # substract self energy 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, sigma_minus_dc[icrsh][sig], gf) S -= tmp S.invert() @@ -786,7 +722,7 @@ class SumkDFTTools(SumkDFT): volcc, volpc = self.cellvolume(self.latticetype, self.latticeconstants, self.latticeangles) - n_direction,n_q,n_w= self.Pw_optic.shape + n_direction, n_q, n_w= self.Pw_optic.shape omegaT = self.omega * beta A0 = numpy.empty((n_direction,n_q), dtype=numpy.float_) q_0 = False diff --git a/python/symmetry.py b/python/symmetry.py index 03842fc3..8b50e457 100644 --- a/python/symmetry.py +++ b/python/symmetry.py @@ -39,7 +39,7 @@ class Symmetry: SP: Flag for spin polarisation. """ - assert type(hdf_file) == StringType, "hdf_file must be a filename" + assert type(hdf_file) == StringType, "Symmetry: hdf_file must be a filename." 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) @@ -73,8 +73,8 @@ class Symmetry: def symmetrize(self,obj): - assert isinstance(obj,list), "symmetry: obj has to be a list of objects." - assert len(obj) == self.n_orbits, "symmetry: obj has to be a list of the same length as defined in the init." + assert isinstance(obj,list), "symmetrize: obj has to be a list of objects." + assert len(obj) == self.n_orbits, "symmetrize: obj has to be a list of the same length as defined in the init." 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!