LCPQ/dft_tools
3
0
Fork 0
mirror of https://github.com/triqs/dft_tools synced 2025-01-03 01:55:56 +01:00
Code Issues Releases Wiki Activity

More minor changes

Browse Source
This commit is contained in:
Priyanka Seth 2014-11-17 10:48:04 +01:00
parent 5eaa27a946
commit 5ad2acfee6
2 changed files with 270 additions and 289 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

469
python/sumk_lda.py
View File

@ -68,12 +68,12 @@ class SumkLDA:
# convert spin_block_names to indices -- if spin polarized, differentiate up and down blocks # convert spin_block_names to indices -- if spin polarized, differentiate up and down blocks
self.spin_names_to_ind = [{}, {}] self.spin_names_to_ind = [{}, {}]
for iso in range(2): # SO = 0 or 1 for iso in range(2): # SO = 0 or 1
for ibl in range(self.n_spin_blocks[iso]): for isp in range(self.n_spin_blocks[iso]):
self.spin_names_to_ind[iso][self.spin_block_names[iso][ibl]] = ibl * self.SP 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 # GF structure used for the local things in the k sums
# Most general form allowing for all hybridisation, i.e. largest blocks possible # 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) ] for icrsh in range(self.n_corr_shells) ]
#----- #-----
@ -83,8 +83,8 @@ class SumkLDA:
optional_things = optional_things) optional_things = optional_things)
if (not self.subgroup_present) or (not self.value_read['gf_struct_solver']): 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: # 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]) ) self.gf_struct_solver = [ dict([ (sp, 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]] ]) for sp in self.spin_block_names[self.corr_shells[self.inequiv_to_corr[ish]][4]] ])
for ish in range(self.n_inequiv_shells) for ish in range(self.n_inequiv_shells)
] ]
# Set standard (identity) maps from gf_struct_sumk <-> gf_struct_solver # Set standard (identity) maps from gf_struct_sumk <-> gf_struct_solver
@ -247,7 +247,7 @@ class SumkLDA:
and chemical potential mu.""" and chemical potential mu."""
ntoi = self.spin_names_to_ind[self.SO] 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 not hasattr(self,"Sigma_imp"): with_Sigma = False
@ -261,14 +261,14 @@ class SumkLDA:
set_up_G_upfold = True set_up_G_upfold = True
else: # yes if inconsistencies present in existing G_upfold else: # yes if inconsistencies present in existing G_upfold
GFsize = [ gf.N1 for bname,gf in self.G_upfold] GFsize = [ gf.N1 for bname,gf in self.G_upfold]
unchangedsize = all( [ self.n_orbitals[ik,ntoi[bln[ibl]]] == GFsize[ibl] unchangedsize = all( [ self.n_orbitals[ik,ntoi[spn[isp]]] == GFsize[isp]
for ibl in range(self.n_spin_blocks[self.SO]) ] ) 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 if (not unchangedsize) or (self.G_upfold.mesh.beta != beta): set_up_G_upfold = True
# Set up G_upfold # Set up G_upfold
if set_up_G_upfold: if set_up_G_upfold:
block_structure = [ range(self.n_orbitals[ik,ntoi[b]]) for b in bln ] block_structure = [ range(self.n_orbitals[ik,ntoi[sp]]) for sp in spn ]
gf_struct = [ (bln[ibl], block_structure[ibl]) for ibl in range(self.n_spin_blocks[self.SO]) ] 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] block_ind_list = [block for block,inner in gf_struct]
if with_Sigma: if with_Sigma:
glist = lambda : [ GfImFreq(indices = inner, mesh = self.Sigma_imp[0].mesh) for block,inner in gf_struct] glist = lambda : [ GfImFreq(indices = inner, mesh = self.Sigma_imp[0].mesh) for block,inner in gf_struct]
@ -278,12 +278,12 @@ class SumkLDA:
self.G_upfold.zero() self.G_upfold.zero()
self.G_upfold << iOmega_n 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) M = copy.deepcopy(idmat)
for ibl in range(self.n_spin_blocks[self.SO]): for isp in range(self.n_spin_blocks[self.SO]):
ind = ntoi[bln[ibl]] ind = ntoi[spn[isp]]
n_orb = self.n_orbitals[ik,ind] 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 self.G_upfold -= M
if with_Sigma: if with_Sigma:
@ -295,81 +295,94 @@ class SumkLDA:
return self.G_upfold return self.G_upfold
def density_matrix(self, method = 'using_gf', beta=40.0): def put_Sigma(self, Sigma_imp):
"""Calculate density matrices in one of two ways: """Puts the impurity self energies for inequivalent atoms into the class, respects the multiplicity of the atoms."""
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. 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!"
No assumption on the hopping structure is made (ie diagonal or not). assert len(Sigma_imp) == self.n_inequiv_shells, "give exactly one Sigma for each inequivalent corr. shell!"
if 'using_point_integration': Only works for diagonal hopping matrix (true in wien2k).
""" # init self.Sigma_imp:
dens_mat = [ {} for icrsh in range(self.n_corr_shells)] 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 icrsh in range(self.n_corr_shells):
for bl in self.spin_block_names[self.corr_shells[icrsh][4]]: ish = self.corr_to_inequiv[icrsh] # ish is the index of the inequivalent shell corresponding to icrsh
dens_mat[icrsh][bl] = numpy.zeros([self.corr_shells[icrsh][3],self.corr_shells[icrsh][3]], numpy.complex_) 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)) ikarray = numpy.array(range(self.n_k))
for ik in mpi.slice_array(ikarray): for ik in mpi.slice_array(ikarray):
if method == "using_gf": S = self.lattice_gf_matsubara(ik = ik, mu = mu, with_Sigma = with_Sigma, beta = beta)
S *= self.bz_weights[ik]
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
for icrsh in range(self.n_corr_shells): for icrsh in range(self.n_corr_shells):
for ibl, bn in enumerate(self.spin_block_names[self.corr_shells[icrsh][4]]): tmp = Gloc[icrsh].copy() # init temporary storage
isp = self.spin_names_to_ind[self.corr_shells[icrsh][4]][bn] for bname,gf in tmp: tmp[bname] << self.downfold(ik,icrsh,bname,S[bname],gf)
dim = self.corr_shells[icrsh][3] Gloc[icrsh] += tmp
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() )
# get data from nodes: # collect data from mpi
for icrsh in range(self.n_corr_shells): for icrsh in range(self.n_corr_shells): Gloc[icrsh] << mpi.all_reduce(mpi.world, Gloc[icrsh], lambda x,y : x+y)
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() 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: if self.use_rotations:
for icrsh in range(self.n_corr_shells): for icrsh in range(self.n_corr_shells):
for bn in dens_mat[icrsh]: for bname,gf in Gloc[icrsh]: Gloc[icrsh][bname] << self.rotloc(icrsh,gf,direction = 'toLocal')
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] )
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): 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.gf_struct_solver = [ {} for ish in range(self.n_inequiv_shells) ]
self.sumk_to_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) if include_shells is None: include_shells = range(self.n_inequiv_shells)
for ish in include_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: for sp in self.spin_block_names[self.corr_shells[self.inequiv_to_corr[ish]][4]]:
dm = dens_mat[ish][block] dmbool = (abs(dens_mat[ish][sp]) > threshold) # gives an index list of entries larger that threshold
dmbool = (abs(dm) > threshold) # gives an index list of entries larger that threshold
# Determine off-diagonal entries in upper triangular part of density matrix # Determine off-diagonal entries in upper triangular part of density matrix
offdiag = [] offdiag = []
@ -393,9 +405,9 @@ class SumkLDA:
for j in range(i+1,len(dmbool)): for j in range(i+1,len(dmbool)):
if dmbool[i,j]: offdiag.append([i,j]) if dmbool[i,j]: offdiag.append([i,j])
# Determine the number of non-hybridising blocks in the gf
num_blocs = len(dmbool) num_blocs = len(dmbool)
blocs = [ [i] for i in range(num_blocs) ] blocs = [ [i] for i in range(num_blocs) ]
for i in range(len(offdiag)): 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]) for j in range(len(blocs[offdiag[i][1]])): blocs[offdiag[i][0]].append(blocs[offdiag[i][1]][j])
del blocs[offdiag[i][1]] del blocs[offdiag[i][1]]
@ -407,17 +419,18 @@ class SumkLDA:
offdiag[j].sort() offdiag[j].sort()
num_blocs -= 1 num_blocs -= 1
# Set the gf_struct for the solver accordingly
for i in range(num_blocs): for i in range(num_blocs):
blocs[i].sort() 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) # 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) # and solver_to_sumk taking (solver_block, solver_inner) --> (sumk_block, sumk_index)
for i in range(num_blocs): for i in range(num_blocs):
for j in range(len(blocs[i])): for j in range(len(blocs[i])):
block_sumk = block block_sumk = sp
inner_sumk = blocs[i][j] inner_sumk = blocs[i][j]
block_solv = '%s_%s'%(block,i) block_solv = '%s_%s'%(sp,i)
inner_solv = j inner_solv = j
self.sumk_to_solver[ish][(block_sumk,inner_sumk)] = (block_solv,inner_solv) 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) self.solver_to_sumk[ish][(block_solv,inner_solv)] = (block_sumk,inner_sumk)
@ -457,16 +470,78 @@ class SumkLDA:
return dens_mat return dens_mat
def symm_deg_gf(self,gf_to_symm,orb): def density_matrix(self, method = 'using_gf', beta = 40.0):
"""Symmetrises a GF for the given degenerate shells self.deg_shells""" """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. # For simple dft input, get crystal field splittings.
def eff_atomic_levels(self): def eff_atomic_levels(self):
@ -475,8 +550,8 @@ class SumkLDA:
# define matrices for inequivalent shells: # define matrices for inequivalent shells:
eff_atlevels = [ {} for ish in range(self.n_inequiv_shells) ] eff_atlevels = [ {} for ish in range(self.n_inequiv_shells) ]
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]]: for sp 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_) eff_atlevels[ish][sp] = numpy.identity(self.corr_shells[self.inequiv_to_corr[ish]][3], numpy.complex_)
# Chemical Potential: # Chemical Potential:
for ish in range(self.n_inequiv_shells): 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: # 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) ] self.Hsumk = [ {} for icrsh in range(self.n_corr_shells) ]
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]) 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): for icrsh in range(self.n_corr_shells):
dim = self.corr_shells[icrsh][3] dim = self.corr_shells[icrsh][3]
for ibl, bn in enumerate(self.spin_block_names[self.corr_shells[icrsh][4]]): 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]][bn] isp = self.spin_names_to_ind[self.corr_shells[icrsh][4]][sp]
for ik in range(self.n_k): for ik in range(self.n_k):
n_orb = self.n_orbitals[ik,isp] n_orb = self.n_orbitals[ik,isp]
MMat = numpy.identity(n_orb, numpy.complex_) 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] 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() ) projmat.conjugate().transpose() )
# symmetrisation: # symmetrisation:
@ -514,34 +589,32 @@ class SumkLDA:
# Rotate to local coordinate system: # Rotate to local coordinate system:
if self.use_rotations: if self.use_rotations:
for icrsh in range(self.n_corr_shells): 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() if self.rot_mat_time_inv[icrsh] == 1: self.Hsumk[icrsh][bl] = self.Hsumk[icrsh][bl].conjugate()
self.Hsumk[icrsh][bn] = numpy.dot( numpy.dot(self.rot_mat[icrsh].conjugate().transpose(),self.Hsumk[icrsh][bn]) , self.Hsumk[icrsh][bl] = numpy.dot( numpy.dot(self.rot_mat[icrsh].conjugate().transpose(),self.Hsumk[icrsh][bl]) ,
self.rot_mat[icrsh] ) self.rot_mat[icrsh] )
# add to matrix: # add to matrix:
for ish in range(self.n_inequiv_shells): for ish in range(self.n_inequiv_shells):
for bn in eff_atlevels[ish]: for bl in eff_atlevels[ish]:
eff_atlevels[ish][bn] += self.Hsumk[self.inequiv_to_corr[ish]][bn] eff_atlevels[ish][bl] += self.Hsumk[self.inequiv_to_corr[ish]][bl]
return eff_atlevels return eff_atlevels
def __init_dc(self): def __init_dc(self):
# construct the density matrix dm_imp and double counting arrays # construct the density matrix dm_imp and double counting arrays
self.dc_imp = [ {} for icrsh in range(self.n_corr_shells)] self.dc_imp = [ {} for icrsh in range(self.n_corr_shells)]
for icrsh in range(self.n_corr_shells): for icrsh in range(self.n_corr_shells):
dim = self.corr_shells[icrsh][3] dim = self.corr_shells[icrsh][3]
for j in range(len(self.gf_struct_sumk[icrsh])): spn = self.spin_block_names[self.corr_shells[icrsh][4]]
self.dc_imp[icrsh]['%s'%self.gf_struct_sumk[icrsh][j][0]] = numpy.zeros([dim,dim],numpy.float_) 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)] 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): 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: """Sets the double counting term for inequiv orbital orb:
use_dc_formula = 0: LDA+U FLL double counting, use_dc_formula = 0: LDA+U FLL double counting,
@ -549,7 +622,6 @@ class SumkLDA:
use_dc_formula = 2: AMF. use_dc_formula = 2: AMF.
Be sure that you are using the correct interaction Hamiltonian!""" Be sure that you are using the correct interaction Hamiltonian!"""
for icrsh in range(self.n_corr_shells): 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 iorb = self.corr_to_inequiv[icrsh] # iorb is the index of the inequivalent shell corresponding to icrsh
@ -559,110 +631,70 @@ class SumkLDA:
Ncr = {} Ncr = {}
dim = self.corr_shells[icrsh][3] #*(1+self.corr_shells[icrsh][4]) dim = self.corr_shells[icrsh][3] #*(1+self.corr_shells[icrsh][4])
for j in range(len(self.gf_struct_sumk[icrsh])): spn = self.spin_block_names[self.corr_shells[icrsh][4]]
self.dc_imp[icrsh]['%s'%self.gf_struct_sumk[icrsh][j][0]] = numpy.identity(dim,numpy.float_) for sp in spn:
blname = self.gf_struct_sumk[icrsh][j][0] self.dc_imp[icrsh][sp] = numpy.identity(dim,numpy.float_)
Ncr[blname] = 0.0 Ncr[sp] = 0.0
for block,inner in self.gf_struct_solver[iorb].iteritems(): for block,inner in self.gf_struct_solver[iorb].iteritems():
bl = self.solver_to_sumk_block[iorb][block] bl = self.solver_to_sumk_block[iorb][block]
Ncr[bl] += dens_mat[block].real.trace() Ncr[bl] += dens_mat[block].real.trace()
Ncrtot = 0.0 Ncrtot = 0.0
block_ind_list = [block for block,inner in self.gf_struct_sumk[icrsh]]
for bl in block_ind_list: spn = self.spin_block_names[self.corr_shells[icrsh][4]]
Ncrtot += Ncr[bl] for sp in spn:
Ncrtot += Ncr[sp]
# average the densities if there is no SP: # average the densities if there is no SP:
if self.SP == 0: if self.SP == 0:
for bl in block_ind_list: for sp in spn: Ncr[sp] = Ncrtot / len(spn)
Ncr[bl] = Ncrtot / len(block_ind_list)
# correction for SO: we have only one block in this case, but in DC we need N/2 # 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: elif self.SP == 1 and self.SO == 1:
for bl in block_ind_list: for sp in spn: Ncr[sp] = Ncrtot / 2.0
Ncr[bl] = Ncrtot / 2.0
if use_val is None: if use_val is None:
if use_dc_formula == 0: # FLL if use_dc_formula == 0: # FLL
self.dc_energ[icrsh] = U_interact / 2.0 * Ncrtot * (Ncrtot-1.0) self.dc_energ[icrsh] = U_interact / 2.0 * Ncrtot * (Ncrtot-1.0)
for bl in block_ind_list: for sp in spn:
Uav = U_interact*(Ncrtot-0.5) - J_hund*(Ncr[bl] - 0.5) Uav = U_interact*(Ncrtot-0.5) - J_hund*(Ncr[sp] - 0.5)
self.dc_imp[icrsh][bl] *= Uav self.dc_imp[icrsh][sp] *= Uav
self.dc_energ[icrsh] -= J_hund / 2.0 * (Ncr[bl]) * (Ncr[bl]-1.0) self.dc_energ[icrsh] -= J_hund / 2.0 * (Ncr[sp]) * (Ncr[sp]-1.0)
mpi.report("DC for shell %(icrsh)i and block %(bl)s = %(Uav)f"%locals()) 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 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) 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) 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 self.dc_imp[icrsh][sp] *= Uav
mpi.report("DC for shell %(icrsh)i and block %(bl)s = %(Uav)f"%locals()) mpi.report("DC for shell %(icrsh)i and block %(sp)s = %(Uav)f"%locals())
elif use_dc_formula == 2: # AMF elif use_dc_formula == 2: # AMF
self.dc_energ[icrsh] = 0.5 * U_interact * Ncrtot * Ncrtot self.dc_energ[icrsh] = 0.5 * U_interact * Ncrtot * Ncrtot
for bl in block_ind_list: for sp in spn:
Uav = U_interact*(Ncrtot - Ncr[bl]/dim) - J_hund * (Ncr[bl] - Ncr[bl]/dim) Uav = U_interact*(Ncrtot - Ncr[sp]/dim) - J_hund * (Ncr[sp] - Ncr[sp]/dim)
self.dc_imp[icrsh][bl] *= Uav self.dc_imp[icrsh][sp] *= Uav
self.dc_energ[icrsh] -= (U_interact + (dim-1)*J_hund)/dim * 0.5 * Ncr[bl] * Ncr[bl] 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 %(bl)s = %(Uav)f"%locals()) mpi.report("DC for shell %(icrsh)i and block %(sp)s = %(Uav)f"%locals())
# output: # output:
mpi.report("DC energy for shell %s = %s"%(icrsh,self.dc_energ[icrsh])) mpi.report("DC energy for shell %s = %s"%(icrsh,self.dc_energ[icrsh]))
else: 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 self.dc_energ[icrsh] = use_val * Ncrtot
for sp in spn:
self.dc_imp[icrsh][sp] *= use_val
# output: # output:
mpi.report("DC for shell %(icrsh)i = %(use_val)f"%locals()) mpi.report("DC for shell %(icrsh)i = %(use_val)f"%locals())
mpi.report("DC energy = %s"%self.dc_energ[icrsh]) 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): def add_dc(self):
"""Substracts the double counting term from the impurity self energy.""" """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 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): for degsh in self.deg_shells[orb]:
"""Sets a new chemical potential""" #loop over degenerate shells:
ss = gf_to_symm[degsh[0]].copy()
self.chemical_potential = mu 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): def total_density(self, mu):
@ -695,7 +732,6 @@ class SumkLDA:
dens = 0.0 dens = 0.0
ikarray = numpy.array(range(self.n_k)) ikarray = numpy.array(range(self.n_k))
for ik in mpi.slice_array(ikarray): 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)
@ -708,6 +744,12 @@ class SumkLDA:
return dens return dens
def set_mu(self,mu):
"""Sets a new chemical potential"""
self.chemical_potential = mu
def find_mu(self, precision = 0.01): def find_mu(self, precision = 0.01):
""" """
Searches for mu in order to give the desired charge Searches for mu in order to give the desired charge
@ -715,7 +757,6 @@ class SumkLDA:
""" """
F = lambda mu : self.total_density(mu = mu) F = lambda mu : self.total_density(mu = mu)
density = self.density_required - self.charge_below density = self.density_required - self.charge_below
self.chemical_potential = dichotomy.dichotomy(function = F, self.chemical_potential = dichotomy.dichotomy(function = F,
@ -727,80 +768,21 @@ class SumkLDA:
return self.chemical_potential 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'): def calc_density_correction(self,filename = 'dens_mat.dat'):
""" Calculates the density correction in order to feed it back to the DFT calculations.""" """ 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, "filename has to be a string!"
ntoi = self.spin_names_to_ind[self.SO] 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: # Set up deltaN:
deltaN = {} deltaN = {}
for b in bln: for sp in spn:
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)] 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)) ikarray = numpy.array(range(self.n_k))
dens = {}
for b in bln:
dens[b] = 0.0
for ik in mpi.slice_array(ikarray): 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: for bname,gf in S:
@ -847,7 +829,7 @@ class SumkLDA:
if self.SO == 0: to_write = {f: (0, 'up'), f1: (1, 'down')} if self.SO == 0: to_write = {f: (0, 'up'), f1: (1, 'down')}
if self.SO == 1: to_write = {f: (0, 'ud'), f1: (0, 'ud')} if self.SO == 1: to_write = {f: (0, 'ud'), f1: (0, 'ud')}
for fout in to_write.iterkeys(): for fout in to_write.iterkeys():
isp, bn = to_write[fout] isp, sp = to_write[fout]
for ik in range(self.n_k): for ik in range(self.n_k):
fout.write("%s\n"%self.n_orbitals[ik,isp]) fout.write("%s\n"%self.n_orbitals[ik,isp])
for inu in range(self.n_orbitals[ik,isp]): for inu in range(self.n_orbitals[ik,isp]):
@ -859,7 +841,6 @@ class SumkLDA:
return deltaN, dens return deltaN, dens
################ ################
# FIXME LEAVE UNDOCUMENTED # FIXME LEAVE UNDOCUMENTED
################ ################
@ -912,7 +893,7 @@ class SumkLDA:
dcnew = dichotomy.dichotomy(function = F, dcnew = dichotomy.dichotomy(function = F,
x_init = guess, y_value = density, x_init = guess, y_value = density,
precision_on_y = precision, delta_x = 0.5, precision_on_y = precision, delta_x = 0.5,
max_loops = 100, x_name="Double-Counting", y_name= "Total Density", max_loops = 100, x_name = "Double Counting", y_name= "Total Density",
verbosity = 3)[0] verbosity = 3)[0]
return dcnew return dcnew

58
python/sumk_lda_tools.py
View File

@ -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.""" 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] 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 not hasattr(self,"Sigma_imp"): with_Sigma=False
if with_Sigma: if with_Sigma:
assert all(type(gf) == GfReFreq for bname,gf in self.Sigma_imp[0]), "Real frequency Sigma needed for lattice_gf_realfreq!" 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() stmp = self.add_dc()
else: 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: if self.G_upfold_refreq is None:
# first setting up of G_upfold_refreq # first setting up of G_upfold_refreq
block_structure = [ range(self.n_orbitals[ik,ntoi[b]]) for b in bln ] block_structure = [ range(self.n_orbitals[ik,ntoi[sp]]) for sp in spn ]
gf_struct = [ (bln[ibl], block_structure[ibl]) for ibl in range(self.n_spin_blocks[self.SO]) ] 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] block_ind_list = [block for block,inner in gf_struct]
if with_Sigma: if with_Sigma:
glist = lambda : [ GfReFreq(indices = inner, mesh=self.Sigma_imp[0].mesh) for block,inner in gf_struct] 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() self.G_upfold_refreq.zero()
GFsize = [ gf.N1 for bname,gf in self.G_upfold_refreq] GFsize = [ gf.N1 for bname,gf in self.G_upfold_refreq]
unchangedsize = all( [ self.n_orbitals[ik,ntoi[bln[ibl]]] == GFsize[ibl] unchangedsize = all( [ self.n_orbitals[ik,ntoi[spn[isp]]] == GFsize[isp]
for ibl in range(self.n_spin_blocks[self.SO]) ] ) for isp in range(self.n_spin_blocks[self.SO]) ] )
if not unchangedsize: if not unchangedsize:
block_structure = [ range(self.n_orbitals[ik,ntoi[b]]) for b in bln ] block_structure = [ range(self.n_orbitals[ik,ntoi[sp]]) for sp in spn ]
gf_struct = [ (bln[ibl], block_structure[ibl]) for ibl in range(self.n_spin_blocks[self.SO]) ] 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] block_ind_list = [block for block,inner in gf_struct]
if with_Sigma: if with_Sigma:
glist = lambda : [ GfReFreq(indices = inner, mesh =self.Sigma_imp[0].mesh) for block,inner in gf_struct] 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 = BlockGf(name_list = block_ind_list, block_list = glist(),make_copies=False)
self.G_upfold_refreq.zero() 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 self.G_upfold_refreq << Omega + 1j*broadening
M = copy.deepcopy(idmat) M = copy.deepcopy(idmat)
for ibl in range(self.n_spin_blocks[self.SO]): for isp in range(self.n_spin_blocks[self.SO]):
ind = ntoi[bln[ibl]] ind = ntoi[spn[isp]]
n_orb = self.n_orbitals[ik,ind] 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 self.G_upfold_refreq -= M
if with_Sigma: if with_Sigma:
@ -166,9 +166,9 @@ class SumkLDATools(SumkLDA):
# init: # init:
Gloc = [] Gloc = []
for icrsh in range(self.n_corr_shells): 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]] 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 icrsh in range(self.n_corr_shells): Gloc[icrsh].zero() # initialize to zero
for ik in range(self.n_k): for ik in range(self.n_k):
@ -246,7 +246,7 @@ class SumkLDATools(SumkLDA):
mu = self.chemical_potential 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 ) 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)]
for ish in range(self.n_shells): Gproj[ish].zero() for ish in range(self.n_shells): Gproj[ish].zero()
@ -366,7 +366,7 @@ class SumkLDATools(SumkLDA):
# calculate A(k,w): # calculate A(k,w):
mu = self.chemical_potential mu = self.chemical_potential
bln = self.spin_block_names[self.SO] spn = self.spin_block_names[self.SO]
# init DOS: # init DOS:
M = [x.real for x in self.Sigma_imp[0].mesh] M = [x.real for x in self.Sigma_imp[0].mesh]
@ -381,19 +381,19 @@ class SumkLDATools(SumkLDA):
if ishell is None: if ishell is None:
Akw = {} 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: else:
Akw = {} 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: if fermi_surface:
om_minplot = -2.0*broadening om_minplot = -2.0*broadening
om_maxplot = 2.0*broadening om_maxplot = 2.0*broadening
Akw = {} 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): if not ishell is None:
GFStruct_proj = [ (b, range(self.shells[ishell][3])) for b in bln ] 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 = 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() Gproj.zero()
@ -428,7 +428,7 @@ class SumkLDATools(SumkLDA):
for iom in range(n_om): for iom in range(n_om):
if (M[iom] > om_minplot) and (M[iom] < om_maxplot): if (M[iom] > om_minplot) and (M[iom] < om_maxplot):
for ish in range(self.shells[ishell][3]): for ish in range(self.shells[ishell][3]):
for ibn in bln: for ibn in spn:
Akw[ibn][ish,ik,iom] = Gproj[ibn].data[iom,ish,ish].imag/(-3.1415926535) 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 mpi.is_master_node():
if ishell is None: if ishell is None:
for ibn in bln: for ibn in spn:
# loop over GF blocs: # loop over GF blocs:
if invert_Akw: if invert_Akw:
@ -470,7 +470,7 @@ class SumkLDATools(SumkLDA):
f.close() f.close()
else: else:
for ibn in bln: for ibn in spn:
for ish in range(self.shells[ishell][3]): for ish in range(self.shells[ishell][3]):
if invert_Akw: if invert_Akw:
@ -503,13 +503,13 @@ class SumkLDATools(SumkLDA):
if self.symm_op: self.symmpar = Symmetry(self.hdf_file,subgroup=self.symmpar_data) if self.symm_op: self.symmpar = Symmetry(self.hdf_file,subgroup=self.symmpar_data)
# Density matrix in the window # 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] ntoi = self.spin_names_to_ind[self.SO]
self.dens_mat_window = [ [numpy.zeros([self.shells[ish][3],self.shells[ish][3]],numpy.complex_) for ish in range(self.n_shells)] self.dens_mat_window = [ [numpy.zeros([self.shells[ish][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 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"): 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) 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)] for ish in range(self.n_shells)]
@ -553,7 +553,7 @@ class SumkLDATools(SumkLDA):
isp+=1 isp+=1
# add Density matrices to get the total: # 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)] 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(bln)) ] for isp in range(len(spn)) ]
return dens_mat return dens_mat