LCPQ/dft_tools
3
0
Fork 0
mirror of https://github.com/triqs/dft_tools synced 2024-12-22 04:13:47 +01:00
Code Issues Releases Wiki Activity

Variable name changes for clarity and consistency

Browse Source
This commit is contained in:
Priyanka Seth 2014-11-14 09:43:28 +01:00
parent 0114562baa
commit 6708788ed7
5 changed files with 272 additions and 270 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

9
python/TODOFIX
View File

@ -8,6 +8,15 @@ Substitutions:
* read_symmetry_input -> convert_symmetry_input * read_symmetry_input -> convert_symmetry_input
* Symm_corr -> symmcorr * Symm_corr -> symmcorr
internal substitutions:
Symm_par --> symmpar
sig -> bname
names_to_ind -> spin_names_to_ind
n_spin_blocks_gf -> n_spin_blocks
block_names -> spin_block_names
a_list -> block_ind_list
a,al -> block,inner
********** **********
* changed default h5 subgroup names * changed default h5 subgroup names

347
python/sumk_lda.py
View File

@ -49,8 +49,6 @@ class SumkLDA:
self.symmpar_data = symmpar_data self.symmpar_data = symmpar_data
self.bands_data = bands_data self.bands_data = bands_data
self.lda_output = lda_output self.lda_output = lda_output
self.block_names = [ ['up','down'], ['ud'] ]
self.n_spin_blocks_gf = [2,1]
self.G_upfold = None self.G_upfold = None
self.h_field = h_field self.h_field = h_field
@ -64,42 +62,44 @@ class SumkLDA:
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!")
# FIXME -- REMOVE THIS, WRITE DATA IN CONVERTER
# determine the number of inequivalent correlated shells (self.n_inequiv_corr_shells) # determine the number of inequivalent correlated shells (self.n_inequiv_corr_shells)
# and related maps (self.shellmap, self.invshellmap) # and related maps (self.shellmap, self.invshellmap)
self.inequiv_shells(self.corr_shells) self.inequiv_shells(self.corr_shells)
# field to convert block_names to indices self.spin_block_names = [ ['up','down'], ['ud'] ]
self.names_to_ind = [{}, {}] self.n_spin_blocks = [2,1]
for ibl in range(2): # convert spin_block_names to indices -- if spin polarized, differentiate up and down blocks
for inm in range(self.n_spin_blocks_gf[ibl]): self.spin_names_to_ind = [{}, {}]
self.names_to_ind[ibl][self.block_names[ibl][inm]] = inm * self.SP for iso in range(2): # SO = 0 or 1
for ibl in range(self.n_spin_blocks[iso]):
self.spin_names_to_ind[iso][self.spin_block_names[iso][ibl]] = ibl * self.SP
# 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_corr = [ [ (al, range( self.corr_shells[i][3])) for al in self.block_names[self.corr_shells[i][4]] ] self.gf_struct_corr = [ [ (b, range( self.corr_shells[i][3])) for b in self.spin_block_names[self.corr_shells[i][4]] ]
for i in xrange(self.n_corr_shells) ] for i in xrange(self.n_corr_shells) ]
#----- #-----
# If these quantities are not in HDF, set them up and save into HDF # If these quantities are not in HDF, set them up
optional_things = ['gf_struct_solver','map_inv','map','chemical_potential','dc_imp','dc_energ','deg_shells'] optional_things = ['gf_struct_solver','map_inv','map','chemical_potential','dc_imp','dc_energ','deg_shells']
self.subgroup_present, self.value_read = self.read_input_from_hdf(subgrp = self.lda_output, things_to_read = [], self.subgroup_present, self.value_read = self.read_input_from_hdf(subgrp = self.lda_output, things_to_read = [],
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([ (al, range(self.corr_shells[self.invshellmap[i]][3]) ) self.gf_struct_solver = [ dict([ (b, range(self.corr_shells[self.invshellmap[i]][3]) )
for al in self.block_names[self.corr_shells[self.invshellmap[i]][4]] ]) for b in self.spin_block_names[self.corr_shells[self.invshellmap[i]][4]] ])
for i in range(self.n_inequiv_corr_shells) for i in range(self.n_inequiv_corr_shells)
] ]
self.map = [ {} for i in xrange(self.n_inequiv_corr_shells) ] self.map = [ {} for i in xrange(self.n_inequiv_corr_shells) ]
self.map_inv = [ {} for i in xrange(self.n_inequiv_corr_shells) ] self.map_inv = [ {} for i in xrange(self.n_inequiv_corr_shells) ]
for i in xrange(self.n_inequiv_corr_shells): for i in xrange(self.n_inequiv_corr_shells):
for al in self.block_names[self.corr_shells[self.invshellmap[i]][4]]: for b in self.spin_block_names[self.corr_shells[self.invshellmap[i]][4]]:
self.map[i][al] = [al for j in range( self.corr_shells[self.invshellmap[i]][3] ) ] self.map[i][b] = [b for j in range( self.corr_shells[self.invshellmap[i]][3] ) ]
self.map_inv[i][al] = al self.map_inv[i][b] = b
if (not self.subgroup_present) or (not self.value_read['dc_imp']): if (not self.subgroup_present) or (not self.value_read['dc_imp']):
# init the double counting: self.__init_dc() # initialise the double counting
self.__init_dc()
if (not self.subgroup_present) or (not self.value_read['chemical_potential']): if (not self.subgroup_present) or (not self.value_read['chemical_potential']):
self.chemical_potential = mu self.chemical_potential = mu
@ -113,9 +113,9 @@ class SumkLDA:
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 blocs, if desired # Analyse the block structure and determine the smallest blocs, if desired
if (use_lda_blocks): dm=self.analyse_BS() if (use_lda_blocks): dm=self.analyse_block_structure()
# Now save things again 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)
@ -130,7 +130,7 @@ class SumkLDA:
""" """
value_read = True value_read = True
# init variables on all nodes to ensure mpi broadcast works at the end # initialise variables on all nodes to ensure mpi broadcast works at the end
for it in things_to_read: setattr(self,it,0) for it in things_to_read: setattr(self,it,0)
for it in optional_things: setattr(self,it,0) for it in optional_things: setattr(self,it,0)
@ -147,7 +147,7 @@ class SumkLDA:
value_read = False value_read = False
if (value_read and (len(optional_things)>0)): if (value_read and (len(optional_things)>0)):
# if necessary things worked, now 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:
if it in ar[subgrp]: if it in ar[subgrp]:
@ -172,7 +172,7 @@ class SumkLDA:
def save(self,things_to_save): def save(self,things_to_save):
"""Saves some quantities into an HDF5 arxiv""" """Saves some quantities into an HDF5 archive"""
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')
@ -188,11 +188,11 @@ class SumkLDA:
# CORE FUNCTIONS # CORE FUNCTIONS
################ ################
def downfold(self,ik,icrsh,sig,gf_to_downfold,gf_inp): def downfold(self,ik,icrsh,bname,gf_to_downfold,gf_inp):
"""Downfolding a block of the Greens function""" """Downfolding a block of the Greens function"""
gf_downfolded = gf_inp.copy() gf_downfolded = gf_inp.copy()
isp = self.names_to_ind[self.SO][sig] # get spin index for proj. matrices isp = self.spin_names_to_ind[self.SO][bname] # get spin index for proj. matrices
dim = self.corr_shells[icrsh][3] dim = self.corr_shells[icrsh][3]
n_orb = self.n_orbitals[ik,isp] n_orb = self.n_orbitals[ik,isp]
projmat = self.proj_mat[ik,isp,icrsh,0:dim,0:n_orb] projmat = self.proj_mat[ik,isp,icrsh,0:dim,0:n_orb]
@ -202,11 +202,11 @@ class SumkLDA:
return gf_downfolded return gf_downfolded
def upfold(self,ik,icrsh,sig,gf_to_upfold,gf_inp): def upfold(self,ik,icrsh,bname,gf_to_upfold,gf_inp):
"""Upfolding a block of the Greens function""" """Upfolding a block of the Greens function"""
gf_upfolded = gf_inp.copy() gf_upfolded = gf_inp.copy()
isp = self.names_to_ind[self.SO][sig] # get spin index for proj. matrices isp = self.spin_names_to_ind[self.SO][bname] # get spin index for proj. matrices
dim = self.corr_shells[icrsh][3] dim = self.corr_shells[icrsh][3]
n_orb = self.n_orbitals[ik,isp] n_orb = self.n_orbitals[ik,isp]
projmat = self.proj_mat[ik,isp,icrsh,0:dim,0:n_orb] projmat = self.proj_mat[ik,isp,icrsh,0:dim,0:n_orb]
@ -247,8 +247,8 @@ class SumkLDA:
"""Calculates the lattice Green function from the LDA hopping and the self energy at k-point number ik """Calculates the lattice Green function from the LDA hopping and the self energy at k-point number ik
and chemical potential mu.""" and chemical potential mu."""
ntoi = self.names_to_ind[self.SO] ntoi = self.spin_names_to_ind[self.SO]
bln = self.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
@ -261,27 +261,27 @@ class SumkLDA:
if self.G_upfold is None: # yes if not G_upfold provided if self.G_upfold is None: # yes if not G_upfold provided
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 sig,gf in self.G_upfold] GFsize = [ gf.N1 for bname,gf in self.G_upfold]
unchangedsize = all( [ self.n_orbitals[ik,ntoi[bln[ib]]]==GFsize[ib] unchangedsize = all( [ self.n_orbitals[ik,ntoi[bln[ibl]]]==GFsize[ibl]
for ib in range(self.n_spin_blocks_gf[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:
BS = [ range(self.n_orbitals[ik,ntoi[ib]]) for ib in bln ] block_structure = [ range(self.n_orbitals[ik,ntoi[b]]) for b in bln ]
gf_struct = [ (bln[ib], BS[ib]) for ib in range(self.n_spin_blocks_gf[self.SO]) ] gf_struct = [ (bln[ibl], block_structure[ibl]) for ibl in range(self.n_spin_blocks[self.SO]) ]
a_list = [a for a,al in gf_struct] block_ind_list = [block for block,inner in gf_struct]
if (with_Sigma): if (with_Sigma):
glist = lambda : [ GfImFreq(indices = al, mesh = self.Sigma_imp[0].mesh) for a,al 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 = al, beta = beta) for a,al in gf_struct] glist = lambda : [ GfImFreq(indices = inner, beta = beta) for block,inner in gf_struct]
self.G_upfold = BlockGf(name_list = a_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
idmat = [numpy.identity(self.n_orbitals[ik,ntoi[bl]],numpy.complex_) for bl in bln] idmat = [numpy.identity(self.n_orbitals[ik,ntoi[bl]],numpy.complex_) for bl in bln]
M = copy.deepcopy(idmat) M = copy.deepcopy(idmat)
for ibl in range(self.n_spin_blocks_gf[self.SO]): for ibl in range(self.n_spin_blocks[self.SO]):
ind = ntoi[bln[ibl]] ind = ntoi[bln[ibl]]
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[ibl] = self.hopping[ik,ind,0:n_orb,0:n_orb] - (idmat[ibl]*mu) - (idmat[ibl] * self.h_field * (1-2*ibl))
@ -289,7 +289,7 @@ class SumkLDA:
if (with_Sigma): if (with_Sigma):
for icrsh in xrange(self.n_corr_shells): for icrsh in xrange(self.n_corr_shells):
for sig,gf in self.G_upfold: gf -= self.upfold(ik,icrsh,sig,stmp[icrsh][sig],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()
@ -301,27 +301,27 @@ class SumkLDA:
def simple_point_dens_mat(self): def simple_point_dens_mat(self):
ntoi = self.names_to_ind[self.SO] ntoi = self.spin_names_to_ind[self.SO]
bln = self.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 xrange(self.n_corr_shells)]
for icrsh in xrange(self.n_corr_shells): for icrsh in xrange(self.n_corr_shells):
for bl in self.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[ib]]]==len(MMat[ib]) unchangedsize = all( [ self.n_orbitals[ik,ntoi[bln[ibl]]]==len(MMat[ibl])
for ib in range(self.n_spin_blocks_gf[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)
@ -331,22 +331,20 @@ class SumkLDA:
for icrsh in range(self.n_corr_shells): for icrsh in range(self.n_corr_shells):
for ibn,bn in enumerate(self.block_names[self.corr_shells[icrsh][4]]): for ibl, bn in enumerate(self.spin_block_names[self.corr_shells[icrsh][4]]):
isp = self.names_to_ind[self.corr_shells[icrsh][4]][bn] isp = self.spin_names_to_ind[self.corr_shells[icrsh][4]][bn]
dim = self.corr_shells[icrsh][3] dim = self.corr_shells[icrsh][3]
n_orb = self.n_orbitals[ik,isp] n_orb = self.n_orbitals[ik,isp]
projmat = self.proj_mat[ik,isp,icrsh,0:dim,0:n_orb]
#print ik, bn, isp dens_mat[icrsh][bn] += self.bz_weights[ik] * numpy.dot( numpy.dot(projmat,MMat[ibl]) ,
dens_mat[icrsh][bn] += self.bz_weights[ik] * numpy.dot( numpy.dot(self.proj_mat[ik,isp,icrsh,0:dim,0:n_orb],MMat[ibn]) , projmat.transpose().conjugate() )
self.proj_mat[ik,isp,icrsh,0:dim,0:n_orb].transpose().conjugate() )
# 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 sig in dens_mat[icrsh]: for bname in dens_mat[icrsh]:
dens_mat[icrsh][sig] = mpi.all_reduce(mpi.world,dens_mat[icrsh][sig],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:
@ -357,7 +355,6 @@ class SumkLDA:
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 assumptions on structure (ie diagonal or not)
@ -366,7 +363,7 @@ class SumkLDA:
dens_mat = [ {} for icrsh in xrange(self.n_corr_shells)] dens_mat = [ {} for icrsh in xrange(self.n_corr_shells)]
for icrsh in xrange(self.n_corr_shells): for icrsh in xrange(self.n_corr_shells):
for bl in self.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))
@ -376,24 +373,23 @@ class SumkLDA:
G_upfold = self.lattice_gf_matsubara(ik=ik, beta=beta, mu=self.chemical_potential) G_upfold = self.lattice_gf_matsubara(ik=ik, beta=beta, mu=self.chemical_potential)
G_upfold *= self.bz_weights[ik] G_upfold *= self.bz_weights[ik]
dm = G_upfold.density() dm = G_upfold.density()
MMat = [dm[bl] for bl in self.block_names[self.SO]] MMat = [dm[bl] for bl in self.spin_block_names[self.SO]]
for icrsh in range(self.n_corr_shells): for icrsh in range(self.n_corr_shells):
for ibn,bn in enumerate(self.block_names[self.corr_shells[icrsh][4]]): for ibl, bn in enumerate(self.spin_block_names[self.corr_shells[icrsh][4]]):
isp = self.names_to_ind[self.corr_shells[icrsh][4]][bn] isp = self.spin_names_to_ind[self.corr_shells[icrsh][4]][bn]
dim = self.corr_shells[icrsh][3] dim = self.corr_shells[icrsh][3]
n_orb = self.n_orbitals[ik,isp] n_orb = self.n_orbitals[ik,isp]
#print ik, bn, isp projmat = self.proj_mat[ik,isp,icrsh,0:dim,0:n_orb]
dens_mat[icrsh][bn] += numpy.dot( numpy.dot(self.proj_mat[ik,isp,icrsh,0:dim,0:n_orb],MMat[ibn]), dens_mat[icrsh][bn] += numpy.dot( numpy.dot(projmat,MMat[ibl]),
self.proj_mat[ik,isp,icrsh,0:dim,0:n_orb].transpose().conjugate() ) projmat.transpose().conjugate() )
# 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 sig in dens_mat[icrsh]: for bname in dens_mat[icrsh]:
dens_mat[icrsh][sig] = mpi.all_reduce(mpi.world,dens_mat[icrsh][sig],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:
@ -408,7 +404,7 @@ class SumkLDA:
def analyse_BS(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() if dm is None: dm = self.simple_point_dens_mat()
@ -421,9 +417,9 @@ class SumkLDA:
self.gf_struct_solver[ish] = {} self.gf_struct_solver[ish] = {}
gf_struct_temp = [] gf_struct_temp = []
a_list = [a for a,al in self.gf_struct_corr[self.invshellmap[ish]] ] block_ind_list = [block for block,inner in self.gf_struct_corr[self.invshellmap[ish]] ]
for a in a_list: for block in block_ind_list:
dm = dens_mat[ish][a] 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
offdiag = [] offdiag = []
@ -448,16 +444,16 @@ class SumkLDA:
for i in range(NBlocs): for i in range(NBlocs):
blocs[i].sort() blocs[i].sort()
self.gf_struct_solver[ish].update( [('%s_%s'%(a,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'%(a,i),blocs[i]) ) gf_struct_temp.append( ('%s_%s'%(block,i),blocs[i]) )
# map is the mapping of the blocs from the SK blocs to the CTQMC blocs: # map is the mapping of the blocs from the SK blocs to the CTQMC blocs:
self.map[ish][a] = range(len(dmbool)) self.map[ish][block] = range(len(dmbool))
for ibl in range(NBlocs): for ibl in range(NBlocs):
for j in range(len(blocs[ibl])): for j in range(len(blocs[ibl])):
self.map[ish][a][blocs[ibl][j]] = '%s_%s'%(a,ibl) self.map[ish][block][blocs[ibl][j]] = '%s_%s'%(block,ibl)
self.map_inv[ish]['%s_%s'%(a,ibl)] = a self.map_inv[ish]['%s_%s'%(block,ibl)] = block
# now calculate degeneracies of orbitals: # now calculate degeneracies of orbitals:
@ -512,7 +508,7 @@ class SumkLDA:
# define matrices for inequivalent shells: # define matrices for inequivalent shells:
eff_atlevels = [ {} for ish in range(self.n_inequiv_corr_shells) ] eff_atlevels = [ {} for ish in range(self.n_inequiv_corr_shells) ]
for ish in range(self.n_inequiv_corr_shells): for ish in range(self.n_inequiv_corr_shells):
for bn in self.block_names[self.corr_shells[self.invshellmap[ish]][4]]: for bn in self.spin_block_names[self.corr_shells[self.invshellmap[ish]][4]]:
eff_atlevels[ish][bn] = numpy.identity(self.corr_shells[self.invshellmap[ish]][3], numpy.complex_) eff_atlevels[ish][bn] = numpy.identity(self.corr_shells[self.invshellmap[ish]][3], numpy.complex_)
# Chemical Potential: # Chemical Potential:
@ -530,20 +526,21 @@ 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 ish in range(self.n_corr_shells) ] self.Hsumk = [ {} for ish 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.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])
self.Hsumk[icrsh][bn] = numpy.zeros([dim,dim],numpy.complex_) self.Hsumk[icrsh][bn] = 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 ibn, bn in enumerate(self.block_names[self.corr_shells[icrsh][4]]): for ibl, bn in enumerate(self.spin_block_names[self.corr_shells[icrsh][4]]):
isp = self.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 xrange(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*ibn) * self.h_field * MMat MMat = self.hopping[ik,isp,0:n_orb,0:n_orb] - (1-2*ibl) * self.h_field * MMat
self.Hsumk[icrsh][bn] += self.bz_weights[ik] * numpy.dot( numpy.dot(self.proj_mat[ik,isp,icrsh,0:dim,0:n_orb],MMat), projmat = self.proj_mat[ik,isp,icrsh,0:dim,0:n_orb]
self.proj_mat[ik,isp,icrsh,0:dim,0:n_orb].conjugate().transpose() ) self.Hsumk[icrsh][bn] += self.bz_weights[ik] * numpy.dot( numpy.dot(projmat,MMat),
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)
@ -604,24 +601,24 @@ class SumkLDA:
blname = self.gf_struct_corr[icrsh][j][0] blname = self.gf_struct_corr[icrsh][j][0]
Ncr[blname] = 0.0 Ncr[blname] = 0.0
for a,al in self.gf_struct_solver[iorb].iteritems(): for block,inner in self.gf_struct_solver[iorb].iteritems():
bl = self.map_inv[iorb][a] bl = self.map_inv[iorb][block]
Ncr[bl] += dens_mat[a].real.trace() Ncr[bl] += dens_mat[block].real.trace()
M = self.corr_shells[icrsh][3] M = self.corr_shells[icrsh][3]
Ncrtot = 0.0 Ncrtot = 0.0
a_list = [a for a,al in self.gf_struct_corr[icrsh]] block_ind_list = [block for block,inner in self.gf_struct_corr[icrsh]]
for bl in a_list: for bl in block_ind_list:
Ncrtot += Ncr[bl] Ncrtot += Ncr[bl]
# 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 a_list: for bl in block_ind_list:
Ncr[bl] = Ncrtot / len(a_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 # 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 a_list: for bl in block_ind_list:
Ncr[bl] = Ncrtot / 2.0 Ncr[bl] = Ncrtot / 2.0
if (use_val is None): if (use_val is None):
@ -629,7 +626,7 @@ class SumkLDA:
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 a_list: for bl in block_ind_list:
Uav = U_interact*(Ncrtot-0.5) - J_hund*(Ncr[bl] - 0.5) Uav = U_interact*(Ncrtot-0.5) - J_hund*(Ncr[bl] - 0.5)
self.dc_imp[icrsh][bl] *= Uav self.dc_imp[icrsh][bl] *= Uav
self.dc_energ[icrsh] -= J_hund / 2.0 * (Ncr[bl]) * (Ncr[bl]-1.0) self.dc_energ[icrsh] -= J_hund / 2.0 * (Ncr[bl]) * (Ncr[bl]-1.0)
@ -638,7 +635,7 @@ class SumkLDA:
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 + (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) self.dc_energ[icrsh] = (U_interact + (M-1)*(U_interact-2.0*J_hund) + (M-1)*(U_interact-3.0*J_hund))/(2*M-1) / 2.0 * Ncrtot * (Ncrtot-1.0)
for bl in a_list: 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) 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 self.dc_imp[icrsh][bl] *= Uav
mpi.report("DC for shell %(icrsh)i and block %(bl)s = %(Uav)f"%locals()) mpi.report("DC for shell %(icrsh)i and block %(bl)s = %(Uav)f"%locals())
@ -646,7 +643,7 @@ class SumkLDA:
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 a_list: for bl in block_ind_list:
Uav = U_interact*(Ncrtot - Ncr[bl]/M) - J_hund * (Ncr[bl] - Ncr[bl]/M) Uav = U_interact*(Ncrtot - Ncr[bl]/M) - J_hund * (Ncr[bl] - Ncr[bl]/M)
self.dc_imp[icrsh][bl] *= Uav self.dc_imp[icrsh][bl] *= Uav
self.dc_energ[icrsh] -= (U_interact + (M-1)*J_hund)/M * 0.5 * Ncr[bl] * Ncr[bl] self.dc_energ[icrsh] -= (U_interact + (M-1)*J_hund)/M * 0.5 * Ncr[bl] * Ncr[bl]
@ -657,8 +654,8 @@ class SumkLDA:
else: else:
a_list = [a for a,al in self.gf_struct_corr[icrsh]] block_ind_list = [block for block,inner in self.gf_struct_corr[icrsh]]
for bl in a_list: for bl in block_ind_list:
self.dc_imp[icrsh][bl] *= use_val self.dc_imp[icrsh][bl] *= use_val
self.dc_energ[icrsh] = use_val * Ncrtot self.dc_energ[icrsh] = use_val * Ncrtot
@ -676,13 +673,13 @@ class SumkLDA:
assert len(Sigma_imp)==self.n_inequiv_corr_shells, "give exactly one Sigma for each inequivalent corr. shell!" assert len(Sigma_imp)==self.n_inequiv_corr_shells, "give exactly one Sigma for each inequivalent corr. shell!"
# init self.Sigma_imp: # init self.Sigma_imp:
if all(type(g) == GfImFreq for name,g 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 = [ (a,GfImFreq(indices = al, mesh = Sigma_imp[0].mesh)) for a,al in self.gf_struct_corr[i] ], self.Sigma_imp = [ BlockGf( name_block_generator = [ (block,GfImFreq(indices = inner, mesh = Sigma_imp[0].mesh)) for block,inner in self.gf_struct_corr[i] ],
make_copies = False) for i in xrange(self.n_corr_shells) ] make_copies = False) for i in xrange(self.n_corr_shells) ]
elif all(type(g) == GfReFreq for name,g 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 = [ (a,GfReFreq(indices = al, mesh = Sigma_imp[0].mesh)) for a,al in self.gf_struct_corr[i] ], self.Sigma_imp = [ BlockGf( name_block_generator = [ (block,GfReFreq(indices = inner, mesh = Sigma_imp[0].mesh)) for block,inner in self.gf_struct_corr[i] ],
make_copies = False) for i in xrange(self.n_corr_shells) ] make_copies = False) for i in xrange(self.n_corr_shells) ]
else: else:
raise ValueError, "This type of Sigma is not handled." raise ValueError, "This type of Sigma is not handled."
@ -697,26 +694,26 @@ class SumkLDA:
for blname in self.map[s]: for blname in self.map[s]:
cnt[blname] = 0 cnt[blname] = 0
for a,al in self.gf_struct_solver[s].iteritems(): for block,inner in self.gf_struct_solver[s].iteritems():
blname = self.map_inv[s][a] blname = self.map_inv[s][block]
map_ind[a] = range(len(al)) map_ind[block] = range(len(inner))
for i in al: for i in inner:
map_ind[a][i] = cnt[blname] map_ind[block][i] = cnt[blname]
cnt[blname]+=1 cnt[blname]+=1
for bl, orblist in self.gf_struct_solver[s].iteritems(): for block,inner in self.gf_struct_solver[s].iteritems():
for i in range(len(orblist)): for i in range(len(inner)):
for j in range(len(orblist)): for j in range(len(inner)):
ind1 = orblist[i] ind1 = inner[i]
ind2 = orblist[j] ind2 = inner[j]
ind1_imp = map_ind[bl][ind1] ind1_imp = map_ind[block][ind1]
ind2_imp = map_ind[bl][ind2] ind2_imp = map_ind[block][ind2]
self.Sigma_imp[icrsh][self.map_inv[s][bl]][ind1_imp,ind2_imp] << Sigma_imp[s][bl][ind1,ind2] self.Sigma_imp[icrsh][self.map_inv[s][block]][ind1_imp,ind2_imp] << Sigma_imp[s][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 xrange(self.n_corr_shells):
for sig,gf in self.Sigma_imp[icrsh]: self.Sigma_imp[icrsh][sig] << self.rotloc(icrsh,gf,direction='toGlobal') for bname,gf in self.Sigma_imp[icrsh]: self.Sigma_imp[icrsh][bname] << self.rotloc(icrsh,gf,direction='toGlobal')
@ -726,10 +723,10 @@ 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 xrange(self.n_corr_shells):
for bl,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][bl],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()))
sres[icrsh][bl] -= dccont sres[icrsh][bname] -= dccont
return sres # list of self energies corrected by DC return sres # list of self energies corrected by DC
@ -741,36 +738,6 @@ class SumkLDA:
self.chemical_potential = mu self.chemical_potential = mu
def inequiv_shells(self,lst):
"""
The number of inequivalent shells is calculated from lst, and a mapping is given as
map(i_corr_shells) = i_inequiv_corr_shells
invmap(i_inequiv_corr_shells) = i_corr_shells
in order to put the Self energies to all equivalent shells, and for extracting Gloc
"""
tmp = []
self.shellmap = [0 for i in range(len(lst))]
self.invshellmap = [0]
self.n_inequiv_corr_shells = 1
tmp.append( lst[0][1:3] )
if (len(lst)>1):
for i in range(len(lst)-1):
fnd = False
for j in range(self.n_inequiv_corr_shells):
if (tmp[j]==lst[i+1][1:3]):
fnd = True
self.shellmap[i+1] = j
if (fnd==False):
self.shellmap[i+1] = self.n_inequiv_corr_shells
self.n_inequiv_corr_shells += 1
tmp.append( lst[i+1][1:3] )
self.invshellmap.append(i+1)
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 mu. Since in general n_orbitals depends on k,
@ -839,7 +806,7 @@ class SumkLDA:
for icrsh in xrange(self.n_corr_shells): for icrsh in xrange(self.n_corr_shells):
tmp = Gloc[icrsh].copy() # init temporary storage tmp = Gloc[icrsh].copy() # init temporary storage
for sig,gf in tmp: tmp[sig] << self.downfold(ik,icrsh,sig,S[sig],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:
@ -855,10 +822,10 @@ class SumkLDA:
# 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 xrange(self.n_corr_shells):
for sig,gf in Gloc[icrsh]: Gloc[icrsh][sig] << 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 = [ (a,GfImFreq(indices = al, mesh = Gloc[0].mesh)) for a,al in self.gf_struct_solver[i].iteritems() ], Glocret = [ BlockGf( name_block_generator = [ (block,GfImFreq(indices = inner, mesh = Gloc[0].mesh)) for block,inner in self.gf_struct_solver[i].iteritems() ],
make_copies = False) for i in xrange(self.n_inequiv_corr_shells) ] make_copies = False) for i in xrange(self.n_inequiv_corr_shells) ]
for ish in xrange(self.n_inequiv_corr_shells): for ish in xrange(self.n_inequiv_corr_shells):
@ -868,21 +835,21 @@ class SumkLDA:
for blname in self.map[ish]: for blname in self.map[ish]:
cnt[blname] = 0 cnt[blname] = 0
for a,al in self.gf_struct_solver[ish].iteritems(): for block,inner in self.gf_struct_solver[ish].iteritems():
blname = self.map_inv[ish][a] blname = self.map_inv[ish][block]
map_ind[a] = range(len(al)) map_ind[block] = range(len(inner))
for i in al: for i in inner:
map_ind[a][i] = cnt[blname] map_ind[block][i] = cnt[blname]
cnt[blname]+=1 cnt[blname]+=1
for bl, orblist in self.gf_struct_solver[ish].iteritems(): for block,inner in self.gf_struct_solver[ish].iteritems():
for i in range(len(orblist)): for i in range(len(inner)):
for j in range(len(orblist)): for j in range(len(inner)):
ind1 = orblist[i] ind1 = inner[i]
ind2 = orblist[j] ind2 = inner[j]
ind1_imp = map_ind[bl][ind1] ind1_imp = map_ind[block][ind1]
ind2_imp = map_ind[bl][ind2] ind2_imp = map_ind[block][ind2]
Glocret[ish][bl][ind1,ind2] << Gloc[self.invshellmap[ish]][self.map_inv[ish][bl]][ind1_imp,ind2_imp] Glocret[ish][block][ind1,ind2] << Gloc[self.invshellmap[ish]][self.map_inv[ish][block]][ind1_imp,ind2_imp]
# return only the inequivalent shells: # return only the inequivalent shells:
@ -894,31 +861,31 @@ class SumkLDA:
assert (type(filename)==StringType), "filename has to be a string!" assert (type(filename)==StringType), "filename has to be a string!"
ntoi = self.names_to_ind[self.SO] ntoi = self.spin_names_to_ind[self.SO]
bln = self.block_names[self.SO] bln = self.spin_block_names[self.SO]
# Set up deltaN: # Set up deltaN:
deltaN = {} deltaN = {}
for ib in bln: for b in bln:
deltaN[ib] = [ numpy.zeros( [self.n_orbitals[ik,ntoi[ib]],self.n_orbitals[ik,ntoi[ib]]], numpy.complex_) for ik in range(self.n_k)] deltaN[b] = [ numpy.zeros( [self.n_orbitals[ik,ntoi[b]],self.n_orbitals[ik,ntoi[b]]], numpy.complex_) for ik in range(self.n_k)]
ikarray=numpy.array(range(self.n_k)) ikarray=numpy.array(range(self.n_k))
dens = {} dens = {}
for ib in bln: for b in bln:
dens[ib] = 0.0 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 sig,g in S: for bname,gf in S:
deltaN[sig][ik] = S[sig].density() deltaN[bname][ik] = S[bname].density()
dens[sig] += self.bz_weights[ik] * S[sig].total_density() dens[bname] += self.bz_weights[ik] * S[bname].total_density()
#put mpi Barrier: #put mpi Barrier:
for sig in deltaN: for bname in deltaN:
for ik in range(self.n_k): for ik in range(self.n_k):
deltaN[sig][ik] = mpi.all_reduce(mpi.world,deltaN[sig][ik],lambda x,y : x+y) deltaN[bname][ik] = mpi.all_reduce(mpi.world,deltaN[bname][ik],lambda x,y : x+y)
dens[sig] = mpi.all_reduce(mpi.world,dens[sig],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:
@ -967,7 +934,6 @@ class SumkLDA:
return deltaN, dens return deltaN, dens
################ ################
# FIXME LEAVE UNDOCUMENTED # FIXME LEAVE UNDOCUMENTED
################ ################
@ -1038,7 +1004,8 @@ class SumkLDA:
for ish in range(self.n_corr_shells): for ish in range(self.n_corr_shells):
dim = self.corr_shells[ish][3] dim = self.corr_shells[ish][3]
n_orb = self.n_orbitals[ik,0] n_orb = self.n_orbitals[ik,0]
dens_mat[ish][:,:] += numpy.dot(self.proj_mat[ik,0,ish,0:dim,0:n_orb],self.proj_mat[ik,0,ish,0:dim,0:n_orb].transpose().conjugate()) * self.bz_weights[ik] projmat = self.proj_mat[ik,0,ish,0:dim,0:n_orb]
dens_mat[ish][:,:] += numpy.dot(projmat, projmat.transpose().conjugate()) * self.bz_weights[ik]
if (self.symm_op!=0): dens_mat = self.symmcorr.symmetrize(dens_mat) if (self.symm_op!=0): dens_mat = self.symmcorr.symmetrize(dens_mat)
@ -1078,3 +1045,35 @@ class SumkLDA:
if atomlst[i+1]>atomlst[i]: atoms += 1 if atomlst[i+1]>atomlst[i]: atoms += 1
return atoms return atoms
##############################
# DUPLICATED, NEED TO REMOVE #
##############################
def inequiv_shells(self,lst):
"""
The number of inequivalent shells is calculated from lst, and a mapping is given as
map(i_corr_shells) = i_inequiv_corr_shells
invmap(i_inequiv_corr_shells) = i_corr_shells
in order to put the Self energies to all equivalent shells, and for extracting Gloc
"""
tmp = []
self.shellmap = [0 for i in range(len(lst))]
self.invshellmap = [0]
self.n_inequiv_corr_shells = 1
tmp.append( lst[0][1:3] )
if (len(lst)>1):
for i in range(len(lst)-1):
fnd = False
for j in range(self.n_inequiv_corr_shells):
if (tmp[j]==lst[i+1][1:3]):
fnd = True
self.shellmap[i+1] = j
if (fnd==False):
self.shellmap[i+1] = self.n_inequiv_corr_shells
self.n_inequiv_corr_shells += 1
tmp.append( lst[i+1][1:3] )
self.invshellmap.append(i+1)

151
python/sumk_lda_tools.py
View File

@ -41,11 +41,11 @@ class SumkLDATools(SumkLDA):
symmpar_data=symmpar_data, bands_data=bands_data) symmpar_data=symmpar_data, bands_data=bands_data)
def downfold_pc(self,ik,ir,ish,sig,gf_to_downfold,gf_inp): def downfold_pc(self,ik,ir,ish,bname,gf_to_downfold,gf_inp):
"""Downfolding a block of the Greens function""" """Downfolding a block of the Greens function"""
gf_downfolded = gf_inp.copy() gf_downfolded = gf_inp.copy()
isp = self.names_to_ind[self.SO][sig] # get spin index for proj. matrices isp = self.spin_names_to_ind[self.SO][bname] # get spin index for proj. matrices
dim = self.shells[ish][3] dim = self.shells[ish][3]
n_orb = self.n_orbitals[ik,isp] n_orb = self.n_orbitals[ik,isp]
L=self.proj_mat_pc[ik,isp,ish,ir,0:dim,0:n_orb] L=self.proj_mat_pc[ik,isp,ish,ir,0:dim,0:n_orb]
@ -85,48 +85,48 @@ class SumkLDATools(SumkLDA):
"""Calculates the lattice Green function on the real frequency axis. If self energy is """Calculates the lattice Green function on the real frequency axis. If self energy is
present and with_Sigma=True, the mesh is taken from Sigma. Otherwise, the mesh has to be given.""" present and with_Sigma=True, the mesh is taken from Sigma. Otherwise, the mesh has to be given."""
ntoi = self.names_to_ind[self.SO] ntoi = self.spin_names_to_ind[self.SO]
bln = self.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(g) == GfReFreq for name,g 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
BS = [ range(self.n_orbitals[ik,ntoi[ib]]) for ib in bln ] block_structure = [ range(self.n_orbitals[ik,ntoi[b]]) for b in bln ]
gf_struct = [ (bln[ib], BS[ib]) for ib in range(self.n_spin_blocks_gf[self.SO]) ] gf_struct = [ (bln[ibl], block_structure[ibl]) for ibl in range(self.n_spin_blocks[self.SO]) ]
a_list = [a for a,al in gf_struct] block_ind_list = [block for block,inner in gf_struct]
if (with_Sigma): if (with_Sigma):
glist = lambda : [ GfReFreq(indices = al, mesh =self.Sigma_imp[0].mesh) for a,al 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 = al, window=(mesh[0],mesh[1]),n_points=mesh[2]) for a,al 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 = a_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 sig,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[ib]]]==GFsize[ib] unchangedsize = all( [ self.n_orbitals[ik,ntoi[bln[ibl]]]==GFsize[ibl]
for ib in range(self.n_spin_blocks_gf[self.SO]) ] ) for ibl in range(self.n_spin_blocks[self.SO]) ] )
if (not unchangedsize): if (not unchangedsize):
BS = [ range(self.n_orbitals[ik,ntoi[ib]]) for ib in bln ] block_structure = [ range(self.n_orbitals[ik,ntoi[b]]) for b in bln ]
gf_struct = [ (bln[ib], BS[ib]) for ib in range(self.n_spin_blocks_gf[self.SO]) ] gf_struct = [ (bln[ibl], block_structure[ibl]) for ibl in range(self.n_spin_blocks[self.SO]) ]
a_list = [a for a,al in gf_struct] block_ind_list = [block for block,inner in gf_struct]
if (with_Sigma): if (with_Sigma):
glist = lambda : [ GfReFreq(indices = al, mesh =self.Sigma_imp[0].mesh) for a,al 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 = al, window=(mesh[0],mesh[1]),n_points=mesh[2]) for a,al 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 = a_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[bl]],numpy.complex_) for bl in bln] idmat = [numpy.identity(self.n_orbitals[ik,ntoi[b]],numpy.complex_) for b in bln]
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_gf[self.SO]): for ibl in range(self.n_spin_blocks[self.SO]):
ind = ntoi[bln[ibl]] ind = ntoi[bln[ibl]]
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[ibl] = self.hopping[ik,ind,0:n_orb,0:n_orb] - (idmat[ibl]*mu) - (idmat[ibl] * self.h_field * (1-2*ibl))
@ -135,7 +135,7 @@ class SumkLDATools(SumkLDA):
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 xrange(self.n_corr_shells):
for sig,gf in tmp: tmp[sig] << self.upfold(ik,icrsh,sig,stmp[icrsh][sig],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
self.G_upfold_refreq.invert() self.G_upfold_refreq.invert()
@ -152,13 +152,13 @@ class SumkLDATools(SumkLDA):
for i in range(n_om): om_mesh[i] = om_min + delta_om * i for i in range(n_om): om_mesh[i] = om_min + delta_om * i
DOS = {} DOS = {}
for bn in self.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_corr_shells) ] DOSproj = [ {} for icrsh in range(self.n_inequiv_corr_shells) ]
DOSproj_orb = [ {} for icrsh in range(self.n_inequiv_corr_shells) ] DOSproj_orb = [ {} for icrsh in range(self.n_inequiv_corr_shells) ]
for icrsh in range(self.n_inequiv_corr_shells): for icrsh in range(self.n_inequiv_corr_shells):
for bn in self.block_names[self.corr_shells[self.invshellmap[icrsh]][4]]: for bn in self.spin_block_names[self.corr_shells[self.invshellmap[icrsh]][4]]:
dl = self.corr_shells[self.invshellmap[icrsh]][3] dl = self.corr_shells[self.invshellmap[icrsh]][3]
DOSproj[icrsh][bn] = numpy.zeros([n_om],numpy.float_) DOSproj[icrsh][bn] = numpy.zeros([n_om],numpy.float_)
DOSproj_orb[icrsh][bn] = numpy.zeros([n_om,dl,dl],numpy.float_) DOSproj_orb[icrsh][bn] = numpy.zeros([n_om,dl,dl],numpy.float_)
@ -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 = [a for a,al in self.gf_struct_corr[icrsh]] b_list = [block for block,inner in self.gf_struct_corr[icrsh]]
#glist = lambda : [ GfReFreq(indices = al, beta = beta, mesh_array = mesh) for a,al in self.gf_struct_corr[icrsh]] #glist = lambda : [ GfReFreq(indices = inner, beta = beta, mesh_array = mesh) for block,inner in self.gf_struct_corr[icrsh]]
glist = lambda : [ GfReFreq(indices = al, window = (om_min,om_max), n_points = n_om) for a,al in self.gf_struct_corr[icrsh]] glist = lambda : [ GfReFreq(indices = inner, window = (om_min,om_max), n_points = n_om) for block,inner in self.gf_struct_corr[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 xrange(self.n_corr_shells): Gloc[icrsh].zero() # initialize to zero
@ -179,13 +179,13 @@ class SumkLDATools(SumkLDA):
# non-projected DOS # non-projected DOS
for iom in range(n_om): for iom in range(n_om):
for sig,gf in G_upfold: for bname,gf in G_upfold:
asd = gf.data[iom,:,:].imag.trace()/(-3.1415926535) asd = gf.data[iom,:,:].imag.trace()/(-3.1415926535)
DOS[sig][iom] += asd DOS[bname][iom] += asd
for icrsh in xrange(self.n_corr_shells): for icrsh in xrange(self.n_corr_shells):
tmp = Gloc[icrsh].copy() tmp = Gloc[icrsh].copy()
for sig,gf in tmp: tmp[sig] << self.downfold(ik,icrsh,sig,G_upfold[sig],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
@ -194,18 +194,18 @@ class SumkLDATools(SumkLDA):
if (self.use_rotations): if (self.use_rotations):
for icrsh in xrange(self.n_corr_shells): for icrsh in xrange(self.n_corr_shells):
for sig,gf in Gloc[icrsh]: Gloc[icrsh][sig] << 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
for ish in range(self.n_inequiv_corr_shells): for ish in range(self.n_inequiv_corr_shells):
for sig,gf in Gloc[self.invshellmap[ish]]: # loop over spins for bname,gf in Gloc[self.invshellmap[ish]]: # loop over spins
for iom in range(n_om): DOSproj[ish][sig][iom] += gf.data[iom,:,:].imag.trace()/(-3.1415926535) for iom in range(n_om): DOSproj[ish][bname][iom] += gf.data[iom,:,:].imag.trace()/(-3.1415926535)
DOSproj_orb[ish][sig][:,:,:] += 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.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]))
f.close() f.close()
@ -241,16 +241,14 @@ class SumkLDATools(SumkLDA):
assert hasattr(self,"Sigma_imp"), "Set Sigma First!!" assert hasattr(self,"Sigma_imp"), "Set Sigma First!!"
#things_to_read = ['Dens_Mat_below','N_parproj','Proj_Mat_pc','rotmat_all']
#value_read = self.read_input_from_HDF(SubGrp=self.parproj_data, things_to_read=things_to_read)
value_read = self.read_parproj_input_from_hdf() value_read = self.read_parproj_input_from_hdf()
if not value_read: return value_read if not value_read: return value_read
if self.symm_op: self.Symm_par = Symmetry(self.hdf_file,subgroup=self.symmpar_data) if self.symm_op: self.symmpar = Symmetry(self.hdf_file,subgroup=self.symmpar_data)
mu = self.chemical_potential mu = self.chemical_potential
gf_struct_proj = [ [ (al, range(self.shells[i][3])) for al in self.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 xrange(self.n_shells) ]
Gproj = [BlockGf(name_block_generator = [ (a,GfReFreq(indices = al, mesh = self.Sigma_imp[0].mesh)) for a,al 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 xrange(self.n_shells)]
for ish in range(self.n_shells): Gproj[ish].zero() for ish in range(self.n_shells): Gproj[ish].zero()
@ -258,13 +256,13 @@ class SumkLDATools(SumkLDA):
n_om = len(Msh) n_om = len(Msh)
DOS = {} DOS = {}
for bn in self.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 ish in range(self.n_shells) ] DOSproj = [ {} for ish in range(self.n_shells) ]
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.block_names[self.SO]: for bn in self.spin_block_names[self.SO]:
dl = self.shells[ish][3] dl = 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,dl,dl],numpy.float_)
@ -278,38 +276,38 @@ class SumkLDATools(SumkLDA):
# non-projected DOS # non-projected DOS
for iom in range(n_om): for iom in range(n_om):
for sig,gf in S: DOS[sig][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 xrange(self.n_shells):
tmp = Gproj[ish].copy() tmp = Gproj[ish].copy()
for ir in xrange(self.n_parproj[ish]): for ir in xrange(self.n_parproj[ish]):
for sig,gf in tmp: tmp[sig] << self.downfold_pc(ik,ir,ish,sig,S[sig],gf) for 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 sig in DOS: for bname in DOS:
DOS[sig] = mpi.all_reduce(mpi.world,DOS[sig],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 xrange(self.n_shells):
Gproj[ish] << mpi.all_reduce(mpi.world,Gproj[ish],lambda x,y : x+y) Gproj[ish] << mpi.all_reduce(mpi.world,Gproj[ish],lambda x,y : x+y)
mpi.barrier() mpi.barrier()
if (self.symm_op!=0): Gproj = self.Symm_par.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 xrange(self.n_shells):
for sig,gf in Gproj[ish]: Gproj[ish][sig] << 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):
for sig,gf in Gproj[ish]: for bname,gf in Gproj[ish]:
for iom in range(n_om): DOSproj[ish][sig][iom] += gf.data[iom,:,:].imag.trace()/(-3.1415926535) for iom in range(n_om): DOSproj[ish][bname][iom] += gf.data[iom,:,:].imag.trace()/(-3.1415926535)
DOSproj_orb[ish][sig][:,:,:] += 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.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')
for i in range(n_om): f.write("%s %s\n"%(Msh[i],DOS[bn][i])) for i in range(n_om): f.write("%s %s\n"%(Msh[i],DOS[bn][i]))
f.close() f.close()
@ -369,7 +367,7 @@ class SumkLDATools(SumkLDA):
# calculate A(k,w): # calculate A(k,w):
mu = self.chemical_potential mu = self.chemical_potential
bln = self.block_names[self.SO] bln = 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]
@ -384,20 +382,20 @@ class SumkLDATools(SumkLDA):
if (ishell is None): if (ishell is None):
Akw = {} Akw = {}
for ibn in bln: Akw[ibn] = 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:
Akw = {} Akw = {}
for ibn in bln: Akw[ibn] = numpy.zeros([self.shells[ishell][3],self.n_k, n_om ],numpy.float_) for b in bln: Akw[b] = 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 ibn in bln: Akw[ibn] = numpy.zeros([self.n_k,1],numpy.float_) for b in bln: Akw[b] = numpy.zeros([self.n_k,1],numpy.float_)
if not (ishell is None): if not (ishell is None):
GFStruct_proj = [ (al, range(self.shells[ishell][3])) for al in bln ] GFStruct_proj = [ (b, range(self.shells[ishell][3])) for b in bln ]
Gproj = BlockGf(name_block_generator = [ (a,GfReFreq(indices = al, mesh = self.Sigma_imp[0].mesh)) for a,al 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 xrange(self.n_k):
@ -408,10 +406,10 @@ 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):
if fermi_surface: if fermi_surface:
for sig,gf in S: Akw[sig][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:
for sig,gf in S: Akw[sig][ik,iom] += gf.data[iom,:,:].imag.trace()/(-3.1415926535) for bname,gf in S: Akw[bname][ik,iom] += gf.data[iom,:,:].imag.trace()/(-3.1415926535)
Akw[sig][ik,iom] += ik*shift # shift Akw for plotting in xmgrace -- REMOVE Akw[bname][ik,iom] += ik*shift # shift Akw for plotting in xmgrace -- REMOVE
else: else:
@ -419,14 +417,14 @@ class SumkLDATools(SumkLDA):
Gproj.zero() Gproj.zero()
tmp = Gproj.copy() tmp = Gproj.copy()
for ir in xrange(self.n_parproj[ishell]): for ir in xrange(self.n_parproj[ishell]):
for sig,gf in tmp: tmp[sig] << self.downfold_pc(ik,ir,ishell,sig,S[sig],gf) for bname,gf in tmp: tmp[bname] << self.downfold_pc(ik,ir,ishell,bname,S[bname],gf)
Gproj += tmp Gproj += tmp
# FIXME NEED TO READ IN ROTMAT_ALL FROM PARPROJ SUBGROUP, REPLACE ROTLOC WITH ROTLOC_ALL # FIXME NEED TO READ IN ROTMAT_ALL FROM PARPROJ SUBGROUP, REPLACE ROTLOC WITH ROTLOC_ALL
# TO BE FIXED: # TO BE FIXED:
# rotate to local frame # rotate to local frame
#if (self.use_rotations): #if (self.use_rotations):
# for sig,gf in Gproj: Gproj[sig] << 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):
@ -501,27 +499,24 @@ class SumkLDATools(SumkLDA):
"""Calculates the orbitally-resolved density matrix for all the orbitals considered in the input. """Calculates the orbitally-resolved density matrix for all the orbitals considered in the input.
The theta-projectors are used, hence case.parproj data is necessary""" The theta-projectors are used, hence case.parproj data is necessary"""
#things_to_read = ['Dens_Mat_below','N_parproj','Proj_Mat_pc','rotmat_all']
#value_read = self.read_input_from_HDF(SubGrp=self.parproj_data,things_to_read=things_to_read)
value_read = self.read_parproj_input_from_hdf() value_read = self.read_parproj_input_from_hdf()
if not value_read: return value_read if not value_read: return value_read
if self.symm_op: self.Symm_par = 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.block_names[self.SO] bln = self.spin_block_names[self.SO]
ntoi = self.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(bln)) ] # init the density matrix
mu = self.chemical_potential mu = self.chemical_potential
GFStruct_proj = [ [ (al, range(self.shells[i][3])) for al in bln ] for i in xrange(self.n_shells) ] GFStruct_proj = [ [ (b, range(self.shells[i][3])) for b in bln ] for i in xrange(self.n_shells) ]
if hasattr(self,"Sigma_imp"): if hasattr(self,"Sigma_imp"):
Gproj = [BlockGf(name_block_generator = [ (a,GfImFreq(indices = al, mesh = self.Sigma_imp[0].mesh)) for a,al 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 xrange(self.n_shells)]
beta = self.Sigma_imp[0].mesh.beta beta = self.Sigma_imp[0].mesh.beta
else: else:
Gproj = [BlockGf(name_block_generator = [ (a,GfImFreq(indices = al, beta = beta)) for a,al 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 xrange(self.n_shells)]
for ish in xrange(self.n_shells): Gproj[ish].zero() for ish in xrange(self.n_shells): Gproj[ish].zero()
@ -535,7 +530,7 @@ class SumkLDATools(SumkLDA):
for ish in xrange(self.n_shells): for ish in xrange(self.n_shells):
tmp = Gproj[ish].copy() tmp = Gproj[ish].copy()
for ir in xrange(self.n_parproj[ish]): for ir in xrange(self.n_parproj[ish]):
for sig,gf in tmp: tmp[sig] << self.downfold_pc(ik,ir,ish,sig,S[sig],gf) for 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:
@ -545,17 +540,17 @@ class SumkLDATools(SumkLDA):
# Symmetrisation: # Symmetrisation:
if (self.symm_op!=0): Gproj = self.Symm_par.symmetrize(Gproj) if (self.symm_op!=0): Gproj = self.symmpar.symmetrize(Gproj)
for ish in xrange(self.n_shells): for ish in xrange(self.n_shells):
# Rotation to local: # Rotation to local:
if (self.use_rotations): if (self.use_rotations):
for sig,gf in Gproj[ish]: Gproj[ish][sig] << 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
for sig,gf in Gproj[ish]: #dmg.append(Gproj[ish].density()[sig]) for bname,gf in Gproj[ish]: #dmg.append(Gproj[ish].density()[bname])
self.dens_mat_window[isp][ish] = Gproj[ish].density()[sig] self.dens_mat_window[isp][ish] = Gproj[ish].density()[bname]
isp+=1 isp+=1
# add Density matrices to get the total: # add Density matrices to get the total:

4
python/symmetry.py
View File

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

31
python/trans_basis.py
View File

@ -73,30 +73,29 @@ class TransBasis:
'''Rotates a given GF into the new basis.''' '''Rotates a given GF into the new basis.'''
# build a full GF # build a full GF
gfrotated = BlockGf( name_block_generator = [ (a,GfImFreq(indices = al, mesh = gf_to_rot.mesh)) for a,al in self.SK.gf_struct_corr[0] ], make_copies = False) gfrotated = BlockGf( name_block_generator = [ (block,GfImFreq(indices = inner, mesh = gf_to_rot.mesh)) for block,inner in self.SK.gf_struct_corr[0] ], make_copies = False)
# transform the CTQMC blocks to the full matrix: # transform the CTQMC blocks to the full matrix:
s = self.SK.shellmap[0] # s is the index of the inequivalent shell corresponding to icrsh s = self.SK.shellmap[0] # s is the index of the inequivalent shell corresponding to icrsh
for bl, orblist in self.gf_struct_solver[s].iteritems(): for block, inner in self.gf_struct_solver[s].iteritems():
for i in range(len(orblist)): for i in range(len(inner)):
for j in range(len(orblist)): for j in range(len(inner)):
ind1 = orblist[i] ind1 = inner[i]
ind2 = orblist[j] ind2 = inner[j]
gfrotated[self.SK.map_inv[s][bl]][ind1,ind2] << gf_to_rot[bl][ind1,ind2] gfrotated[self.SK.map_inv[s][block]][ind1,ind2] << gf_to_rot[block][ind1,ind2]
# Rotate using the matrix w # Rotate using the matrix w
for sig,bn in gfrotated: for bname,gf in gfrotated:
gfrotated[sig].from_L_G_R(self.w.transpose().conjugate(),gfrotated[sig],self.w) gfrotated[bname].from_L_G_R(self.w.transpose().conjugate(),gfrotated[bname],self.w)
gfreturn = gf_to_rot.copy() gfreturn = gf_to_rot.copy()
# Put back into CTQMC basis: # Put back into CTQMC basis:
for bl, orblist in self.gf_struct_solver[s].iteritems(): for block, inner in self.gf_struct_solver[s].iteritems():
for i in range(len(orblist)): for i in range(len(inner)):
for j in range(len(orblist)): for j in range(len(inner)):
ind1 = orblist[i] ind1 = inner[i]
ind2 = orblist[j] ind2 = inner[j]
gfreturn[bl][ind1,ind2] << gfrotated[self.SK.map_inv[0][bl]][ind1,ind2] gfreturn[block][ind1,ind2] << gfrotated[self.SK.map_inv[0][block]][ind1,ind2]
return gfreturn return gfreturn