LCPQ/dft_tools
3
0
Fork 0
mirror of https://github.com/triqs/dft_tools synced 2024-06-29 08:24:54 +02:00
Code Issues Releases Wiki Activity

Tidy up of indices

Browse Source
This commit is contained in:
Priyanka Seth 2014-11-15 20:04:54 +01:00
parent 518cbccd8c
commit 4cb0d67e02
2 changed files with 292 additions and 311 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

446
python/sumk_lda.py
View File

@ -39,7 +39,7 @@ class SumkLDA:
Initialises the class from data previously stored into an HDF5 Initialises the class from data previously stored into an HDF5
""" """
if not (type(hdf_file)==StringType): if not type(hdf_file) == StringType:
mpi.report("Give a string for the HDF5 filename to read the input!") mpi.report("Give a string for the HDF5 filename to read the input!")
else: else:
self.hdf_file = hdf_file self.hdf_file = hdf_file
@ -59,8 +59,8 @@ class SumkLDA:
'n_inequiv_shells', 'corr_to_inequiv', 'inequiv_to_corr'] 'n_inequiv_shells', 'corr_to_inequiv', 'inequiv_to_corr']
self.subgroup_present, self.value_read = self.read_input_from_hdf(subgrp = self.lda_data, things_to_read = things_to_read) self.subgroup_present, self.value_read = self.read_input_from_hdf(subgrp = self.lda_data, things_to_read = things_to_read)
if (self.SO) and (abs(self.h_field)>0.000001): if self.SO and (abs(self.h_field) > 0.000001):
self.h_field=0.0 self.h_field = 0.0
mpi.report("For SO, the external magnetic field is not implemented, setting it to 0!") mpi.report("For SO, the external magnetic field is not implemented, setting it to 0!")
self.spin_block_names = [ ['up','down'], ['ud'] ] self.spin_block_names = [ ['up','down'], ['ud'] ]
@ -73,8 +73,8 @@ class SumkLDA:
# 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[i][3])) for b in self.spin_block_names[self.corr_shells[i][4]] ] self.gf_struct_sumk = [ [ (b, range( self.corr_shells[icrsh][3])) for b in self.spin_block_names[self.corr_shells[icrsh][4]] ]
for i in xrange(self.n_corr_shells) ] for icrsh in range(self.n_corr_shells) ]
#----- #-----
# If these quantities are not in HDF, set them up # If these quantities are not in HDF, set them up
@ -83,9 +83,9 @@ 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[i]][3]) ) 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[i]][4]] ]) for b in self.spin_block_names[self.corr_shells[self.inequiv_to_corr[ish]][4]] ])
for i 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
self.sumk_to_solver = [ {} for ish in range(self.n_inequiv_shells) ] self.sumk_to_solver = [ {} for ish in range(self.n_inequiv_shells) ]
@ -105,18 +105,18 @@ class SumkLDA:
self.chemical_potential = mu self.chemical_potential = mu
if (not self.subgroup_present) or (not self.value_read['deg_shells']): if (not self.subgroup_present) or (not self.value_read['deg_shells']):
self.deg_shells = [ [] for i in range(self.n_inequiv_shells)] self.deg_shells = [ [] for ish in range(self.n_inequiv_shells)]
#----- #-----
if self.symm_op: if self.symm_op:
self.symmcorr = Symmetry(hdf_file,subgroup=self.symmcorr_data) self.symmcorr = Symmetry(hdf_file,subgroup=self.symmcorr_data)
# Analyse the block structure and determine the smallest blocks, if desired # Analyse the block structure and determine the smallest blocks, if desired
if (use_lda_blocks): dm=self.analyse_block_structure() if use_lda_blocks: dm = self.analyse_block_structure()
# Now save new things to HDF5: # Now save new things to HDF5:
# FIXME WHAT HAPPENS TO h_field? INPUT TO __INIT__? ADD TO OPTIONAL_THINGS? # FIXME WHAT HAPPENS TO h_field? INPUT TO __INIT__? ADD TO OPTIONAL_THINGS?
things_to_save=['chemical_potential','dc_imp','dc_energ','h_field'] things_to_save = ['chemical_potential','dc_imp','dc_energ','h_field']
self.save(things_to_save) self.save(things_to_save)
################ ################
@ -134,7 +134,7 @@ class SumkLDA:
for it in optional_things: setattr(self,it,0) for it in optional_things: setattr(self,it,0)
if mpi.is_master_node(): if mpi.is_master_node():
ar=HDFArchive(self.hdf_file,'a') ar = HDFArchive(self.hdf_file,'a')
if subgrp in ar: if subgrp in ar:
subgroup_present = True subgroup_present = True
# first read the necessary things: # first read the necessary things:
@ -145,7 +145,7 @@ class SumkLDA:
mpi.report("Loading %s failed!"%it) mpi.report("Loading %s failed!"%it)
value_read = False value_read = False
if (value_read and (len(optional_things)>0)): if value_read and (len(optional_things) > 0):
# if successfully read necessary items, read optional things: # if successfully read necessary items, read optional things:
value_read = {} value_read = {}
for it in optional_things: for it in optional_things:
@ -175,7 +175,7 @@ class SumkLDA:
if not (mpi.is_master_node()): return # do nothing on nodes if not (mpi.is_master_node()): return # do nothing on nodes
ar = HDFArchive(self.hdf_file,'a') ar = HDFArchive(self.hdf_file,'a')
if not (self.lda_output in ar): ar.create_group(self.lda_output) if not self.lda_output in ar: ar.create_group(self.lda_output)
for it in things_to_save: for it in things_to_save:
try: try:
ar[self.lda_output][it] = getattr(self,it) ar[self.lda_output][it] = getattr(self,it)
@ -219,21 +219,21 @@ class SumkLDA:
"""Local <-> Global rotation of a GF block. """Local <-> Global rotation of a GF block.
direction: 'toLocal' / 'toGlobal' """ direction: 'toLocal' / 'toGlobal' """
assert ((direction=='toLocal')or(direction=='toGlobal')),"Give direction 'toLocal' or 'toGlobal' in rotloc!" assert ((direction == 'toLocal') or (direction == 'toGlobal')),"Give direction 'toLocal' or 'toGlobal' in rotloc!"
gf_rotated = gf_to_rotate.copy() gf_rotated = gf_to_rotate.copy()
if (direction=='toGlobal'): if direction == 'toGlobal':
if ((self.rot_mat_time_inv[icrsh]==1) and (self.SO)): if (self.rot_mat_time_inv[icrsh] == 1) and self.SO:
gf_rotated << gf_rotated.transpose() gf_rotated << gf_rotated.transpose()
gf_rotated.from_L_G_R(self.rot_mat[icrsh].conjugate(),gf_rotated,self.rot_mat[icrsh].transpose()) gf_rotated.from_L_G_R(self.rot_mat[icrsh].conjugate(),gf_rotated,self.rot_mat[icrsh].transpose())
else: else:
gf_rotated.from_L_G_R(self.rot_mat[icrsh],gf_rotated,self.rot_mat[icrsh].conjugate().transpose()) gf_rotated.from_L_G_R(self.rot_mat[icrsh],gf_rotated,self.rot_mat[icrsh].conjugate().transpose())
elif (direction=='toLocal'): elif direction == 'toLocal':
if ((self.rot_mat_time_inv[icrsh]==1) and (self.SO)): if (self.rot_mat_time_inv[icrsh] == 1) and self.SO:
gf_rotated << gf_rotated.transpose() gf_rotated << gf_rotated.transpose()
gf_rotated.from_L_G_R(self.rot_mat[icrsh].transpose(),gf_rotated,self.rot_mat[icrsh].conjugate()) gf_rotated.from_L_G_R(self.rot_mat[icrsh].transpose(),gf_rotated,self.rot_mat[icrsh].conjugate())
else: else:
@ -249,9 +249,9 @@ class SumkLDA:
ntoi = self.spin_names_to_ind[self.SO] ntoi = self.spin_names_to_ind[self.SO]
bln = self.spin_block_names[self.SO] bln = 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:
stmp = self.add_dc() stmp = self.add_dc()
beta = self.Sigma_imp[0].mesh.beta # override beta if Sigma is present beta = self.Sigma_imp[0].mesh.beta # override beta if Sigma is present
@ -261,20 +261,20 @@ 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[bln[ibl]]] == GFsize[ibl]
for ibl in range(self.n_spin_blocks[self.SO]) ] ) for ibl 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[b]]) for b in bln ]
gf_struct = [ (bln[ibl], block_structure[ibl]) for ibl in range(self.n_spin_blocks[self.SO]) ] gf_struct = [ (bln[ibl], block_structure[ibl]) for ibl 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]
else: else:
glist = lambda : [ GfImFreq(indices = inner, beta = beta) for block,inner in gf_struct] 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.zero()
self.G_upfold << iOmega_n self.G_upfold << iOmega_n
@ -286,8 +286,8 @@ class SumkLDA:
M[ibl] = self.hopping[ik,ind,0:n_orb,0:n_orb] - (idmat[ibl]*mu) - (idmat[ibl] * self.h_field * (1-2*ibl)) M[ibl] = self.hopping[ik,ind,0:n_orb,0:n_orb] - (idmat[ibl]*mu) - (idmat[ibl] * self.h_field * (1-2*ibl))
self.G_upfold -= M self.G_upfold -= M
if (with_Sigma): if with_Sigma:
for icrsh in xrange(self.n_corr_shells): 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 self.G_upfold: gf -= self.upfold(ik,icrsh,bname,stmp[icrsh][bname],gf)
self.G_upfold.invert() self.G_upfold.invert()
@ -295,40 +295,36 @@ class SumkLDA:
return self.G_upfold return self.G_upfold
def simple_point_dens_mat(self): def simple_point_dens_mat(self):
ntoi = self.spin_names_to_ind[self.SO] ntoi = self.spin_names_to_ind[self.SO]
bln = self.spin_block_names[self.SO] bln = self.spin_block_names[self.SO]
MMat = [numpy.zeros( [self.n_orbitals[0,ntoi[bl]],self.n_orbitals[0,ntoi[bl]]], numpy.complex_) for bl in bln] MMat = [numpy.zeros( [self.n_orbitals[0,ntoi[bl]],self.n_orbitals[0,ntoi[bl]]], numpy.complex_) for bl in bln]
dens_mat = [ {} for icrsh in xrange(self.n_corr_shells)] dens_mat = [ {} for icrsh in range(self.n_corr_shells)]
for icrsh in xrange(self.n_corr_shells): for icrsh in range(self.n_corr_shells):
for bl in self.spin_block_names[self.corr_shells[icrsh][4]]: 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_) dens_mat[icrsh][bl] = numpy.zeros([self.corr_shells[icrsh][3],self.corr_shells[icrsh][3]], numpy.complex_)
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):
unchangedsize = all( [ self.n_orbitals[ik,ntoi[bln[ibl]]]==len(MMat[ibl]) unchangedsize = all( [ self.n_orbitals[ik,ntoi[bln[ibl]]] == len(MMat[ibl])
for ibl in range(self.n_spin_blocks[self.SO]) ] ) for ibl in range(self.n_spin_blocks[self.SO]) ] )
if (not unchangedsize): if not unchangedsize:
MMat = [numpy.zeros( [self.n_orbitals[ik,ntoi[bl]],self.n_orbitals[ik,ntoi[bl]]], numpy.complex_) for bl in bln] MMat = [numpy.zeros( [self.n_orbitals[ik,ntoi[bl]],self.n_orbitals[ik,ntoi[bl]]], numpy.complex_) for bl in bln]
for ibl, bl in enumerate(bln): for ibl, bl in enumerate(bln):
ind = ntoi[bl] ind = ntoi[bl]
for inu in range(self.n_orbitals[ik,ind]): 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) 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 MMat[ibl][inu,inu] = 1.0
else: else:
MMat[ibl][inu,inu] = 0.0 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]]): 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] isp = self.spin_names_to_ind[self.corr_shells[icrsh][4]][bn]
@ -341,31 +337,31 @@ class SumkLDA:
# get data from nodes: # get data from nodes:
for icrsh in range(self.n_corr_shells): for icrsh in range(self.n_corr_shells):
for bname in dens_mat[icrsh]: for bname in dens_mat[icrsh]:
dens_mat[icrsh][bname] = mpi.all_reduce(mpi.world,dens_mat[icrsh][bname],lambda x,y : x+y) 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) if self.symm_op != 0: dens_mat = self.symmcorr.symmetrize(dens_mat)
# Rotate to local coordinate system: # Rotate to local coordinate system:
if (self.use_rotations): if self.use_rotations:
for icrsh in xrange(self.n_corr_shells): for icrsh in range(self.n_corr_shells):
for bn in dens_mat[icrsh]: for bn in dens_mat[icrsh]:
if (self.rot_mat_time_inv[icrsh]==1): dens_mat[icrsh][bn] = dens_mat[icrsh][bn].conjugate() 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]) , dens_mat[icrsh][bn] = numpy.dot( numpy.dot(self.rot_mat[icrsh].conjugate().transpose(),dens_mat[icrsh][bn]) ,
self.rot_mat[icrsh]) self.rot_mat[icrsh] )
return dens_mat return dens_mat
# calculate upfolded gf, then density matrix -- no assumptions on structure (ie diagonal or not) # Calculate upfolded gf, then density matrix. No assumption on the structure made here (ie diagonal or not).
def density_gf(self,beta): def density_gf(self,beta):
"""Calculates the density without setting up Gloc. It is useful for Hubbard I, and very fast.""" """Calculates the density without setting up Gloc. It is useful for Hubbard I, and very quick."""
dens_mat = [ {} for icrsh in xrange(self.n_corr_shells)] dens_mat = [ {} for icrsh in range(self.n_corr_shells)]
for icrsh in xrange(self.n_corr_shells): for icrsh in range(self.n_corr_shells):
for bl in self.spin_block_names[self.corr_shells[icrsh][4]]: 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_) dens_mat[icrsh][bl] = numpy.zeros([self.corr_shells[icrsh][3],self.corr_shells[icrsh][3]], numpy.complex_)
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):
@ -386,16 +382,16 @@ class SumkLDA:
# get data from nodes: # get data from nodes:
for icrsh in range(self.n_corr_shells): for icrsh in range(self.n_corr_shells):
for bname in dens_mat[icrsh]: for bname in dens_mat[icrsh]:
dens_mat[icrsh][bname] = mpi.all_reduce(mpi.world,dens_mat[icrsh][bname],lambda x,y : x+y) 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) if self.symm_op != 0: dens_mat = self.symmcorr.symmetrize(dens_mat)
# Rotate to local coordinate system: # Rotate to local coordinate system:
if (self.use_rotations): if self.use_rotations:
for icrsh in xrange(self.n_corr_shells): for icrsh in range(self.n_corr_shells):
for bn in dens_mat[icrsh]: for bn in dens_mat[icrsh]:
if (self.rot_mat_time_inv[icrsh]==1): dens_mat[icrsh][bn] = dens_mat[icrsh][bn].conjugate() 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]) , dens_mat[icrsh][bn] = numpy.dot( numpy.dot(self.rot_mat[icrsh].conjugate().transpose(),dens_mat[icrsh][bn]) ,
self.rot_mat[icrsh] ) self.rot_mat[icrsh] )
@ -406,53 +402,49 @@ class SumkLDA:
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 function block structure from simple point integration."""
if dm is None: dm = self.simple_point_dens_mat() self.gf_struct_solver = [ {} for ish in range(self.n_inequiv_shells) ]
dens_mat = [dm[self.inequiv_to_corr[ish]] for ish in xrange(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) ]
self.solver_to_sumk = [ {} for ish in range(self.n_inequiv_shells) ] self.solver_to_sumk = [ {} for ish in range(self.n_inequiv_shells) ]
self.solver_to_sumk_block = [ {} for ish in range(self.n_inequiv_shells) ] self.solver_to_sumk_block = [ {} for ish in range(self.n_inequiv_shells) ]
if include_shells is None: include_shells=range(self.n_inequiv_shells) if dm is None: dm = self.simple_point_dens_mat()
dens_mat = [ dm[self.inequiv_to_corr[ish]] for ish in 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:
self.gf_struct_solver[ish] = {} block_ind_list = [ block for block,inner in self.gf_struct_sumk[self.inequiv_to_corr[ish]] ]
gf_struct_temp = []
block_ind_list = [block for block,inner in self.gf_struct_sumk[self.inequiv_to_corr[ish]] ]
for block in block_ind_list: for block in block_ind_list:
dm = dens_mat[ish][block] dm = dens_mat[ish][block]
dmbool = (abs(dm) > 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
offdiag = [] offdiag = []
for i in xrange(len(dmbool)): for i in range(len(dmbool)):
for j in xrange(i,len(dmbool)): for j in range(i+1,len(dmbool)):
if ((dmbool[i,j])&(i!=j)): offdiag.append([i,j]) if dmbool[i,j]: offdiag.append([i,j])
NBlocs = len(dmbool) num_blocs = len(dmbool)
blocs = [ [i] for i in range(NBlocs) ] blocs = [ [i] for i in range(num_blocs) ]
for i in range(len(offdiag)): for i in range(len(offdiag)):
if (offdiag[i][0]!=offdiag[i][1]): 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]] for j in range(i+1,len(offdiag)):
for j in range(i+1,len(offdiag)): if offdiag[j][0] == offdiag[i][1]: offdiag[j][0] = offdiag[i][0]
if (offdiag[j][0]==offdiag[i][1]): offdiag[j][0]=offdiag[i][0] if offdiag[j][1] == offdiag[i][1]: offdiag[j][1] = offdiag[i][0]
if (offdiag[j][1]==offdiag[i][1]): offdiag[j][1]=offdiag[i][0] if offdiag[j][0] > offdiag[i][1]: offdiag[j][0] -= 1
if (offdiag[j][0]>offdiag[i][1]): offdiag[j][0] -= 1 if offdiag[j][1] > offdiag[i][1]: offdiag[j][1] -= 1
if (offdiag[j][1]>offdiag[i][1]): offdiag[j][1] -= 1 offdiag[j].sort()
offdiag[j].sort() num_blocs -= 1
NBlocs-=1
for i in range(NBlocs): 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'%(block,i),range(len(blocs[i])))] )
gf_struct_temp.append( ('%s_%s'%(block,i),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(NBlocs): 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 = block
inner_sumk = blocs[i][j] inner_sumk = blocs[i][j]
@ -464,31 +456,31 @@ class SumkLDA:
# now calculate degeneracies of orbitals: # now calculate degeneracies of orbitals:
dm = {} dm = {}
for bl in gf_struct_temp: for block,inner in self.gf_struct_solver[ish].iteritems():
bln = bl[0]
ind = bl[1]
# get dm for the blocks: # get dm for the blocks:
dm[bln] = numpy.zeros([len(ind),len(ind)],numpy.complex_) dm[block] = numpy.zeros([len(inner),len(inner)],numpy.complex_)
for i in range(len(ind)): for ind1 in inner:
for j in range(len(ind)): for ind2 in inner:
dm[bln][i,j] = dens_mat[ish][self.solver_to_sumk_block[ish][bln]][ind[i],ind[j]] block_sumk,ind1_sumk = self.solver_to_sumk[ish][(block,ind1)]
block_sumk,ind2_sumk = self.solver_to_sumk[ish][(block,ind2)]
dm[block][ind1,ind2] = dens_mat[ish][block_sumk][ind1_sumk,ind2_sumk]
for bl in gf_struct_temp: for block1 in self.gf_struct_solver[ish].iterkeys():
for bl2 in gf_struct_temp: for block2 in self.gf_struct_solver[ish].iterkeys():
if (dm[bl[0]].shape==dm[bl2[0]].shape) : if dm[block1].shape == dm[block2].shape:
if ( ( (abs(dm[bl[0]]-dm[bl2[0]])<threshold).all() ) and (bl[0]!=bl2[0]) ): if ( (abs(dm[block1] - dm[block2]) < threshold).all() ) and (block1 != block2):
# check if it was already there: # check if it was already there:
ind1=-1 ind1 = -1
ind2=-2 ind2 = -2
for n,ind in enumerate(self.deg_shells[ish]): for n,ind in enumerate(self.deg_shells[ish]):
if (bl[0] in ind): ind1 = n if block1 in ind: ind1 = n
if (bl2[0] in ind): ind2 = n if block2 in ind: ind2 = n
if ((ind1 < 0) and (ind2 >= 0)): if (ind1 < 0) and (ind2 >= 0):
self.deg_shells[ish][ind2].append(bl[0]) self.deg_shells[ish][ind2].append(block1)
elif ((ind1 >= 0) and (ind2 < 0)): elif (ind1 >= 0) and (ind2 < 0):
self.deg_shells[ish][ind1].append(bl2[0]) self.deg_shells[ish][ind1].append(block2)
elif ((ind1 < 0) and (ind2 < 0)): elif (ind1 < 0) and (ind2 < 0):
self.deg_shells[ish].append([bl[0],bl2[0]]) self.deg_shells[ish].append([block1,block2])
things_to_save = ['gf_struct_solver','sumk_to_solver','solver_to_sumk','solver_to_sumk_block','deg_shells'] things_to_save = ['gf_struct_solver','sumk_to_solver','solver_to_sumk','solver_to_sumk_block','deg_shells']
self.save(things_to_save) self.save(things_to_save)
@ -518,19 +510,18 @@ class SumkLDA:
eff_atlevels[ish][bn] = numpy.identity(self.corr_shells[self.inequiv_to_corr[ish]][3], numpy.complex_) eff_atlevels[ish][bn] = numpy.identity(self.corr_shells[self.inequiv_to_corr[ish]][3], numpy.complex_)
# Chemical Potential: # Chemical Potential:
for ish in xrange(self.n_inequiv_shells): for ish in range(self.n_inequiv_shells):
for ii in eff_atlevels[ish]: eff_atlevels[ish][ii] *= -self.chemical_potential for ii in eff_atlevels[ish]: eff_atlevels[ish][ii] *= -self.chemical_potential
# double counting term: # double counting term:
#if hasattr(self,"dc_imp"): for ish in range(self.n_inequiv_shells):
for ish in xrange(self.n_inequiv_shells):
for ii in eff_atlevels[ish]: for ii in eff_atlevels[ish]:
eff_atlevels[ish][ii] -= self.dc_imp[self.inequiv_to_corr[ish]][ii] eff_atlevels[ish][ii] -= self.dc_imp[self.inequiv_to_corr[ish]][ii]
# sum over k: # sum over k:
if not hasattr(self,"Hsumk"): if not hasattr(self,"Hsumk"):
# 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 ish 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 bn 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])
@ -540,7 +531,7 @@ class SumkLDA:
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 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] isp = self.spin_names_to_ind[self.corr_shells[icrsh][4]][bn]
for ik in xrange(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*ibl) * self.h_field * MMat
@ -549,21 +540,19 @@ class SumkLDA:
projmat.conjugate().transpose() ) projmat.conjugate().transpose() )
# symmetrisation: # symmetrisation:
if (self.symm_op!=0): self.Hsumk = self.symmcorr.symmetrize(self.Hsumk) if self.symm_op != 0: self.Hsumk = self.symmcorr.symmetrize(self.Hsumk)
# Rotate to local coordinate system: # Rotate to local coordinate system:
if (self.use_rotations): if self.use_rotations:
for icrsh in xrange(self.n_corr_shells): for icrsh in range(self.n_corr_shells):
for bn in self.Hsumk[icrsh]: for bn 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][bn] = self.Hsumk[icrsh][bn].conjugate()
#if (self.corr_shells[icrsh][4]==0): 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]) , self.Hsumk[icrsh][bn] = numpy.dot( numpy.dot(self.rot_mat[icrsh].conjugate().transpose(),self.Hsumk[icrsh][bn]) ,
self.rot_mat[icrsh] ) self.rot_mat[icrsh] )
# add to matrix: # add to matrix:
for ish in xrange(self.n_inequiv_shells): for ish in range(self.n_inequiv_shells):
for bn in eff_atlevels[ish]: for bn in eff_atlevels[ish]:
eff_atlevels[ish][bn] += self.Hsumk[self.inequiv_to_corr[ish]][bn] eff_atlevels[ish][bn] += self.Hsumk[self.inequiv_to_corr[ish]][bn]
@ -575,13 +564,12 @@ class SumkLDA:
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.dm_imp = [ {} for i in xrange(self.n_corr_shells)] self.dc_imp = [ {} for icrsh in range(self.n_corr_shells)]
self.dc_imp = [ {} for i in xrange(self.n_corr_shells)] for icrsh in range(self.n_corr_shells):
for i in xrange(self.n_corr_shells): dim = self.corr_shells[icrsh][3]
l = self.corr_shells[i][3] for j in range(len(self.gf_struct_sumk[icrsh])):
for j in xrange(len(self.gf_struct_sumk[i])): self.dc_imp[icrsh]['%s'%self.gf_struct_sumk[icrsh][j][0]] = numpy.zeros([dim,dim],numpy.float_)
self.dc_imp[i]['%s'%self.gf_struct_sumk[i][j][0]] = numpy.zeros([l,l],numpy.float_) self.dc_energ = [0.0 for icrsh in range(self.n_corr_shells)]
self.dc_energ = [0.0 for i in xrange(self.n_corr_shells)]
@ -593,83 +581,80 @@ class SumkLDA:
Be sure that you are using the correct interaction Hamiltonian!""" Be sure that you are using the correct interaction Hamiltonian!"""
for icrsh in xrange(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
if (iorb==orb): if iorb != orb: continue # ignore this orbital
# do this orbital
Ncr = {}
l = self.corr_shells[icrsh][3] #*(1+self.corr_shells[icrsh][4])
for j in xrange(len(self.gf_struct_sumk[icrsh])): Ncr = {}
self.dc_imp[icrsh]['%s'%self.gf_struct_sumk[icrsh][j][0]] = numpy.identity(l,numpy.float_) dim = self.corr_shells[icrsh][3] #*(1+self.corr_shells[icrsh][4])
blname = self.gf_struct_sumk[icrsh][j][0]
Ncr[blname] = 0.0
for block,inner in self.gf_struct_solver[iorb].iteritems(): for j in range(len(self.gf_struct_sumk[icrsh])):
bl = self.solver_to_sumk_block[iorb][block] self.dc_imp[icrsh]['%s'%self.gf_struct_sumk[icrsh][j][0]] = numpy.identity(dim,numpy.float_)
Ncr[bl] += dens_mat[block].real.trace() blname = self.gf_struct_sumk[icrsh][j][0]
Ncr[blname] = 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()
M = self.corr_shells[icrsh][3] 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]
# 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)
# 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
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())
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:
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())
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())
# output:
mpi.report("DC energy for shell %s = %s"%(icrsh,self.dc_energ[icrsh]))
else:
Ncrtot = 0.0
block_ind_list = [block for block,inner in self.gf_struct_sumk[icrsh]] block_ind_list = [block for block,inner in self.gf_struct_sumk[icrsh]]
for bl in block_ind_list: for bl in block_ind_list:
Ncrtot += Ncr[bl] self.dc_imp[icrsh][bl] *= use_val
# average the densities if there is no SP: self.dc_energ[icrsh] = use_val * Ncrtot
if (self.SP==0):
for bl in block_ind_list:
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
elif (self.SP==1 and self.SO==1):
for bl in block_ind_list:
Ncr[bl] = Ncrtot / 2.0
if (use_val is None): # output:
mpi.report("DC for shell %(icrsh)i = %(use_val)f"%locals())
if (use_dc_formula==0): # FLL mpi.report("DC energy = %s"%self.dc_energ[icrsh])
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())
elif (use_dc_formula==1): # Held's formula, with U_interact the interorbital onsite interaction
self.dc_energ[icrsh] = (U_interact + (M-1)*(U_interact-2.0*J_hund) + (M-1)*(U_interact-3.0*J_hund))/(2*M-1) / 2.0 * Ncrtot * (Ncrtot-1.0)
for bl in block_ind_list:
Uav =(U_interact + (M-1)*(U_interact-2.0*J_hund) + (M-1)*(U_interact-3.0*J_hund))/(2*M-1) * (Ncrtot-0.5)
self.dc_imp[icrsh][bl] *= Uav
mpi.report("DC for shell %(icrsh)i and block %(bl)s = %(Uav)f"%locals())
elif (use_dc_formula==2): # AMF
self.dc_energ[icrsh] = 0.5 * U_interact * Ncrtot * Ncrtot
for bl in block_ind_list:
Uav = U_interact*(Ncrtot - Ncr[bl]/M) - J_hund * (Ncr[bl] - Ncr[bl]/M)
self.dc_imp[icrsh][bl] *= Uav
self.dc_energ[icrsh] -= (U_interact + (M-1)*J_hund)/M * 0.5 * Ncr[bl] * Ncr[bl]
mpi.report("DC for shell %(icrsh)i and block %(bl)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
# output:
mpi.report("DC for shell %(icrsh)i = %(use_val)f"%locals())
mpi.report("DC energy = %s"%self.dc_energ[icrsh])
@ -677,34 +662,34 @@ class SumkLDA:
"""Puts the impurity self energies for inequivalent atoms into the class, respects the multiplicity of the atoms.""" """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 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 len(Sigma_imp) == self.n_inequiv_shells, "give exactly one Sigma for each inequivalent corr. shell!"
# init self.Sigma_imp: # init self.Sigma_imp:
if all(type(gf) == GfImFreq for bname,gf in Sigma_imp[0]): if all(type(gf) == GfImFreq for bname,gf in Sigma_imp[0]):
# Imaginary frequency Sigma: # 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[i] ], 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 i in xrange(self.n_corr_shells) ] make_copies = False) for icrsh in range(self.n_corr_shells) ]
elif all(type(gf) == GfReFreq for bname,gf in Sigma_imp[0]): elif all(type(gf) == GfReFreq for bname,gf in Sigma_imp[0]):
# Real frequency Sigma: # 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[i] ], 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 i in xrange(self.n_corr_shells) ] make_copies = False) for icrsh in range(self.n_corr_shells) ]
else: else:
raise ValueError, "This type of Sigma is not handled." raise ValueError, "This type of Sigma is not handled."
# transform the CTQMC blocks to the full matrix: # transform the CTQMC blocks to the full matrix:
for icrsh in xrange(self.n_corr_shells): 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 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 block,inner in self.gf_struct_solver[ish].iteritems():
for ind1 in inner: for ind1 in inner:
for ind2 in inner: for ind2 in inner:
block_imp,ind1_imp = self.solver_to_sumk[ish][(block,ind1)] block_sumk,ind1_sumk = self.solver_to_sumk[ish][(block,ind1)]
block_imp,ind2_imp = self.solver_to_sumk[ish][(block,ind2)] block_sumk,ind2_sumk = self.solver_to_sumk[ish][(block,ind2)]
self.Sigma_imp[icrsh][block_imp][ind1_imp,ind2_imp] << Sigma_imp[ish][block][ind1,ind2] self.Sigma_imp[icrsh][block_sumk][ind1_sumk,ind2_sumk] << Sigma_imp[ish][block][ind1,ind2]
# rotation from local to global coordinate system: # rotation from local to global coordinate system:
if (self.use_rotations): if self.use_rotations:
for icrsh in xrange(self.n_corr_shells): 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 self.Sigma_imp[icrsh]: self.Sigma_imp[icrsh][bname] << self.rotloc(icrsh,gf,direction='toGlobal')
@ -714,7 +699,7 @@ class SumkLDA:
# Be careful: Sigma_imp is already in the global coordinate system!! # Be careful: Sigma_imp is already in the global coordinate system!!
sres = [s.copy() for s in self.Sigma_imp] sres = [s.copy() for s in self.Sigma_imp]
for icrsh in xrange(self.n_corr_shells): for icrsh in range(self.n_corr_shells):
for bname,gf in sres[icrsh]: for bname,gf in sres[icrsh]:
# Transform dc_imp to global coordinate system # 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())) dccont = numpy.dot(self.rot_mat[icrsh],numpy.dot(self.dc_imp[icrsh][bname],self.rot_mat[icrsh].conjugate().transpose()))
@ -732,12 +717,10 @@ class SumkLDA:
def total_density(self, mu): def total_density(self, mu):
""" """
Calculates the total charge for the energy window for a given mu. Since in general n_orbitals depends on k, Calculates the total charge for the energy window for a given chemical potential mu.
the calculation is done in the following order: Since in general n_orbitals depends on k, the calculation is done in the following order:
G_aa'(k,iw) -> n(k) = Tr G_aa'(k,iw) -> sum_k n_k G_aa'(k,iw) -> n(k) = Tr G_aa'(k,iw) -> sum_k n_k
mu: chemical potential
The calculation is done in the global coordinate system, if distinction is made between local/global! The calculation is done in the global coordinate system, if distinction is made between local/global!
""" """
@ -750,7 +733,7 @@ class SumkLDA:
dens += self.bz_weights[ik] * S.total_density() dens += self.bz_weights[ik] * S.total_density()
# collect data from mpi: # collect data from mpi:
dens = mpi.all_reduce(mpi.world,dens,lambda x,y : x+y) dens = mpi.all_reduce(mpi.world, dens, lambda x,y : x+y)
mpi.barrier() mpi.barrier()
return dens return dens
@ -766,7 +749,6 @@ class SumkLDA:
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,
x_init = self.chemical_potential, y_value = density, x_init = self.chemical_potential, y_value = density,
precision_on_y = precision, delta_x = 0.5, max_loops = 100, precision_on_y = precision, delta_x = 0.5, max_loops = 100,
@ -783,10 +765,10 @@ class SumkLDA:
if with_Sigma = False: Sigma is not included => non-interacting local GF if with_Sigma = False: Sigma is not included => non-interacting local GF
""" """
if (mu is None): mu = self.chemical_potential if mu is None: mu = self.chemical_potential
Gloc = [ self.Sigma_imp[icrsh].copy() for icrsh in xrange(self.n_corr_shells) ] # this list will be returned Gloc = [ self.Sigma_imp[icrsh].copy() for icrsh in range(self.n_corr_shells) ] # this list will be returned
for icrsh in xrange(self.n_corr_shells): Gloc[icrsh].zero() # initialize to zero for icrsh in range(self.n_corr_shells): Gloc[icrsh].zero() # initialize to zero
beta = Gloc[0].mesh.beta beta = Gloc[0].mesh.beta
ikarray=numpy.array(range(self.n_k)) ikarray=numpy.array(range(self.n_k))
@ -796,36 +778,36 @@ class SumkLDA:
S = self.lattice_gf_matsubara(ik=ik,mu=mu,with_Sigma = with_Sigma, beta = beta) S = self.lattice_gf_matsubara(ik=ik,mu=mu,with_Sigma = with_Sigma, beta = beta)
S *= self.bz_weights[ik] S *= self.bz_weights[ik]
for icrsh in xrange(self.n_corr_shells): for icrsh in range(self.n_corr_shells):
tmp = Gloc[icrsh].copy() # init temporary storage 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,S[bname],gf)
Gloc[icrsh] += tmp Gloc[icrsh] += tmp
#collect data from mpi: #collect data from mpi:
for icrsh in xrange(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) Gloc[icrsh] << mpi.all_reduce(mpi.world, Gloc[icrsh], lambda x,y : x+y)
mpi.barrier() mpi.barrier()
# Gloc[:] is now the sum over k projected to the local orbitals. # Gloc[:] is now the sum over k projected to the local orbitals.
# here comes the symmetrisation, if needed: # here comes the symmetrisation, if needed:
if (self.symm_op!=0): Gloc = self.symmcorr.symmetrize(Gloc) if self.symm_op != 0: Gloc = self.symmcorr.symmetrize(Gloc)
# Gloc is rotated to the local coordinate system: # Gloc is rotated to the local coordinate system:
if (self.use_rotations): if self.use_rotations:
for icrsh in xrange(self.n_corr_shells): for icrsh in range(self.n_corr_shells):
for bname,gf in Gloc[icrsh]: Gloc[icrsh][bname] << self.rotloc(icrsh,gf,direction='toLocal') for bname,gf in Gloc[icrsh]: Gloc[icrsh][bname] << self.rotloc(icrsh,gf,direction='toLocal')
# transform to CTQMC blocks: # 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() ], 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 xrange(self.n_inequiv_shells) ] make_copies = False) for ish in range(self.n_inequiv_shells) ]
for ish in xrange(self.n_inequiv_shells): for ish in range(self.n_inequiv_shells):
for block,inner in self.gf_struct_solver[ish].iteritems(): for block,inner in self.gf_struct_solver[ish].iteritems():
for ind1 in inner: for ind1 in inner:
for ind2 in inner: for ind2 in inner:
block_imp,ind1_imp = self.solver_to_sumk[ish][(block,ind1)] block_sumk,ind1_sumk = self.solver_to_sumk[ish][(block,ind1)]
block_imp,ind2_imp = self.solver_to_sumk[ish][(block,ind2)] block_sumk,ind2_sumk = self.solver_to_sumk[ish][(block,ind2)]
Glocret[ish][block][ind1,ind2] << Gloc[self.inequiv_to_corr[ish]][block_imp][ind1_imp,ind2_imp] Glocret[ish][block][ind1,ind2] << Gloc[self.inequiv_to_corr[ish]][block_sumk][ind1_sumk,ind2_sumk]
# return only the inequivalent shells: # return only the inequivalent shells:
return Glocret return Glocret
@ -834,7 +816,7 @@ class SumkLDA:
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] bln = self.spin_block_names[self.SO]
@ -859,25 +841,25 @@ class SumkLDA:
#put mpi Barrier: #put mpi Barrier:
for bname in deltaN: for bname in deltaN:
for ik in range(self.n_k): for ik in range(self.n_k):
deltaN[bname][ik] = mpi.all_reduce(mpi.world,deltaN[bname][ik],lambda x,y : x+y) deltaN[bname][ik] = mpi.all_reduce(mpi.world, deltaN[bname][ik], lambda x,y : x+y)
dens[bname] = mpi.all_reduce(mpi.world,dens[bname],lambda x,y : x+y) dens[bname] = mpi.all_reduce(mpi.world, dens[bname], lambda x,y : x+y)
mpi.barrier() mpi.barrier()
# now save to file: # now save to file:
if (mpi.is_master_node()): if mpi.is_master_node():
if (self.SP==0): if self.SP == 0:
f=open(filename,'w') f=open(filename,'w')
else: else:
f=open(filename+'up','w') f=open(filename+'up','w')
f1=open(filename+'dn','w') f1=open(filename+'dn','w')
# write chemical potential (in Rydberg): # write chemical potential (in Rydberg):
f.write("%.14f\n"%(self.chemical_potential/self.energy_unit)) f.write("%.14f\n"%(self.chemical_potential/self.energy_unit))
if (self.SP!=0): f1.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 # write beta in ryderg-1
f.write("%.14f\n"%(S.mesh.beta*self.energy_unit)) f.write("%.14f\n"%(S.mesh.beta*self.energy_unit))
if (self.SP!=0): f1.write("%.14f\n"%(S.mesh.beta*self.energy_unit)) if self.SP != 0: f1.write("%.14f\n"%(S.mesh.beta*self.energy_unit))
if (self.SP==0): # no spin-polarization if self.SP == 0: # no spin-polarization
for ik in range(self.n_k): for ik in range(self.n_k):
f.write("%s\n"%self.n_orbitals[ik,0]) f.write("%s\n"%self.n_orbitals[ik,0])
@ -890,11 +872,11 @@ class SumkLDA:
f.write("\n") f.write("\n")
f.close() f.close()
elif (self.SP==1): # with spin-polarization elif self.SP == 1: # with spin-polarization
# dict of filename : (spin index, block_name) # dict of filename: (spin index, block_name)
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, bn = to_write[fout]
for ik in range(self.n_k): for ik in range(self.n_k):
@ -919,7 +901,7 @@ class SumkLDA:
def F(mu): 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): if orb is None:
dens = 0.0 dens = 0.0
for ish in range(self.n_inequiv_shells): for ish in range(self.n_inequiv_shells):
dens += gnonint[ish].total_density() dens += gnonint[ish].total_density()
@ -947,13 +929,13 @@ class SumkLDA:
def F(dc): 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): if dens_req is None:
return self.total_density(mu=mu) return self.total_density(mu=mu)
else: else:
return self.extract_G_loc()[orb].total_density() return self.extract_G_loc()[orb].total_density()
if (dens_req is None): if dens_req is None:
density = self.density_required - self.charge_below density = self.density_required - self.charge_below
else: else:
density = dens_req density = dens_req

157
python/sumk_lda_tools.py
View File

@ -59,19 +59,19 @@ class SumkLDATools(SumkLDA):
"""Local <-> Global rotation of a GF block. """Local <-> Global rotation of a GF block.
direction: 'toLocal' / 'toGlobal' """ direction: 'toLocal' / 'toGlobal' """
assert ((direction=='toLocal')or(direction=='toGlobal')),"Give direction 'toLocal' or 'toGlobal' in rotloc!" assert (direction == 'toLocal' or direction == 'toGlobal'),"Give direction 'toLocal' or 'toGlobal' in rotloc!"
gf_rotated = gf_to_rotate.copy() gf_rotated = gf_to_rotate.copy()
if (direction=='toGlobal'): if direction == 'toGlobal':
if ((self.rot_mat_all_time_inv[ish]==1) and (self.SO)): if (self.rot_mat_all_time_inv[ish] == 1) and self.SO:
gf_rotated << gf_rotated.transpose() gf_rotated << gf_rotated.transpose()
gf_rotated.from_L_G_R(self.rot_mat_all[ish].conjugate(),gf_rotated,self.rot_mat_all[ish].transpose()) gf_rotated.from_L_G_R(self.rot_mat_all[ish].conjugate(),gf_rotated,self.rot_mat_all[ish].transpose())
else: else:
gf_rotated.from_L_G_R(self.rot_mat_all[ish],gf_rotated,self.rot_mat_all[ish].conjugate().transpose()) gf_rotated.from_L_G_R(self.rot_mat_all[ish],gf_rotated,self.rot_mat_all[ish].conjugate().transpose())
elif (direction=='toLocal'): elif direction == 'toLocal':
if ((self.rot_mat_all_time_inv[ish]==1)and(self.SO)): if (self.rot_mat_all_time_inv[ish] == 1) and self.SO:
gf_rotated << gf_rotated.transpose() gf_rotated << gf_rotated.transpose()
gf_rotated.from_L_G_R(self.rot_mat_all[ish].transpose(),gf_rotated,self.rot_mat_all[ish].conjugate()) gf_rotated.from_L_G_R(self.rot_mat_all[ish].transpose(),gf_rotated,self.rot_mat_all[ish].conjugate())
else: else:
@ -88,37 +88,37 @@ class SumkLDATools(SumkLDA):
ntoi = self.spin_names_to_ind[self.SO] ntoi = self.spin_names_to_ind[self.SO]
bln = self.spin_block_names[self.SO] bln = 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[b]]) for b in bln ]
gf_struct = [ (bln[ibl], block_structure[ibl]) for ibl in range(self.n_spin_blocks[self.SO]) ] gf_struct = [ (bln[ibl], block_structure[ibl]) for ibl 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]
else: else:
glist = lambda : [ GfReFreq(indices = inner, window=(mesh[0],mesh[1]),n_points=mesh[2]) for block,inner in gf_struct] 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 = BlockGf(name_list = block_ind_list, block_list = glist(),make_copies=False)
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[bln[ibl]]] == GFsize[ibl]
for ibl in range(self.n_spin_blocks[self.SO]) ] ) for ibl 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[b]]) for b in bln ]
gf_struct = [ (bln[ibl], block_structure[ibl]) for ibl in range(self.n_spin_blocks[self.SO]) ] gf_struct = [ (bln[ibl], block_structure[ibl]) for ibl 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]
else: else:
glist = lambda : [ GfReFreq(indices = inner, window=(mesh[0],mesh[1]),n_points=mesh[2]) for block,inner in gf_struct] 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 = BlockGf(name_list = block_ind_list, block_list = glist(),make_copies=False)
self.G_upfold_refreq.zero() self.G_upfold_refreq.zero()
@ -132,9 +132,9 @@ class SumkLDATools(SumkLDA):
M[ibl] = self.hopping[ik,ind,0:n_orb,0:n_orb] - (idmat[ibl]*mu) - (idmat[ibl] * self.h_field * (1-2*ibl)) M[ibl] = self.hopping[ik,ind,0:n_orb,0:n_orb] - (idmat[ibl]*mu) - (idmat[ibl] * self.h_field * (1-2*ibl))
self.G_upfold_refreq -= M self.G_upfold_refreq -= M
if (with_Sigma): if with_Sigma:
tmp = self.G_upfold_refreq.copy() # init temporary storage tmp = self.G_upfold_refreq.copy() # init temporary storage
for icrsh in xrange(self.n_corr_shells): for icrsh in range(self.n_corr_shells):
for bname,gf in tmp: tmp[bname] << self.upfold(ik,icrsh,bname,stmp[icrsh][bname],gf) 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 -= tmp # adding to the upfolded GF
@ -155,24 +155,23 @@ class SumkLDATools(SumkLDA):
for bn in self.spin_block_names[self.SO]: for bn in self.spin_block_names[self.SO]:
DOS[bn] = numpy.zeros([n_om],numpy.float_) DOS[bn] = numpy.zeros([n_om],numpy.float_)
DOSproj = [ {} for icrsh in range(self.n_inequiv_shells) ] DOSproj = [ {} for ish in range(self.n_inequiv_shells) ]
DOSproj_orb = [ {} for icrsh in range(self.n_inequiv_shells) ] DOSproj_orb = [ {} for ish in range(self.n_inequiv_shells) ]
for icrsh 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[icrsh]][4]]: for bn in self.spin_block_names[self.corr_shells[self.inequiv_to_corr[ish]][4]]:
dl = self.corr_shells[self.inequiv_to_corr[icrsh]][3] dim = self.corr_shells[self.inequiv_to_corr[ish]][3]
DOSproj[icrsh][bn] = numpy.zeros([n_om],numpy.float_) DOSproj[ish][bn] = numpy.zeros([n_om],numpy.float_)
DOSproj_orb[icrsh][bn] = numpy.zeros([n_om,dl,dl],numpy.float_) DOSproj_orb[ish][bn] = numpy.zeros([n_om,dim,dim],numpy.float_)
# 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]] b_list = [block for block,inner in self.gf_struct_sumk[icrsh]]
#glist = lambda : [ GfReFreq(indices = inner, beta = beta, mesh_array = mesh) 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]] 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 = b_list, block_list = glist(),make_copies=False))
for icrsh in xrange(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 xrange(self.n_k): 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.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_upfold *= self.bz_weights[ik]
@ -183,17 +182,17 @@ class SumkLDATools(SumkLDA):
asd = gf.data[iom,:,:].imag.trace()/(-3.1415926535) asd = gf.data[iom,:,:].imag.trace()/(-3.1415926535)
DOS[bname][iom] += asd DOS[bname][iom] += asd
for icrsh in xrange(self.n_corr_shells): for icrsh in range(self.n_corr_shells):
tmp = Gloc[icrsh].copy() 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_upfold[bname],gf) # downfolding G
Gloc[icrsh] += tmp Gloc[icrsh] += tmp
if (self.symm_op!=0): Gloc = self.symmcorr.symmetrize(Gloc) if self.symm_op != 0: Gloc = self.symmcorr.symmetrize(Gloc)
if (self.use_rotations): if self.use_rotations:
for icrsh in xrange(self.n_corr_shells): for icrsh in range(self.n_corr_shells):
for bname,gf in Gloc[icrsh]: Gloc[icrsh][bname] << self.rotloc(icrsh,gf,direction='toLocal') for bname,gf in Gloc[icrsh]: Gloc[icrsh][bname] << self.rotloc(icrsh,gf,direction='toLocal')
# Gloc can now also be used to look at orbitally resolved quantities # Gloc can now also be used to look at orbitally resolved quantities
@ -204,7 +203,7 @@ class SumkLDATools(SumkLDA):
DOSproj_orb[ish][bname][:,:,:] += gf.data[:,:,:].imag/(-3.1415926535) DOSproj_orb[ish][bname][:,:,:] += gf.data[:,:,:].imag/(-3.1415926535)
# output: # output:
if (mpi.is_master_node()): if mpi.is_master_node():
for bn in self.spin_block_names[self.SO]: for bn in self.spin_block_names[self.SO]:
f=open('DOS%s.dat'%bn, 'w') 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 i in range(n_om): f.write("%s %s\n"%(om_mesh[i],DOS[bn][i]))
@ -247,9 +246,9 @@ 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 xrange(self.n_shells) ] 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) ]
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 xrange(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()
Msh = [x.real for x in self.Sigma_imp[0].mesh] Msh = [x.real for x in self.Sigma_imp[0].mesh]
@ -263,9 +262,9 @@ class SumkLDATools(SumkLDA):
DOSproj_orb = [ {} for ish in range(self.n_shells) ] DOSproj_orb = [ {} for ish in range(self.n_shells) ]
for ish in range(self.n_shells): for ish in range(self.n_shells):
for bn in self.spin_block_names[self.SO]: for bn in self.spin_block_names[self.SO]:
dl = self.shells[ish][3] dim = self.shells[ish][3]
DOSproj[ish][bn] = numpy.zeros([n_om],numpy.float_) DOSproj[ish][bn] = numpy.zeros([n_om],numpy.float_)
DOSproj_orb[ish][bn] = numpy.zeros([n_om,dl,dl],numpy.float_) DOSproj_orb[ish][bn] = numpy.zeros([n_om,dim,dim],numpy.float_)
ikarray=numpy.array(range(self.n_k)) ikarray=numpy.array(range(self.n_k))
@ -279,24 +278,24 @@ class SumkLDATools(SumkLDA):
for bname,gf in S: DOS[bname][iom] += gf.data[iom,:,:].imag.trace()/(-3.1415926535) for bname,gf in S: DOS[bname][iom] += gf.data[iom,:,:].imag.trace()/(-3.1415926535)
#projected DOS: #projected DOS:
for ish in xrange(self.n_shells): for ish in range(self.n_shells):
tmp = Gproj[ish].copy() tmp = Gproj[ish].copy()
for ir in xrange(self.n_parproj[ish]): 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,S[bname],gf)
Gproj[ish] += tmp Gproj[ish] += tmp
# collect data from mpi: # collect data from mpi:
for bname in DOS: for bname in DOS:
DOS[bname] = mpi.all_reduce(mpi.world,DOS[bname],lambda x,y : x+y) DOS[bname] = mpi.all_reduce(mpi.world, DOS[bname], lambda x,y : x+y)
for ish in xrange(self.n_shells): for ish in range(self.n_shells):
Gproj[ish] << mpi.all_reduce(mpi.world,Gproj[ish],lambda x,y : x+y) Gproj[ish] << mpi.all_reduce(mpi.world, Gproj[ish], lambda x,y : x+y)
mpi.barrier() mpi.barrier()
if (self.symm_op!=0): Gproj = self.symmpar.symmetrize(Gproj) if self.symm_op != 0: Gproj = self.symmpar.symmetrize(Gproj)
# rotation to local coord. system: # rotation to local coord. system:
if (self.use_rotations): if self.use_rotations:
for ish in xrange(self.n_shells): for ish in range(self.n_shells):
for bname,gf in Gproj[ish]: Gproj[ish][bname] << self.rotloc_all(ish,gf,direction='toLocal') for bname,gf in Gproj[ish]: Gproj[ish][bname] << self.rotloc_all(ish,gf,direction='toLocal')
for ish in range(self.n_shells): for ish in range(self.n_shells):
@ -305,7 +304,7 @@ class SumkLDATools(SumkLDA):
DOSproj_orb[ish][bname][:,:,:] += gf.data[:,:,:].imag / (-3.1415926535) DOSproj_orb[ish][bname][:,:,:] += gf.data[:,:,:].imag / (-3.1415926535)
if (mpi.is_master_node()): if mpi.is_master_node():
# output to files # output to files
for bn in self.spin_block_names[self.SO]: for bn in self.spin_block_names[self.SO]:
f=open('./DOScorr%s.dat'%bn, 'w') f=open('./DOScorr%s.dat'%bn, 'w')
@ -341,14 +340,14 @@ class SumkLDATools(SumkLDA):
# FIXME CAN REMOVE? # FIXME CAN REMOVE?
# print hamiltonian for checks: # print hamiltonian for checks:
if ((self.SP==1)and(self.SO==0)): if self.SP == 1 and self.SO == 0:
f1=open('hamup.dat','w') f1=open('hamup.dat','w')
f2=open('hamdn.dat','w') f2=open('hamdn.dat','w')
for ik in xrange(self.n_k): for ik in range(self.n_k):
for i in xrange(self.n_orbitals[ik,0]): for i in range(self.n_orbitals[ik,0]):
f1.write('%s %s\n'%(ik,self.hopping[ik,0,i,i].real)) f1.write('%s %s\n'%(ik,self.hopping[ik,0,i,i].real))
for i in xrange(self.n_orbitals[ik,1]): for i in range(self.n_orbitals[ik,1]):
f2.write('%s %s\n'%(ik,self.hopping[ik,1,i,i].real)) f2.write('%s %s\n'%(ik,self.hopping[ik,1,i,i].real))
f1.write('\n') f1.write('\n')
f2.write('\n') f2.write('\n')
@ -356,8 +355,8 @@ class SumkLDATools(SumkLDA):
f2.close() f2.close()
else: else:
f=open('ham.dat','w') f=open('ham.dat','w')
for ik in xrange(self.n_k): for ik in range(self.n_k):
for i in xrange(self.n_orbitals[ik,0]): for i in range(self.n_orbitals[ik,0]):
f.write('%s %s\n'%(ik,self.hopping[ik,0,i,i].real)) f.write('%s %s\n'%(ik,self.hopping[ik,0,i,i].real))
f.write('\n') f.write('\n')
f.close() f.close()
@ -380,7 +379,7 @@ class SumkLDATools(SumkLDA):
om_minplot = plot_range[0] om_minplot = plot_range[0]
om_maxplot = plot_range[1] om_maxplot = plot_range[1]
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 b in bln: Akw[b] = numpy.zeros([self.n_k, n_om ],numpy.float_)
else: else:
@ -398,13 +397,13 @@ class SumkLDATools(SumkLDA):
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()
for ik in xrange(self.n_k): for ik in range(self.n_k):
S = self.lattice_gf_realfreq(ik=ik,mu=mu,broadening=broadening) S = self.lattice_gf_realfreq(ik=ik,mu=mu,broadening=broadening)
if (ishell is None): if ishell is None:
# non-projected A(k,w) # non-projected A(k,w)
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):
if fermi_surface: 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 S: Akw[bname][ik,0] += gf.data[iom,:,:].imag.trace()/(-3.1415926535) * (M[1]-M[0])
else: else:
@ -416,7 +415,7 @@ class SumkLDATools(SumkLDA):
# projected A(k,w): # projected A(k,w):
Gproj.zero() Gproj.zero()
tmp = Gproj.copy() tmp = Gproj.copy()
for ir in xrange(self.n_parproj[ishell]): 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,S[bname],gf)
Gproj += tmp Gproj += tmp
@ -427,20 +426,20 @@ class SumkLDATools(SumkLDA):
# for bname,gf in Gproj: Gproj[bname] << self.rotloc(0,gf,direction='toLocal') # for bname,gf in Gproj: Gproj[bname] << self.rotloc(0,gf,direction='toLocal')
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 bln:
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)
# END k-LOOP # END k-LOOP
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 bln:
# loop over GF blocs: # loop over GF blocs:
if (invert_Akw): if invert_Akw:
maxAkw=Akw[ibn].max() maxAkw=Akw[ibn].max()
minAkw=Akw[ibn].min() minAkw=Akw[ibn].min()
@ -453,15 +452,15 @@ class SumkLDATools(SumkLDA):
for ik in range(self.n_k): for ik in range(self.n_k):
if fermi_surface: if fermi_surface:
if (invert_Akw): if invert_Akw:
Akw[ibn][ik,0] = 1.0/(minAkw-maxAkw)*(Akw[ibn][ik,0] - maxAkw) Akw[ibn][ik,0] = 1.0/(minAkw-maxAkw)*(Akw[ibn][ik,0] - maxAkw)
f.write('%s %s\n'%(ik,Akw[ibn][ik,0])) f.write('%s %s\n'%(ik,Akw[ibn][ik,0]))
else: else:
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):
if (invert_Akw): if invert_Akw:
Akw[ibn][ik,iom] = 1.0/(minAkw-maxAkw)*(Akw[ibn][ik,iom] - maxAkw) Akw[ibn][ik,iom] = 1.0/(minAkw-maxAkw)*(Akw[ibn][ik,iom] - maxAkw)
if (shift>0.0001): if shift > 0.0001:
f.write('%s %s\n'%(M[iom],Akw[ibn][ik,iom])) f.write('%s %s\n'%(M[iom],Akw[ibn][ik,iom]))
else: else:
f.write('%s %s %s\n'%(ik,M[iom],Akw[ibn][ik,iom])) f.write('%s %s %s\n'%(ik,M[iom],Akw[ibn][ik,iom]))
@ -474,7 +473,7 @@ class SumkLDATools(SumkLDA):
for ibn in bln: for ibn in bln:
for ish in range(self.shells[ishell][3]): for ish in range(self.shells[ishell][3]):
if (invert_Akw): if invert_Akw:
maxAkw=Akw[ibn][ish,:,:].max() maxAkw=Akw[ibn][ish,:,:].max()
minAkw=Akw[ibn][ish,:,:].min() minAkw=Akw[ibn][ish,:,:].min()
@ -482,10 +481,10 @@ class SumkLDATools(SumkLDA):
for ik in range(self.n_k): for ik in range(self.n_k):
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):
if (invert_Akw): if invert_Akw:
Akw[ibn][ish,ik,iom] = 1.0/(minAkw-maxAkw)*(Akw[ibn][ish,ik,iom] - maxAkw) Akw[ibn][ish,ik,iom] = 1.0/(minAkw-maxAkw)*(Akw[ibn][ish,ik,iom] - maxAkw)
if (shift>0.0001): if shift > 0.0001:
f.write('%s %s\n'%(M[iom],Akw[ibn][ish,ik,iom])) f.write('%s %s\n'%(M[iom],Akw[ibn][ish,ik,iom]))
else: else:
f.write('%s %s %s\n'%(ik,M[iom],Akw[ibn][ish,ik,iom])) f.write('%s %s %s\n'%(ik,M[iom],Akw[ibn][ish,ik,iom]))
@ -510,16 +509,16 @@ class SumkLDATools(SumkLDA):
for isp in range(len(bln)) ] # init the density matrix for isp in range(len(bln)) ] # 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 xrange(self.n_shells) ] GFStruct_proj = [ [ (b, range(self.shells[i][3])) for b in bln ] 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 xrange(self.n_shells)] for ish in range(self.n_shells)]
beta = self.Sigma_imp[0].mesh.beta beta = self.Sigma_imp[0].mesh.beta
else: else:
Gproj = [BlockGf(name_block_generator = [ (block,GfImFreq(indices = inner, beta = beta)) for block,inner in GFStruct_proj[ish] ], make_copies = False) Gproj = [BlockGf(name_block_generator = [ (block,GfImFreq(indices = inner, beta = beta)) for block,inner in GFStruct_proj[ish] ], make_copies = False)
for ish in xrange(self.n_shells)] for ish in range(self.n_shells)]
for ish in xrange(self.n_shells): Gproj[ish].zero() for ish in range(self.n_shells): Gproj[ish].zero()
ikarray=numpy.array(range(self.n_k)) ikarray=numpy.array(range(self.n_k))
@ -527,25 +526,25 @@ class SumkLDATools(SumkLDA):
S = self.lattice_gf_matsubara(ik=ik,mu=mu,beta=beta) S = self.lattice_gf_matsubara(ik=ik,mu=mu,beta=beta)
S *= self.bz_weights[ik] S *= self.bz_weights[ik]
for ish in xrange(self.n_shells): for ish in range(self.n_shells):
tmp = Gproj[ish].copy() tmp = Gproj[ish].copy()
for ir in xrange(self.n_parproj[ish]): 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,S[bname],gf)
Gproj[ish] += tmp Gproj[ish] += tmp
#collect data from mpi: #collect data from mpi:
for ish in xrange(self.n_shells): for ish in range(self.n_shells):
Gproj[ish] << mpi.all_reduce(mpi.world,Gproj[ish],lambda x,y : x+y) Gproj[ish] << mpi.all_reduce(mpi.world, Gproj[ish], lambda x,y : x+y)
mpi.barrier() mpi.barrier()
# Symmetrisation: # Symmetrisation:
if (self.symm_op!=0): Gproj = self.symmpar.symmetrize(Gproj) if self.symm_op != 0: Gproj = self.symmpar.symmetrize(Gproj)
for ish in xrange(self.n_shells): for ish in range(self.n_shells):
# Rotation to local: # Rotation to local:
if (self.use_rotations): if self.use_rotations:
for bname,gf in Gproj[ish]: Gproj[ish][bname] << self.rotloc_all(ish,gf,direction='toLocal') for bname,gf in Gproj[ish]: Gproj[ish][bname] << self.rotloc_all(ish,gf,direction='toLocal')
isp = 0 isp = 0