From 5ad2acfee6755a2b4fae16b97a361363f50e18c2 Mon Sep 17 00:00:00 2001 From: Priyanka Seth Date: Mon, 17 Nov 2014 10:48:04 +0100 Subject: [PATCH] More minor changes --- python/sumk_lda.py | 501 +++++++++++++++++++-------------------- python/sumk_lda_tools.py | 58 ++--- 2 files changed, 270 insertions(+), 289 deletions(-) diff --git a/python/sumk_lda.py b/python/sumk_lda.py index 41f53824..fb8965bd 100644 --- a/python/sumk_lda.py +++ b/python/sumk_lda.py @@ -68,12 +68,12 @@ class SumkLDA: # convert spin_block_names to indices -- if spin polarized, differentiate up and down blocks self.spin_names_to_ind = [{}, {}] for iso in range(2): # SO = 0 or 1 - for ibl in range(self.n_spin_blocks[iso]): - self.spin_names_to_ind[iso][self.spin_block_names[iso][ibl]] = ibl * self.SP + for isp in range(self.n_spin_blocks[iso]): + self.spin_names_to_ind[iso][self.spin_block_names[iso][isp]] = isp * self.SP # GF structure used for the local things in the k sums # Most general form allowing for all hybridisation, i.e. largest blocks possible - self.gf_struct_sumk = [ [ (b, range( self.corr_shells[icrsh][3])) for b in self.spin_block_names[self.corr_shells[icrsh][4]] ] + self.gf_struct_sumk = [ [ (sp, range( self.corr_shells[icrsh][3])) for sp in self.spin_block_names[self.corr_shells[icrsh][4]] ] for icrsh in range(self.n_corr_shells) ] #----- @@ -83,8 +83,8 @@ class SumkLDA: optional_things = optional_things) if (not self.subgroup_present) or (not self.value_read['gf_struct_solver']): # No gf_struct was stored in HDF, so first set a standard one: - self.gf_struct_solver = [ dict([ (b, range(self.corr_shells[self.inequiv_to_corr[ish]][3]) ) - for b in self.spin_block_names[self.corr_shells[self.inequiv_to_corr[ish]][4]] ]) + self.gf_struct_solver = [ dict([ (sp, range(self.corr_shells[self.inequiv_to_corr[ish]][3]) ) + for sp in self.spin_block_names[self.corr_shells[self.inequiv_to_corr[ish]][4]] ]) for ish in range(self.n_inequiv_shells) ] # Set standard (identity) maps from gf_struct_sumk <-> gf_struct_solver @@ -247,7 +247,7 @@ class SumkLDA: and chemical potential mu.""" ntoi = self.spin_names_to_ind[self.SO] - bln = self.spin_block_names[self.SO] + spn = self.spin_block_names[self.SO] if not hasattr(self,"Sigma_imp"): with_Sigma = False @@ -261,29 +261,29 @@ class SumkLDA: 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[bln[ibl]]] == GFsize[ibl] - for ibl in range(self.n_spin_blocks[self.SO]) ] ) + 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 # Set up G_upfold if set_up_G_upfold: - block_structure = [ range(self.n_orbitals[ik,ntoi[b]]) for b in bln ] - gf_struct = [ (bln[ibl], block_structure[ibl]) for ibl in range(self.n_spin_blocks[self.SO]) ] + 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 = BlockGf(name_list = block_ind_list, block_list = glist(), make_copies = False) self.G_upfold.zero() self.G_upfold << iOmega_n - idmat = [numpy.identity(self.n_orbitals[ik,ntoi[bl]],numpy.complex_) for bl in bln] + 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[bln[ibl]] + for isp in range(self.n_spin_blocks[self.SO]): + ind = ntoi[spn[isp]] 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)) + 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 if with_Sigma: @@ -295,81 +295,94 @@ class SumkLDA: return self.G_upfold - 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. - 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). - """ - dens_mat = [ {} for icrsh in range(self.n_corr_shells)] + 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_shells, "give exactly one Sigma for each inequivalent corr. shell!" + + # init self.Sigma_imp: + 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) ] + 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) ] + else: + raise ValueError, "This type of Sigma is not handled." + + # transform the CTQMC blocks to the full matrix: for icrsh in range(self.n_corr_shells): - for bl in self.spin_block_names[self.corr_shells[icrsh][4]]: - dens_mat[icrsh][bl] = numpy.zeros([self.corr_shells[icrsh][3],self.corr_shells[icrsh][3]], numpy.complex_) + ish = self.corr_to_inequiv[icrsh] # ish is the index of the inequivalent shell corresponding to icrsh + 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)] + self.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') + + + def extract_G_loc(self, mu = None, with_Sigma = True): + """ + Extracts the local downfolded Green function at the chemical potential of the class. + At the end, the local G is rotated from the global coordinate system to the local system. + if with_Sigma = False: Sigma is not included => non-interacting local GF + """ + + 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 + beta = Gloc[0].mesh.beta ikarray = numpy.array(range(self.n_k)) for ik in mpi.slice_array(ikarray): - 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() - MMat = [dm[bl] for bl in self.spin_block_names[self.SO]] - - elif method == "using_point_integration": - - ntoi = self.spin_names_to_ind[self.SO] - bln = self.spin_block_names[self.SO] - unchangedsize = all( [self.n_orbitals[ik,ntoi[bl]] == self.n_orbitals[0,ntoi[bl]] for bl in bln] ) - if unchangedsize: - dim = self.n_orbitals[0,ntoi[bl]] - else: - dim = self.n_orbitals[ik,ntoi[bl]] - MMat = [numpy.zeros( [dim,dim], numpy.complex_) for bl in bln] - - 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: # only works for diagonal hopping matrix (true in wien2k) - MMat[ibl][inu,inu] = 1.0 - else: - MMat[ibl][inu,inu] = 0.0 + S = self.lattice_gf_matsubara(ik = ik, mu = mu, with_Sigma = with_Sigma, beta = beta) + S *= self.bz_weights[ik] for icrsh in range(self.n_corr_shells): - for ibl, bn in enumerate(self.spin_block_names[self.corr_shells[icrsh][4]]): - isp = self.spin_names_to_ind[self.corr_shells[icrsh][4]][bn] - dim = self.corr_shells[icrsh][3] - n_orb = self.n_orbitals[ik,isp] - projmat = self.proj_mat[ik,isp,icrsh,0:dim,0:n_orb] - if method == "using_gf": - dens_mat[icrsh][bn] += numpy.dot( numpy.dot(projmat,MMat[ibl]), - projmat.transpose().conjugate() ) - elif method == "using_point_integration": - dens_mat[icrsh][bn] += self.bz_weights[ik] * numpy.dot( numpy.dot(projmat,MMat[ibl]) , - projmat.transpose().conjugate() ) + tmp = Gloc[icrsh].copy() # init temporary storage + for bname,gf in tmp: tmp[bname] << self.downfold(ik,icrsh,bname,S[bname],gf) + Gloc[icrsh] += tmp - # 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) + # collect data from mpi + for icrsh in range(self.n_corr_shells): Gloc[icrsh] << mpi.all_reduce(mpi.world, Gloc[icrsh], lambda x,y : x+y) mpi.barrier() - if self.symm_op != 0: dens_mat = self.symmcorr.symmetrize(dens_mat) + # Gloc[:] is now the sum over k projected to the local orbitals. + # here comes the symmetrisation, if needed: + if self.symm_op != 0: Gloc = self.symmcorr.symmetrize(Gloc) - # Rotate to local coordinate system: + # Gloc is rotated to the local coordinate system: if self.use_rotations: for icrsh in range(self.n_corr_shells): - for bn in dens_mat[icrsh]: - if self.rot_mat_time_inv[icrsh] == 1: dens_mat[icrsh][bn] = dens_mat[icrsh][bn].conjugate() - dens_mat[icrsh][bn] = numpy.dot( numpy.dot(self.rot_mat[icrsh].conjugate().transpose(),dens_mat[icrsh][bn]), - self.rot_mat[icrsh] ) + for bname,gf in Gloc[icrsh]: Gloc[icrsh][bname] << self.rotloc(icrsh,gf,direction = 'toLocal') - return dens_mat + # 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() ], + 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] + + # return only the inequivalent shells: + return Glocret def analyse_block_structure(self, threshold = 0.00001, include_shells = None, dm = None): - """ Determines the Green function block structure from simple point integration.""" + """ Determines the Green's function block structure from simple point integration.""" self.gf_struct_solver = [ {} for ish in range(self.n_inequiv_shells) ] self.sumk_to_solver = [ {} for ish in range(self.n_inequiv_shells) ] @@ -382,10 +395,9 @@ class SumkLDA: if include_shells is None: include_shells = range(self.n_inequiv_shells) for ish in include_shells: - block_ind_list = [ block for block,inner in self.gf_struct_sumk[self.inequiv_to_corr[ish]] ] - for block in block_ind_list: - dm = dens_mat[ish][block] - dmbool = (abs(dm) > threshold) # gives an index list of entries larger that threshold + + for sp in self.spin_block_names[self.corr_shells[self.inequiv_to_corr[ish]][4]]: + dmbool = (abs(dens_mat[ish][sp]) > threshold) # gives an index list of entries larger that threshold # Determine off-diagonal entries in upper triangular part of density matrix offdiag = [] @@ -393,9 +405,9 @@ class SumkLDA: for j in range(i+1,len(dmbool)): if dmbool[i,j]: offdiag.append([i,j]) + # Determine the number of non-hybridising blocks in the gf num_blocs = len(dmbool) blocs = [ [i] for i in range(num_blocs) ] - for i in range(len(offdiag)): for j in range(len(blocs[offdiag[i][1]])): blocs[offdiag[i][0]].append(blocs[offdiag[i][1]][j]) del blocs[offdiag[i][1]] @@ -407,17 +419,18 @@ class SumkLDA: offdiag[j].sort() num_blocs -= 1 + # Set the gf_struct for the solver accordingly for i in range(num_blocs): blocs[i].sort() - self.gf_struct_solver[ish].update( [('%s_%s'%(block,i),range(len(blocs[i])))] ) + self.gf_struct_solver[ish].update( [('%s_%s'%(sp,i),range(len(blocs[i])))] ) # Construct sumk_to_solver taking (sumk_block, sumk_index) --> (solver_block, solver_inner) # and solver_to_sumk taking (solver_block, solver_inner) --> (sumk_block, sumk_index) for i in range(num_blocs): for j in range(len(blocs[i])): - block_sumk = block + block_sumk = sp inner_sumk = blocs[i][j] - block_solv = '%s_%s'%(block,i) + block_solv = '%s_%s'%(sp,i) inner_solv = j self.sumk_to_solver[ish][(block_sumk,inner_sumk)] = (block_solv,inner_solv) self.solver_to_sumk[ish][(block_solv,inner_solv)] = (block_sumk,inner_sumk) @@ -457,16 +470,78 @@ class SumkLDA: return dens_mat - def symm_deg_gf(self,gf_to_symm,orb): - """Symmetrises a GF for the given degenerate shells self.deg_shells""" + 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. + 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). + """ + dens_mat = [ {} for icrsh in range(self.n_corr_shells)] + for icrsh in range(self.n_corr_shells): + for sp in self.spin_block_names[self.corr_shells[icrsh][4]]: + dens_mat[icrsh][sp] = numpy.zeros([self.corr_shells[icrsh][3],self.corr_shells[icrsh][3]], numpy.complex_) + + ikarray = numpy.array(range(self.n_k)) + for ik in mpi.slice_array(ikarray): + + 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() + MMat = [dm[sp] for sp in self.spin_block_names[self.SO]] + + elif method == "using_point_integration": + + ntoi = self.spin_names_to_ind[self.SO] + spn = self.spin_block_names[self.SO] + unchangedsize = all( [self.n_orbitals[ik,ntoi[sp]] == self.n_orbitals[0,ntoi[sp]] for sp in spn] ) + if unchangedsize: + dim = self.n_orbitals[0,ntoi[sp]] + else: + dim = self.n_orbitals[ik,ntoi[sp]] + MMat = [numpy.zeros( [dim,dim], numpy.complex_) for sp in spn] + + for isp, sp in enumerate(spn): + ind = ntoi[sp] + for inu in range(self.n_orbitals[ik,ind]): + if (self.hopping[ik,ind,inu,inu] - self.h_field*(1-2*isp)) < 0.0: # only works for diagonal hopping matrix (true in wien2k) + MMat[isp][inu,inu] = 1.0 + else: + MMat[isp][inu,inu] = 0.0 + + for icrsh in range(self.n_corr_shells): + for isp, sp in enumerate(self.spin_block_names[self.corr_shells[icrsh][4]]): + isp = self.spin_names_to_ind[self.corr_shells[icrsh][4]][sp] + dim = self.corr_shells[icrsh][3] + n_orb = self.n_orbitals[ik,isp] + projmat = self.proj_mat[ik,isp,icrsh,0:dim,0:n_orb] + if method == "using_gf": + dens_mat[icrsh][sp] += numpy.dot( numpy.dot(projmat,MMat[isp]), + projmat.transpose().conjugate() ) + elif method == "using_point_integration": + dens_mat[icrsh][sp] += self.bz_weights[ik] * numpy.dot( numpy.dot(projmat,MMat[isp]) , + projmat.transpose().conjugate() ) + + # 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) + mpi.barrier() + + if self.symm_op != 0: dens_mat = self.symmcorr.symmetrize(dens_mat) + + # 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]), + self.rot_mat[icrsh] ) + + return dens_mat - for degsh in self.deg_shells[orb]: - #loop over degenerate shells: - ss = gf_to_symm[degsh[0]].copy() - 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 simple dft input, get crystal field splittings. def eff_atomic_levels(self): @@ -475,8 +550,8 @@ class SumkLDA: # define matrices for inequivalent shells: eff_atlevels = [ {} for ish in range(self.n_inequiv_shells) ] for ish in range(self.n_inequiv_shells): - for bn in self.spin_block_names[self.corr_shells[self.inequiv_to_corr[ish]][4]]: - eff_atlevels[ish][bn] = numpy.identity(self.corr_shells[self.inequiv_to_corr[ish]][3], numpy.complex_) + for sp in self.spin_block_names[self.corr_shells[self.inequiv_to_corr[ish]][4]]: + eff_atlevels[ish][sp] = numpy.identity(self.corr_shells[self.inequiv_to_corr[ish]][3], numpy.complex_) # Chemical Potential: for ish in range(self.n_inequiv_shells): @@ -492,20 +567,20 @@ class SumkLDA: # 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 bn in self.spin_block_names[self.corr_shells[icrsh][4]]: + for sp in self.spin_block_names[self.corr_shells[icrsh][4]]: dim = self.corr_shells[icrsh][3] #*(1+self.corr_shells[icrsh][4]) - self.Hsumk[icrsh][bn] = numpy.zeros([dim,dim],numpy.complex_) + self.Hsumk[icrsh][sp] = numpy.zeros([dim,dim],numpy.complex_) for icrsh in range(self.n_corr_shells): dim = self.corr_shells[icrsh][3] - for ibl, bn in enumerate(self.spin_block_names[self.corr_shells[icrsh][4]]): - isp = self.spin_names_to_ind[self.corr_shells[icrsh][4]][bn] + for isp, sp in enumerate(self.spin_block_names[self.corr_shells[icrsh][4]]): + isp = self.spin_names_to_ind[self.corr_shells[icrsh][4]][sp] for ik in range(self.n_k): 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*ibl) * self.h_field * MMat + 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] - self.Hsumk[icrsh][bn] += self.bz_weights[ik] * numpy.dot( numpy.dot(projmat,MMat), + self.Hsumk[icrsh][sp] += self.bz_weights[ik] * numpy.dot( numpy.dot(projmat,MMat), projmat.conjugate().transpose() ) # symmetrisation: @@ -514,42 +589,39 @@ class SumkLDA: # Rotate to local coordinate system: if self.use_rotations: for icrsh in range(self.n_corr_shells): - for bn in self.Hsumk[icrsh]: + for bl in self.Hsumk[icrsh]: - if self.rot_mat_time_inv[icrsh] == 1: self.Hsumk[icrsh][bn] = self.Hsumk[icrsh][bn].conjugate() - self.Hsumk[icrsh][bn] = numpy.dot( numpy.dot(self.rot_mat[icrsh].conjugate().transpose(),self.Hsumk[icrsh][bn]) , + 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]) , self.rot_mat[icrsh] ) # add to matrix: for ish in range(self.n_inequiv_shells): - for bn in eff_atlevels[ish]: - eff_atlevels[ish][bn] += self.Hsumk[self.inequiv_to_corr[ish]][bn] + for bl in eff_atlevels[ish]: + eff_atlevels[ish][bl] += self.Hsumk[self.inequiv_to_corr[ish]][bl] return eff_atlevels - def __init_dc(self): # construct the density matrix dm_imp and double counting arrays self.dc_imp = [ {} for icrsh in range(self.n_corr_shells)] for icrsh in range(self.n_corr_shells): dim = self.corr_shells[icrsh][3] - for j in range(len(self.gf_struct_sumk[icrsh])): - self.dc_imp[icrsh]['%s'%self.gf_struct_sumk[icrsh][j][0]] = numpy.zeros([dim,dim],numpy.float_) + spn = self.spin_block_names[self.corr_shells[icrsh][4]] + for sp in spn: self.dc_imp[icrsh][sp] = numpy.zeros([dim,dim],numpy.float_) self.dc_energ = [0.0 for icrsh in range(self.n_corr_shells)] - def set_dc(self,dens_mat,U_interact,J_hund,orb=0,use_dc_formula=0,use_val=None): """Sets the double counting term for inequiv orbital orb: - use_dc_formula=0: LDA+U FLL double counting, - use_dc_formula=1: Held's formula, - use_dc_formula=2: AMF. + use_dc_formula = 0: LDA+U FLL double counting, + use_dc_formula = 1: Held's formula, + use_dc_formula = 2: AMF. Be sure that you are using the correct interaction Hamiltonian!""" - for icrsh in range(self.n_corr_shells): iorb = self.corr_to_inequiv[icrsh] # iorb is the index of the inequivalent shell corresponding to icrsh @@ -559,110 +631,70 @@ class SumkLDA: Ncr = {} dim = self.corr_shells[icrsh][3] #*(1+self.corr_shells[icrsh][4]) - for j in range(len(self.gf_struct_sumk[icrsh])): - self.dc_imp[icrsh]['%s'%self.gf_struct_sumk[icrsh][j][0]] = numpy.identity(dim,numpy.float_) - blname = self.gf_struct_sumk[icrsh][j][0] - Ncr[blname] = 0.0 + spn = self.spin_block_names[self.corr_shells[icrsh][4]] + for sp in spn: + self.dc_imp[icrsh][sp] = numpy.identity(dim,numpy.float_) + Ncr[sp] = 0.0 for block,inner in self.gf_struct_solver[iorb].iteritems(): bl = self.solver_to_sumk_block[iorb][block] Ncr[bl] += dens_mat[block].real.trace() Ncrtot = 0.0 - block_ind_list = [block for block,inner in self.gf_struct_sumk[icrsh]] - for bl in block_ind_list: - Ncrtot += Ncr[bl] + + spn = self.spin_block_names[self.corr_shells[icrsh][4]] + for sp in spn: + Ncrtot += Ncr[sp] # average the densities if there is no SP: if self.SP == 0: - for bl in block_ind_list: - Ncr[bl] = Ncrtot / len(block_ind_list) + for sp in spn: Ncr[sp] = Ncrtot / len(spn) # correction for SO: we have only one block in this case, but in DC we need N/2 elif self.SP == 1 and self.SO == 1: - for bl in block_ind_list: - Ncr[bl] = Ncrtot / 2.0 + for sp in spn: Ncr[sp] = Ncrtot / 2.0 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 block_ind_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()) + for sp in spn: + Uav = U_interact*(Ncrtot-0.5) - J_hund*(Ncr[sp] - 0.5) + self.dc_imp[icrsh][sp] *= Uav + self.dc_energ[icrsh] -= J_hund / 2.0 * (Ncr[sp]) * (Ncr[sp]-1.0) + mpi.report("DC for shell %(icrsh)i and block %(sp)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 + (dim-1)*(U_interact-2.0*J_hund) + (dim-1)*(U_interact-3.0*J_hund))/(2*dim-1) / 2.0 * Ncrtot * (Ncrtot-1.0) - for bl in block_ind_list: + for sp in spn: Uav =(U_interact + (dim-1)*(U_interact-2.0*J_hund) + (dim-1)*(U_interact-3.0*J_hund))/(2*dim-1) * (Ncrtot-0.5) - self.dc_imp[icrsh][bl] *= Uav - mpi.report("DC for shell %(icrsh)i and block %(bl)s = %(Uav)f"%locals()) + self.dc_imp[icrsh][sp] *= Uav + mpi.report("DC for shell %(icrsh)i and block %(sp)s = %(Uav)f"%locals()) elif use_dc_formula == 2: # AMF self.dc_energ[icrsh] = 0.5 * U_interact * Ncrtot * Ncrtot - for bl in block_ind_list: - Uav = U_interact*(Ncrtot - Ncr[bl]/dim) - J_hund * (Ncr[bl] - Ncr[bl]/dim) - self.dc_imp[icrsh][bl] *= Uav - self.dc_energ[icrsh] -= (U_interact + (dim-1)*J_hund)/dim * 0.5 * Ncr[bl] * Ncr[bl] - mpi.report("DC for shell %(icrsh)i and block %(bl)s = %(Uav)f"%locals()) + for sp in spn: + Uav = U_interact*(Ncrtot - Ncr[sp]/dim) - J_hund * (Ncr[sp] - Ncr[sp]/dim) + self.dc_imp[icrsh][sp] *= Uav + self.dc_energ[icrsh] -= (U_interact + (dim-1)*J_hund)/dim * 0.5 * Ncr[sp] * Ncr[sp] + mpi.report("DC for shell %(icrsh)i and block %(sp)s = %(Uav)f"%locals()) # output: mpi.report("DC energy for shell %s = %s"%(icrsh,self.dc_energ[icrsh])) else: - block_ind_list = [block for block,inner in self.gf_struct_sumk[icrsh]] - for bl in block_ind_list: - self.dc_imp[icrsh][bl] *= use_val - self.dc_energ[icrsh] = use_val * Ncrtot + for sp in spn: + self.dc_imp[icrsh][sp] *= use_val # output: mpi.report("DC for shell %(icrsh)i = %(use_val)f"%locals()) mpi.report("DC energy = %s"%self.dc_energ[icrsh]) - - 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_shells, "give exactly one Sigma for each inequivalent corr. shell!" - - # init self.Sigma_imp: - 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) ] - 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) ] - else: - raise ValueError, "This type of Sigma is not handled." - - # transform the CTQMC blocks to the full matrix: - for icrsh in range(self.n_corr_shells): - ish = self.corr_to_inequiv[icrsh] # ish is the index of the inequivalent shell corresponding to icrsh - - 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)] - self.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') - - - def add_dc(self): """Substracts the double counting term from the impurity self energy.""" @@ -677,11 +709,16 @@ class SumkLDA: return sres # list of self energies corrected by DC + def symm_deg_gf(self,gf_to_symm,orb): + """Symmetrises a GF for the given degenerate shells self.deg_shells""" - def set_mu(self,mu): - """Sets a new chemical potential""" - - self.chemical_potential = mu + for degsh in self.deg_shells[orb]: + #loop over degenerate shells: + ss = gf_to_symm[degsh[0]].copy() + ss.zero() + n_deg = len(degsh) + for bl in degsh: ss += gf_to_symm[bl] / (1.0*n_deg) + for bl in degsh: gf_to_symm[bl] << ss def total_density(self, mu): @@ -694,11 +731,10 @@ class SumkLDA: """ dens = 0.0 - ikarray=numpy.array(range(self.n_k)) - + ikarray = numpy.array(range(self.n_k)) for ik in mpi.slice_array(ikarray): - S = self.lattice_gf_matsubara(ik=ik,mu=mu) + S = self.lattice_gf_matsubara(ik = ik, mu = mu) dens += self.bz_weights[ik] * S.total_density() # collect data from mpi: @@ -708,6 +744,12 @@ class SumkLDA: return dens + def set_mu(self,mu): + """Sets a new chemical potential""" + + self.chemical_potential = mu + + def find_mu(self, precision = 0.01): """ Searches for mu in order to give the desired charge @@ -715,7 +757,6 @@ class SumkLDA: """ F = lambda mu : self.total_density(mu = mu) - density = self.density_required - self.charge_below self.chemical_potential = dichotomy.dichotomy(function = F, @@ -725,62 +766,7 @@ class SumkLDA: verbosity = 3)[0] return self.chemical_potential - - - def extract_G_loc(self, mu=None, with_Sigma = True): - """ - Extracts the local downfolded Green function at the chemical potential of the class. - At the end, the local G is rotated from the global coordinate system to the local system. - if with_Sigma = False: Sigma is not included => non-interacting local GF - """ - - 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 - 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] - - 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) - Gloc[icrsh] += tmp - - #collect data from mpi: - for icrsh in range(self.n_corr_shells): - Gloc[icrsh] << mpi.all_reduce(mpi.world, Gloc[icrsh], lambda x,y : x+y) - mpi.barrier() - - # Gloc[:] is now the sum over k projected to the local orbitals. - # here comes the symmetrisation, if needed: - if self.symm_op != 0: Gloc = self.symmcorr.symmetrize(Gloc) - - # Gloc is rotated to the local coordinate system: - if self.use_rotations: - for icrsh in range(self.n_corr_shells): - 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() ], - 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] - - # return only the inequivalent shells: - return Glocret - + def calc_density_correction(self,filename = 'dens_mat.dat'): """ Calculates the density correction in order to feed it back to the DFT calculations.""" @@ -788,21 +774,17 @@ class SumkLDA: assert type(filename) == StringType, "filename has to be a string!" ntoi = self.spin_names_to_ind[self.SO] - bln = self.spin_block_names[self.SO] + spn = self.spin_block_names[self.SO] + dens = {sp: 0.0 for sp in spn} # Set up deltaN: deltaN = {} - for b in bln: - deltaN[b] = [ numpy.zeros( [self.n_orbitals[ik,ntoi[b]],self.n_orbitals[ik,ntoi[b]]], numpy.complex_) for ik in range(self.n_k)] - - ikarray=numpy.array(range(self.n_k)) - - dens = {} - for b in bln: - dens[b] = 0.0 + for sp in spn: + deltaN[sp] = [numpy.zeros([self.n_orbitals[ik,ntoi[sp]],self.n_orbitals[ik,ntoi[sp]]], numpy.complex_) for ik in range(self.n_k)] + 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) + 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() @@ -817,10 +799,10 @@ class SumkLDA: # now save to file: if mpi.is_master_node(): if self.SP == 0: - f=open(filename,'w') + f = open(filename,'w') else: - f=open(filename+'up','w') - f1=open(filename+'dn','w') + f = open(filename+'up','w') + f1 = open(filename+'dn','w') # 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)) @@ -847,7 +829,7 @@ class SumkLDA: if self.SO == 0: to_write = {f: (0, 'up'), f1: (1, 'down')} if self.SO == 1: to_write = {f: (0, 'ud'), f1: (0, 'ud')} for fout in to_write.iterkeys(): - isp, bn = to_write[fout] + isp, sp = to_write[fout] for ik in range(self.n_k): fout.write("%s\n"%self.n_orbitals[ik,isp]) for inu in range(self.n_orbitals[ik,isp]): @@ -859,7 +841,6 @@ class SumkLDA: return deltaN, dens - ################ # FIXME LEAVE UNDOCUMENTED ################ @@ -868,7 +849,7 @@ class SumkLDA: def find_mu_nonint(self, dens_req, orb = None, precision = 0.01): def F(mu): - gnonint = self.extract_G_loc(mu=mu,with_Sigma=False) + gnonint = self.extract_G_loc(mu = mu, with_Sigma = False) if orb is None: dens = 0.0 @@ -897,9 +878,9 @@ class SumkLDA: mu = self.chemical_potential def F(dc): - self.set_dc(dens_mat=dens_mat,U_interact=0,J_hund=0,orb=orb,use_val=dc) + self.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) + return self.total_density(mu = mu) else: return self.extract_G_loc()[orb].total_density() @@ -911,8 +892,8 @@ class SumkLDA: dcnew = dichotomy.dichotomy(function = F, x_init = guess, y_value = density, - precision_on_y = precision, delta_x=0.5, - max_loops = 100, x_name="Double-Counting", y_name= "Total 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 diff --git a/python/sumk_lda_tools.py b/python/sumk_lda_tools.py index b2e837a8..0816c2ea 100644 --- a/python/sumk_lda_tools.py +++ b/python/sumk_lda_tools.py @@ -86,19 +86,19 @@ class SumkLDATools(SumkLDA): 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] - bln = self.spin_block_names[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!" + 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[b]]) for b in bln ] - gf_struct = [ (bln[ibl], block_structure[ibl]) for ibl in range(self.n_spin_blocks[self.SO]) ] + 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] @@ -108,12 +108,12 @@ class SumkLDATools(SumkLDA): self.G_upfold_refreq.zero() GFsize = [ gf.N1 for bname,gf in self.G_upfold_refreq] - unchangedsize = all( [ self.n_orbitals[ik,ntoi[bln[ibl]]] == GFsize[ibl] - for ibl in range(self.n_spin_blocks[self.SO]) ] ) + 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[b]]) for b in bln ] - gf_struct = [ (bln[ibl], block_structure[ibl]) for ibl in range(self.n_spin_blocks[self.SO]) ] + 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] @@ -122,14 +122,14 @@ class SumkLDATools(SumkLDA): 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[b]],numpy.complex_) for b in bln] + 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 ibl in range(self.n_spin_blocks[self.SO]): - ind = ntoi[bln[ibl]] + for isp in range(self.n_spin_blocks[self.SO]): + ind = ntoi[spn[isp]] 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)) + 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: @@ -166,9 +166,9 @@ class SumkLDATools(SumkLDA): # init: Gloc = [] for icrsh in range(self.n_corr_shells): - b_list = [block for block,inner in self.gf_struct_sumk[icrsh]] + spn = self.spin_block_names[self.corr_shells[icrsh][4]] glist = lambda : [ GfReFreq(indices = inner, window = (om_min,om_max), n_points = n_om) for block,inner in self.gf_struct_sumk[icrsh]] - Gloc.append(BlockGf(name_list = b_list, block_list = glist(),make_copies=False)) + Gloc.append(BlockGf(name_list = spn, block_list = glist(),make_copies=False)) for icrsh in range(self.n_corr_shells): Gloc[icrsh].zero() # initialize to zero for ik in range(self.n_k): @@ -246,7 +246,7 @@ class SumkLDATools(SumkLDA): mu = self.chemical_potential - gf_struct_proj = [ [ (b, range(self.shells[i][3])) for b in self.spin_block_names[self.SO] ] for i in range(self.n_shells) ] + gf_struct_proj = [ [ (sp, range(self.shells[i][3])) for sp in self.spin_block_names[self.SO] ] for i in range(self.n_shells) ] Gproj = [BlockGf(name_block_generator = [ (block,GfReFreq(indices = inner, mesh = self.Sigma_imp[0].mesh)) for block,inner in gf_struct_proj[ish] ], make_copies = False ) for ish in range(self.n_shells)] for ish in range(self.n_shells): Gproj[ish].zero() @@ -366,7 +366,7 @@ class SumkLDATools(SumkLDA): # calculate A(k,w): mu = self.chemical_potential - bln = self.spin_block_names[self.SO] + spn = self.spin_block_names[self.SO] # init DOS: M = [x.real for x in self.Sigma_imp[0].mesh] @@ -381,19 +381,19 @@ class SumkLDATools(SumkLDA): if ishell is None: Akw = {} - for b in bln: Akw[b] = numpy.zeros([self.n_k, n_om ],numpy.float_) + for sp in spn: Akw[sp] = numpy.zeros([self.n_k, n_om ],numpy.float_) else: Akw = {} - for b in bln: Akw[b] = numpy.zeros([self.shells[ishell][3],self.n_k, n_om ],numpy.float_) + for sp in spn: Akw[sp] = numpy.zeros([self.shells[ishell][3],self.n_k, n_om ],numpy.float_) if fermi_surface: om_minplot = -2.0*broadening om_maxplot = 2.0*broadening Akw = {} - for b in bln: Akw[b] = numpy.zeros([self.n_k,1],numpy.float_) + for sp in spn: Akw[sp] = numpy.zeros([self.n_k,1],numpy.float_) - if not (ishell is None): - GFStruct_proj = [ (b, range(self.shells[ishell][3])) for b in bln ] + if not ishell is None: + GFStruct_proj = [ (sp, range(self.shells[ishell][3])) for sp in spn ] Gproj = BlockGf(name_block_generator = [ (block,GfReFreq(indices = inner, mesh = self.Sigma_imp[0].mesh)) for block,inner in GFStruct_proj ], make_copies = False) Gproj.zero() @@ -428,7 +428,7 @@ class SumkLDATools(SumkLDA): for iom in range(n_om): if (M[iom] > om_minplot) and (M[iom] < om_maxplot): for ish in range(self.shells[ishell][3]): - for ibn in bln: + for ibn in spn: Akw[ibn][ish,ik,iom] = Gproj[ibn].data[iom,ish,ish].imag/(-3.1415926535) @@ -436,7 +436,7 @@ class SumkLDATools(SumkLDA): if mpi.is_master_node(): if ishell is None: - for ibn in bln: + for ibn in spn: # loop over GF blocs: if invert_Akw: @@ -470,7 +470,7 @@ class SumkLDATools(SumkLDA): f.close() else: - for ibn in bln: + for ibn in spn: for ish in range(self.shells[ishell][3]): if invert_Akw: @@ -503,13 +503,13 @@ class SumkLDATools(SumkLDA): if self.symm_op: self.symmpar = Symmetry(self.hdf_file,subgroup=self.symmpar_data) # Density matrix in the window - bln = self.spin_block_names[self.SO] + spn = self.spin_block_names[self.SO] ntoi = self.spin_names_to_ind[self.SO] self.dens_mat_window = [ [numpy.zeros([self.shells[ish][3],self.shells[ish][3]],numpy.complex_) for ish in range(self.n_shells)] - for isp in range(len(bln)) ] # init the density matrix + for isp in range(len(spn)) ] # init the density matrix mu = self.chemical_potential - GFStruct_proj = [ [ (b, range(self.shells[i][3])) for b in bln ] for i in range(self.n_shells) ] + GFStruct_proj = [ [ (sp, range(self.shells[i][3])) for sp in spn ] for i in range(self.n_shells) ] if hasattr(self,"Sigma_imp"): Gproj = [BlockGf(name_block_generator = [ (block,GfImFreq(indices = inner, mesh = self.Sigma_imp[0].mesh)) for block,inner in GFStruct_proj[ish] ], make_copies = False) for ish in range(self.n_shells)] @@ -553,7 +553,7 @@ class SumkLDATools(SumkLDA): isp+=1 # add Density matrices to get the total: - dens_mat = [ [ self.dens_mat_below[ntoi[bln[isp]]][ish]+self.dens_mat_window[isp][ish] for ish in range(self.n_shells)] - for isp in range(len(bln)) ] + dens_mat = [ [ self.dens_mat_below[ntoi[spn[isp]]][ish]+self.dens_mat_window[isp][ish] for ish in range(self.n_shells)] + for isp in range(len(spn)) ] return dens_mat