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

Merged lattice_gf_matsubara and lattice_gf_realfreq into single function

Browse Source
This commit is contained in:
Priyanka Seth 2014-12-07 00:29:39 +01:00
parent cc3a9deaa8
commit 259fd64824
4 changed files with 214 additions and 347 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

67
python/lattice_gf.py
View File

@ -1,67 +0,0 @@
def lattice_gf(self, ik, mu, iw_or_w="iw", beta=40, broadening, mesh=None, with_Sigma=True, with_dc=True):
"""Calculates the lattice Green function from the DFT hopping and the self energy at k-point number ik
and chemical potential mu."""
ntoi = self.spin_names_to_ind[self.SO]
spn = self.spin_block_names[self.SO]
if (iw_or_w != "iw") and (iw_or_w != "w"):
raise ValueError, "lattice_gf: implemented only for Re/Im frequency functions."
# Are we including Sigma?
if with_Sigma:
if with_dc: sigma_inc_dc = self.add_dc()
Sigma_imp = getattr(self,"Sigma_imp_"+iw_or_w)
if iw_or_w == "iw": beta = Sigma_imp[0].mesh.beta # override beta if Sigma_iw is present
else:
if (iw_or_w == "w") and (mesh is None):
raise ValueError, "Give the mesh=(om_min,om_max,n_points) for the lattice GfReFreq."
# Check if G_upfold is present
set_up_G_upfold = False # Assume not
if not hasattr(self,"G_upfold_"+iw_or_w)):
set_up_G_upfold = True # Need to create G_upfold_iw
setattr(self,"G_upfold_"+iw_or_w,0) # Set temporarily to zero
else: # Check that existing GF is consistent
GFsize = [ gf.N1 for bname,gf in getattr(self,"G_upfold_"+iw_or_w)]
unchangedsize = all( [ self.n_orbitals[ik,ntoi[spn[isp]]]==GFsize[isp]
for isp in range(self.n_spin_blocks[self.SO]) ] )
if not unchangedsize: set_up_G_upfold = True
if (iw_or_w != "iw") and (self.G_upfold_iw.mesh.beta != beta): # additional check for ImFreq
set_up_G_upfold = True
G_upfold = getattr(self,"G_upfold_"+iw_or_w)
# Set up G_upfold
if set_up_G_upfold:
block_structure = [ range(self.n_orbitals[ik,ntoi[sp]]) for sp in spn ]
gf_struct = [ (spn[isp], block_structure[isp]) for isp in range(self.n_spin_blocks[self.SO]) ]
block_ind_list = [block for block,inner in gf_struct]
if with_Sigma:
glist = lambda : [ GfImFreq(indices=inner,mesh=Sigma_imp[0].mesh) for block,inner in gf_struct]
else:
if iw_or_w == "iw":
glist = lambda : [ GfImFreq(indices=inner,beta=beta) for block,inner in gf_struct]
elif iw_or_w == "w":
glist = lambda : [ GfImFreq(indices=inner,window=(mesh[0],mesh[1]),n_points=mesh[2]) for block,inner in gf_struct]
G_upfold = BlockGf(name_list = block_ind_list, block_list = glist(), make_copies = False)
G_upfold.zero()
if iw_or_w == "iw":
G_upfold << iOmega_n
elif iw_or_w == "w":
G_upfold << Omega + 1j*broadening
idmat = [numpy.identity(self.n_orbitals[ik,ntoi[sp]],numpy.complex_) for sp in spn]
M = copy.deepcopy(idmat)
for ibl in range(self.n_spin_blocks[self.SO]):
ind = ntoi[spn[ibl]]
n_orb = self.n_orbitals[ik,ind]
M[ibl] = self.hopping[ik,ind,0:n_orb,0:n_orb] - (idmat[ibl]*mu) - (idmat[ibl] * self.h_field * (1-2*ibl))
G_upfold -= M
if with_Sigma:
for icrsh in range(self.n_corr_shells):
for bname,gf in G_upfold: gf -= self.upfold(ik,icrsh,bname,sigma_inc_dc[icrsh][bname],gf)
G_upfold.invert()
return G_upfold

270
python/sumk_dft.py
View File

@ -49,7 +49,7 @@ class SumkDFT:
self.symmpar_data = symmpar_data self.symmpar_data = symmpar_data
self.bands_data = bands_data self.bands_data = bands_data
self.transp_data = transp_data self.transp_data = transp_data
self.G_upfold = None #self.G_latt_iw = None # DEBUG -- remove later
self.h_field = h_field self.h_field = h_field
# Read input from HDF: # Read input from HDF:
@ -128,7 +128,6 @@ class SumkDFT:
if (len(things_to_read) != 0): mpi.report("Loading failed: No %s subgroup in HDF5!"%subgrp) if (len(things_to_read) != 0): mpi.report("Loading failed: No %s subgroup in HDF5!"%subgrp)
subgroup_present = False subgroup_present = False
value_read = False value_read = False
del ar del ar
# now do the broadcasting: # now do the broadcasting:
for it in things_to_read: setattr(self,it,mpi.bcast(getattr(self,it))) for it in things_to_read: setattr(self,it,mpi.bcast(getattr(self,it)))
@ -163,7 +162,7 @@ class SumkDFT:
try: try:
list_to_return.append(ar[subgrp][it]) list_to_return.append(ar[subgrp][it])
except: except:
raise ValueError, "%s not found, and so not loaded."%it raise ValueError, "load: %s not found, and so not loaded."%it
del ar del ar
return list_to_return return list_to_return
@ -203,7 +202,7 @@ class SumkDFT:
"""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')),"rotloc: Give direction 'toLocal' or 'toGlobal'."
gf_rotated = gf_to_rotate.copy() gf_rotated = gf_to_rotate.copy()
@ -226,76 +225,99 @@ class SumkDFT:
return gf_rotated return gf_rotated
def lattice_gf_matsubara(self,ik,mu,beta=40,with_Sigma=True): def lattice_gf(self, ik, mu, iw_or_w="iw", beta=40, broadening=None, mesh=None, with_Sigma=True, with_dc=True):
"""Calculates the lattice Green function from the DFT hopping and the self energy at k-point number ik """Calculates the lattice Green function from the DFT hopping and the self energy at k-point number ik
and chemical potential mu.""" and chemical potential mu."""
ntoi = self.spin_names_to_ind[self.SO] ntoi = self.spin_names_to_ind[self.SO]
spn = self.spin_block_names[self.SO] spn = self.spin_block_names[self.SO]
if (iw_or_w != "iw") and (iw_or_w != "w"): raise ValueError, "lattice_gf: Implemented only for Re/Im frequency functions."
if not hasattr(self,"Sigma_imp_"+iw_or_w): with_Sigma = False
if broadening is None:
if mesh is None:
broadening = 0.01
else: # broadening = 2 * \Delta omega, where \Delta omega is the spacing of omega points
broadening = 2.0 * ( (mesh[1]-mesh[0])/(mesh[2]-1) )
if not hasattr(self,"Sigma_imp"): with_Sigma = False # Are we including Sigma?
if with_Sigma: if with_Sigma:
stmp = self.add_dc() if with_dc: sigma_minus_dc = self.add_dc(iw_or_w)
beta = self.Sigma_imp[0].mesh.beta # override beta if Sigma is present Sigma_imp = getattr(self,"Sigma_imp_"+iw_or_w)
if iw_or_w == "iw": beta = Sigma_imp[0].mesh.beta # override beta if Sigma_iw is present
else:
if (iw_or_w == "w") and (mesh is None):
raise ValueError, "lattice_gf: Give the mesh=(om_min,om_max,n_points) for the lattice GfReFreq."
# Do we need to set up G_upfold? # Check if G_latt is present
set_up_G_upfold = False # assume not set_up_G_latt = False # Assume not
if self.G_upfold is None: # yes if not G_upfold provided if not hasattr(self,"G_latt_"+iw_or_w):
set_up_G_upfold = True set_up_G_latt = True # Need to create G_latt_(i)w
else: # yes if inconsistencies present in existing G_upfold else: # Check that existing GF is consistent
GFsize = [ gf.N1 for bname,gf in self.G_upfold] G_latt = getattr(self,"G_latt_"+iw_or_w)
unchangedsize = all( [ self.n_orbitals[ik,ntoi[spn[isp]]] == GFsize[isp] GFsize = [ gf.N1 for bname,gf in G_latt]
for isp in range(self.n_spin_blocks[self.SO]) ] ) unchangedsize = all( [ self.n_orbitals[ik,ntoi[spn[isp]]]==GFsize[isp] for isp in range(self.n_spin_blocks[self.SO]) ] )
if (not unchangedsize) or (self.G_upfold.mesh.beta != beta): set_up_G_upfold = True if not unchangedsize: set_up_G_latt = True
if (iw_or_w == "iw") and (self.G_latt_iw.mesh.beta != beta): set_up_G_latt = True # additional check for ImFreq
# Set up G_upfold # Set up G_latt
if set_up_G_upfold: if set_up_G_latt:
block_structure = [ range(self.n_orbitals[ik,ntoi[sp]]) for sp in spn ] block_structure = [ range(self.n_orbitals[ik,ntoi[sp]]) for sp in spn ]
gf_struct = [ (spn[isp], block_structure[isp]) for isp in range(self.n_spin_blocks[self.SO]) ] gf_struct = [ (spn[isp], block_structure[isp]) for isp in range(self.n_spin_blocks[self.SO]) ]
block_ind_list = [block for block,inner in gf_struct] block_ind_list = [block for block,inner in gf_struct]
if with_Sigma: if iw_or_w == "iw":
glist = lambda : [ GfImFreq(indices = inner, mesh = self.Sigma_imp[0].mesh) for block,inner in gf_struct]
else:
glist = lambda : [ GfImFreq(indices=inner,beta=beta) for block,inner in gf_struct] 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) elif iw_or_w == "w":
self.G_upfold.zero() if with_Sigma:
glist = lambda : [ GfReFreq(indices=inner,mesh=Sigma_imp[0].mesh) for block,inner in gf_struct]
else:
glist = lambda : [ GfReFreq(indices=inner,window=(mesh[0],mesh[1]),n_points=mesh[2]) for block,inner in gf_struct]
G_latt = BlockGf(name_list = block_ind_list, block_list = glist(), make_copies = False)
G_latt.zero()
if iw_or_w == "iw":
G_latt << iOmega_n
elif iw_or_w == "w":
G_latt << Omega + 1j*broadening
self.G_upfold << iOmega_n
idmat = [numpy.identity(self.n_orbitals[ik,ntoi[sp]],numpy.complex_) for sp in spn] idmat = [numpy.identity(self.n_orbitals[ik,ntoi[sp]],numpy.complex_) for sp in spn]
M = copy.deepcopy(idmat) M = copy.deepcopy(idmat)
for isp in range(self.n_spin_blocks[self.SO]): for ibl in range(self.n_spin_blocks[self.SO]):
ind = ntoi[spn[isp]] ind = ntoi[spn[ibl]]
n_orb = self.n_orbitals[ik,ind] n_orb = self.n_orbitals[ik,ind]
M[isp] = self.hopping[ik,ind,0:n_orb,0:n_orb] - (idmat[isp]*mu) - (idmat[isp] * self.h_field * (1-2*isp)) 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 G_latt -= M
if with_Sigma: if with_Sigma:
for icrsh in range(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 G_latt: gf -= self.upfold(ik,icrsh,bname,sigma_minus_dc[icrsh][bname],gf)
self.G_upfold.invert() G_latt.invert()
setattr(self,"G_latt_"+iw_or_w,G_latt)
return self.G_upfold return G_latt
def put_Sigma(self, Sigma_imp): def put_Sigma(self, Sigma_imp):
"""Puts the impurity self energies for inequivalent atoms into the class, respects the multiplicity of the atoms.""" """Sets the impurity self energies for inequivalent atoms into the class, respects the multiplicity of the atoms."""
assert isinstance(Sigma_imp,list), "Sigma_imp has to be a list of Sigmas for the correlated shells, even if it is of length 1!" assert isinstance(Sigma_imp,list), "put_Sigma: Sigma_imp has to be a list of Sigmas for the correlated shells, even if it is of length 1!"
assert len(Sigma_imp) == self.n_inequiv_shells, "give exactly one Sigma for each inequivalent corr. shell!" assert len(Sigma_imp) == self.n_inequiv_shells, "put_Sigma: give exactly one Sigma for each inequivalent corr. shell!"
# init self.Sigma_imp: # init self.Sigma_imp_(i)w:
if all(type(gf) == GfImFreq for bname,gf in Sigma_imp[0]): 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[icrsh] ], self.Sigma_imp_iw = [ BlockGf( name_block_generator = [ (block,GfImFreq(indices = inner, mesh = Sigma_imp[0].mesh))
make_copies = False) for icrsh in range(self.n_corr_shells) ] for block,inner in self.gf_struct_sumk[icrsh] ], make_copies = False)
for icrsh in range(self.n_corr_shells) ]
SK_Sigma_imp = self.Sigma_imp_iw
elif all(type(gf) == GfReFreq for bname,gf in Sigma_imp[0]): 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[icrsh] ], self.Sigma_imp_w = [ BlockGf( name_block_generator = [ (block,GfReFreq(indices = inner, mesh = Sigma_imp[0].mesh))
make_copies = False) for icrsh in range(self.n_corr_shells) ] for block,inner in self.gf_struct_sumk[icrsh] ], make_copies = False)
for icrsh in range(self.n_corr_shells) ]
SK_Sigma_imp = self.Sigma_imp_w
else: else:
raise ValueError, "This type of Sigma is not handled." raise ValueError, "put_Sigma: This type of Sigma is not handled."
# transform the CTQMC blocks to the full matrix: # transform the CTQMC blocks to the full matrix:
for icrsh in range(self.n_corr_shells): for icrsh in range(self.n_corr_shells):
@ -305,12 +327,12 @@ class SumkDFT:
for ind2 in inner: for ind2 in inner:
block_sumk,ind1_sumk = self.solver_to_sumk[ish][(block,ind1)] block_sumk,ind1_sumk = self.solver_to_sumk[ish][(block,ind1)]
block_sumk,ind2_sumk = self.solver_to_sumk[ish][(block,ind2)] block_sumk,ind2_sumk = self.solver_to_sumk[ish][(block,ind2)]
self.Sigma_imp[icrsh][block_sumk][ind1_sumk,ind2_sumk] << Sigma_imp[ish][block][ind1,ind2] SK_Sigma_imp[icrsh][block_sumk][ind1_sumk,ind2_sumk] << Sigma_imp[ish][block][ind1,ind2]
# rotation from local to global coordinate system: # rotation from local to global coordinate system:
if self.use_rotations: if self.use_rotations:
for icrsh in range(self.n_corr_shells): for icrsh in range(self.n_corr_shells):
for bname,gf in self.Sigma_imp[icrsh]: self.Sigma_imp[icrsh][bname] << self.rotloc(icrsh, gf, direction = 'toGlobal') for bname,gf in SK_Sigma_imp[icrsh]: gf << self.rotloc(icrsh, gf, direction = 'toGlobal')
def extract_G_loc(self, mu = None, with_Sigma = True): def extract_G_loc(self, mu = None, with_Sigma = True):
@ -321,19 +343,19 @@ class SumkDFT:
""" """
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 range(self.n_corr_shells) ] # this list will be returned Gloc = [ self.Sigma_imp_iw[icrsh].copy() for icrsh in range(self.n_corr_shells) ] # this list will be returned
for icrsh in range(self.n_corr_shells): Gloc[icrsh].zero() # initialize to zero 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))
for ik in mpi.slice_array(ikarray): for ik in mpi.slice_array(ikarray):
S = self.lattice_gf_matsubara(ik = ik, mu = mu, with_Sigma = with_Sigma, beta = beta) G_latt_iw = self.lattice_gf(ik = ik, mu = mu, iw_or_w = "iw", with_Sigma = with_Sigma, beta = beta)
S *= self.bz_weights[ik] G_latt_iw *= self.bz_weights[ik]
for icrsh in range(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,G_latt_iw[bname],gf)
Gloc[icrsh] += tmp Gloc[icrsh] += tmp
# collect data from mpi # collect data from mpi
@ -350,19 +372,18 @@ class SumkDFT:
for bname,gf in Gloc[icrsh]: Gloc[icrsh][bname] << self.rotloc(icrsh,gf,direction = 'toLocal') 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() ], Gloc_inequiv = [ BlockGf( name_block_generator = [ (block,GfImFreq(indices = inner, mesh = Gloc[0].mesh)) for block,inner in self.gf_struct_solver[ish].iteritems() ],
make_copies = False) for ish in range(self.n_inequiv_shells) ] make_copies = False) for ish in range(self.n_inequiv_shells) ]
for ish in range(self.n_inequiv_shells): for ish in range(self.n_inequiv_shells):
for 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_sumk,ind1_sumk = self.solver_to_sumk[ish][(block,ind1)] block_sumk,ind1_sumk = self.solver_to_sumk[ish][(block,ind1)]
block_sumk,ind2_sumk = 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_sumk][ind1_sumk,ind2_sumk] Gloc_inequiv[ish][block][ind1,ind2] << Gloc[self.inequiv_to_corr[ish]][block_sumk][ind1_sumk,ind2_sumk]
# return only the inequivalent shells: # return only the inequivalent shells:
return Glocret return Gloc_inequiv
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):
@ -451,7 +472,7 @@ class SumkDFT:
def density_matrix(self, method = 'using_gf', beta = 40.0): def density_matrix(self, method = 'using_gf', beta = 40.0):
"""Calculate density matrices in one of two ways: """Calculate density matrices in one of two ways:
if 'using_gf': First get upfolded gf (g_loc is not set up), then density matrix. if 'using_gf': First get lattice gf (g_loc is not set up), then density matrix.
It is useful for Hubbard I, and very quick. It is useful for Hubbard I, and very quick.
No assumption on the hopping structure is made (ie diagonal or not). No assumption on the hopping structure is made (ie diagonal or not).
if 'using_point_integration': Only works for diagonal hopping matrix (true in wien2k). if 'using_point_integration': Only works for diagonal hopping matrix (true in wien2k).
@ -466,9 +487,9 @@ class SumkDFT:
if method == "using_gf": if method == "using_gf":
G_upfold = self.lattice_gf_matsubara(ik = ik, beta = beta, mu = self.chemical_potential) G_latt_iw = self.lattice_gf(ik = ik, mu = self.chemical_potential, iw_or_w = "iw", beta = beta)
G_upfold *= self.bz_weights[ik] G_latt_iw *= self.bz_weights[ik]
dm = G_upfold.density() dm = G_latt_iw.density()
MMat = [dm[sp] for sp in self.spin_block_names[self.SO]] MMat = [dm[sp] for sp in self.spin_block_names[self.SO]]
elif method == "using_point_integration": elif method == "using_point_integration":
@ -494,10 +515,10 @@ class SumkDFT:
for icrsh in range(self.n_corr_shells): for icrsh in range(self.n_corr_shells):
for isp, sp in enumerate(self.spin_block_names[self.corr_shells[icrsh]['SO']]): for isp, sp in enumerate(self.spin_block_names[self.corr_shells[icrsh]['SO']]):
isp = self.spin_names_to_ind[self.corr_shells[icrsh]['SO']][sp] ind = self.spin_names_to_ind[self.corr_shells[icrsh]['SO']][sp]
dim = self.corr_shells[icrsh]['dim'] dim = self.corr_shells[icrsh]['dim']
n_orb = self.n_orbitals[ik,isp] n_orb = self.n_orbitals[ik,ind]
projmat = self.proj_mat[ik,isp,icrsh,0:dim,0:n_orb] projmat = self.proj_mat[ik,ind,icrsh,0:dim,0:n_orb]
if method == "using_gf": if method == "using_gf":
dens_mat[icrsh][sp] += numpy.dot( numpy.dot(projmat,MMat[isp]), dens_mat[icrsh][sp] += numpy.dot( numpy.dot(projmat,MMat[isp]),
projmat.transpose().conjugate() ) projmat.transpose().conjugate() )
@ -507,8 +528,8 @@ class SumkDFT:
# 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 sp 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][sp] = mpi.all_reduce(mpi.world, dens_mat[icrsh][sp], 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)
@ -516,9 +537,9 @@ class SumkDFT:
# Rotate to local coordinate system: # Rotate to local coordinate system:
if self.use_rotations: if self.use_rotations:
for icrsh in range(self.n_corr_shells): for icrsh in range(self.n_corr_shells):
for bl in dens_mat[icrsh]: for sp in dens_mat[icrsh]:
if self.rot_mat_time_inv[icrsh] == 1: dens_mat[icrsh][bl] = dens_mat[icrsh][bl].conjugate() if self.rot_mat_time_inv[icrsh] == 1: dens_mat[icrsh][sp] = dens_mat[icrsh][sp].conjugate()
dens_mat[icrsh][bl] = numpy.dot( numpy.dot(self.rot_mat[icrsh].conjugate().transpose(),dens_mat[icrsh][bl]), dens_mat[icrsh][sp] = numpy.dot( numpy.dot(self.rot_mat[icrsh].conjugate().transpose(),dens_mat[icrsh][sp]),
self.rot_mat[icrsh] ) self.rot_mat[icrsh] )
return dens_mat return dens_mat
@ -533,54 +554,41 @@ class SumkDFT:
for ish in range(self.n_inequiv_shells): for ish in range(self.n_inequiv_shells):
for sp in self.spin_block_names[self.corr_shells[self.inequiv_to_corr[ish]]['SO']]: for sp in self.spin_block_names[self.corr_shells[self.inequiv_to_corr[ish]]['SO']]:
eff_atlevels[ish][sp] = numpy.identity(self.corr_shells[self.inequiv_to_corr[ish]]['dim'], numpy.complex_) eff_atlevels[ish][sp] = numpy.identity(self.corr_shells[self.inequiv_to_corr[ish]]['dim'], numpy.complex_)
eff_atlevels[ish][sp] *= -self.chemical_potential
# Chemical Potential: eff_atlevels[ish][sp] -= self.dc_imp[self.inequiv_to_corr[ish]][sp]
for ish in range(self.n_inequiv_shells):
for ii in eff_atlevels[ish]:
eff_atlevels[ish][ii] *= -self.chemical_potential
# double counting term:
for ish in range(self.n_inequiv_shells):
for ii in eff_atlevels[ish]:
eff_atlevels[ish][ii] -= self.dc_imp[self.inequiv_to_corr[ish]][ii]
# 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 icrsh in range(self.n_corr_shells) ] self.Hsumk = [ {} for icrsh in range(self.n_corr_shells) ]
for icrsh in range(self.n_corr_shells):
for sp in self.spin_block_names[self.corr_shells[icrsh]['SO']]:
dim = self.corr_shells[icrsh]['dim'] #*(1+self.corr_shells[icrsh]['SO'])
self.Hsumk[icrsh][sp] = numpy.zeros([dim,dim],numpy.complex_)
for icrsh in range(self.n_corr_shells): for icrsh in range(self.n_corr_shells):
dim = self.corr_shells[icrsh]['dim'] dim = self.corr_shells[icrsh]['dim']
for sp in self.spin_block_names[self.corr_shells[icrsh]['SO']]:
self.Hsumk[icrsh][sp] = numpy.zeros([dim,dim],numpy.complex_)
for isp, sp in enumerate(self.spin_block_names[self.corr_shells[icrsh]['SO']]): for isp, sp in enumerate(self.spin_block_names[self.corr_shells[icrsh]['SO']]):
isp = self.spin_names_to_ind[self.corr_shells[icrsh]['SO']][sp] ind = self.spin_names_to_ind[self.corr_shells[icrsh]['SO']][sp]
for ik in range(self.n_k): for ik in range(self.n_k):
n_orb = self.n_orbitals[ik,isp] n_orb = self.n_orbitals[ik,ind]
MMat = numpy.identity(n_orb, numpy.complex_) MMat = numpy.identity(n_orb, numpy.complex_)
MMat = self.hopping[ik,isp,0:n_orb,0:n_orb] - (1-2*isp) * self.h_field * MMat MMat = self.hopping[ik,ind,0:n_orb,0:n_orb] - (1-2*isp) * self.h_field * MMat
projmat = self.proj_mat[ik,isp,icrsh,0:dim,0:n_orb] projmat = self.proj_mat[ik,ind,icrsh,0:dim,0:n_orb]
self.Hsumk[icrsh][sp] += self.bz_weights[ik] * numpy.dot( numpy.dot(projmat,MMat), self.Hsumk[icrsh][sp] += self.bz_weights[ik] * numpy.dot( numpy.dot(projmat,MMat),
projmat.conjugate().transpose() ) projmat.conjugate().transpose() )
# symmetrisation: # symmetrisation:
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 range(self.n_corr_shells): for icrsh in range(self.n_corr_shells):
for bl in self.Hsumk[icrsh]: for sp in self.Hsumk[icrsh]:
if self.rot_mat_time_inv[icrsh] == 1: self.Hsumk[icrsh][sp] = self.Hsumk[icrsh][sp].conjugate()
if self.rot_mat_time_inv[icrsh] == 1: self.Hsumk[icrsh][bl] = self.Hsumk[icrsh][bl].conjugate() self.Hsumk[icrsh][sp] = numpy.dot( numpy.dot(self.rot_mat[icrsh].conjugate().transpose(),self.Hsumk[icrsh][sp]) ,
self.Hsumk[icrsh][bl] = numpy.dot( numpy.dot(self.rot_mat[icrsh].conjugate().transpose(),self.Hsumk[icrsh][bl]) ,
self.rot_mat[icrsh] ) self.rot_mat[icrsh] )
# add to matrix: # add to matrix:
for ish in range(self.n_inequiv_shells): for ish in range(self.n_inequiv_shells):
for bl in eff_atlevels[ish]: for sp in eff_atlevels[ish]:
eff_atlevels[ish][bl] += self.Hsumk[self.inequiv_to_corr[ish]][bl] eff_atlevels[ish][sp] += self.Hsumk[self.inequiv_to_corr[ish]][sp]
return eff_atlevels return eff_atlevels
@ -588,6 +596,7 @@ class SumkDFT:
def init_dc(self): def init_dc(self):
""" Initialise the double counting terms to have the correct shape.""" """ Initialise the double counting terms to have the correct shape."""
self.dc_imp = [ {} for icrsh in range(self.n_corr_shells)] self.dc_imp = [ {} for icrsh in range(self.n_corr_shells)]
for icrsh in range(self.n_corr_shells): for icrsh in range(self.n_corr_shells):
dim = self.corr_shells[icrsh]['dim'] dim = self.corr_shells[icrsh]['dim']
@ -672,25 +681,24 @@ class SumkDFT:
mpi.report("DC energy = %s"%self.dc_energ[icrsh]) mpi.report("DC energy = %s"%self.dc_energ[icrsh])
def add_dc(self): def add_dc(self,iw_or_w="iw"):
"""Substracts the double counting term from the impurity self energy.""" """Substracts the double counting term from the impurity self energy."""
# 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] sigma_minus_dc = [s.copy() for s in getattr(self,"Sigma_imp_"+iw_or_w)]
for icrsh in range(self.n_corr_shells): for icrsh in range(self.n_corr_shells):
for bname,gf in sres[icrsh]: for bname,gf in sigma_minus_dc[icrsh]:
# Transform dc_imp to global coordinate system # 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()))
sres[icrsh][bname] -= dccont sigma_minus_dc[icrsh][bname] -= dccont
return sres # list of self energies corrected by DC return sigma_minus_dc
def symm_deg_gf(self,gf_to_symm,orb): def symm_deg_gf(self,gf_to_symm,orb):
"""Symmetrises a GF for the given degenerate shells self.deg_shells""" """Symmetrises a GF for the given degenerate shells self.deg_shells"""
for degsh in self.deg_shells[orb]: for degsh in self.deg_shells[orb]:
#loop over degenerate shells:
ss = gf_to_symm[degsh[0]].copy() ss = gf_to_symm[degsh[0]].copy()
ss.zero() ss.zero()
n_deg = len(degsh) n_deg = len(degsh)
@ -710,10 +718,8 @@ class SumkDFT:
dens = 0.0 dens = 0.0
ikarray = numpy.array(range(self.n_k)) ikarray = numpy.array(range(self.n_k))
for ik in mpi.slice_array(ikarray): for ik in mpi.slice_array(ikarray):
G_latt_iw = self.lattice_gf(ik = ik, mu = mu, iw_or_w = "iw")
S = self.lattice_gf_matsubara(ik = ik, mu = mu) dens += self.bz_weights[ik] * G_latt_iw.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()
@ -748,7 +754,7 @@ class SumkDFT:
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, "calc_density_correction: filename has to be a string!"
ntoi = self.spin_names_to_ind[self.SO] ntoi = self.spin_names_to_ind[self.SO]
spn = self.spin_block_names[self.SO] spn = self.spin_block_names[self.SO]
@ -761,12 +767,12 @@ class SumkDFT:
ikarray = numpy.array(range(self.n_k)) ikarray = numpy.array(range(self.n_k))
for ik in mpi.slice_array(ikarray): for ik in mpi.slice_array(ikarray):
S = self.lattice_gf_matsubara(ik = ik, mu = self.chemical_potential) G_latt_iw = self.lattice_gf(ik = ik, mu = self.chemical_potential, iw_or_w = "iw")
for bname,gf in S: for bname,gf in G_latt_iw:
deltaN[bname][ik] = S[bname].density() deltaN[bname][ik] = G_latt_iw[bname].density()
dens[bname] += self.bz_weights[ik] * S[bname].total_density() dens[bname] += self.bz_weights[ik] * G_latt_iw[bname].total_density()
#put mpi Barrier: # mpi reduce:
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)
@ -783,9 +789,9 @@ class SumkDFT:
# 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 rydberg-1
f.write("%.14f\n"%(S.mesh.beta*self.energy_unit)) f.write("%.14f\n"%(G_latt_iw.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"%(G_latt_iw.mesh.beta*self.energy_unit))
if self.SP == 0: # no spin-polarization if self.SP == 0: # no spin-polarization
@ -847,12 +853,10 @@ class SumkDFT:
return self.chemical_potential return self.chemical_potential
def find_dc(self,orb,guess,dens_mat,dens_req=None,precision=0.01): def calc_dc_for_density(self,orb,dc_init,dens_mat,density=None,precision=0.01):
"""Searches for DC in order to fulfill charge neutrality. """Searches for DC in order to fulfill charge neutrality.
If dens_req is given, then DC is set such that the LOCAL charge of orbital If density is given, then DC is set such that the LOCAL charge of orbital
orb coincides with dens_req.""" orb coincides with the given density."""
mu = self.chemical_potential
def F(dc): def F(dc):
self.calc_dc(dens_mat = dens_mat, U_interact = 0, J_hund = 0, orb = orb, use_dc_value = dc) self.calc_dc(dens_mat = dens_mat, U_interact = 0, J_hund = 0, orb = orb, use_dc_value = dc)
@ -861,24 +865,20 @@ class SumkDFT:
else: else:
return self.extract_G_loc()[orb].total_density() return self.extract_G_loc()[orb].total_density()
if density is None: density = self.density_required - self.charge_below
if dens_req is None: dc = dichotomy.dichotomy(function = F,
density = self.density_required - self.charge_below x_init = dc_init, y_value = density,
else:
density = dens_req
dcnew = dichotomy.dichotomy(function = F,
x_init = guess, y_value = density,
precision_on_y = precision, delta_x = 0.5, precision_on_y = precision, delta_x = 0.5,
max_loops = 100, x_name = "Double Counting", y_name= "Total Density", max_loops = 100, x_name = "Double Counting", y_name= "Total Density",
verbosity = 3)[0] verbosity = 3)[0]
return dcnew return dc
# Check that the density matrix from projectors (DM = P Pdagger) is correct (ie matches DFT)
def check_projectors(self): def check_projectors(self):
"""Calculated the density matrix from projectors (DM = P Pdagger) to check that it is correct and
specifically that it matches DFT."""
dens_mat = [numpy.zeros([self.corr_shells[icrsh]['dim'],self.corr_shells[icrsh]['dim']],numpy.complex_) dens_mat = [numpy.zeros([self.corr_shells[icrsh]['dim'],self.corr_shells[icrsh]['dim']],numpy.complex_)
for icrsh in range(self.n_corr_shells)] for icrsh in range(self.n_corr_shells)]
@ -901,21 +901,19 @@ class SumkDFT:
return dens_mat return dens_mat
# Determine the number of equivalent shells def sorts_of_atoms(self,shells):
def sorts_of_atoms(self,lst):
""" """
This routine should determine the number of sorts in the double list lst Determine the number of inequivalent sorts.
""" """
sortlst = [ lst[i][1] for i in range(len(lst)) ] sortlst = [ shells[i]['sort'] for i in range(len(shells)) ]
sorts = len(set(sortlst)) n_sorts = len(set(sortlst))
return sorts return n_sorts
# Determine the number of atoms def number_of_atoms(self,shells):
def number_of_atoms(self,lst):
""" """
This routine should determine the number of atoms in the double list lst Determine the number of inequivalent atoms.
""" """
atomlst = [ lst[i][0] for i in range(len(lst)) ] atomlst = [ shells[i]['atom'] for i in range(len(shells)) ]
atoms = len(set(atomlst)) n_atoms = len(set(atomlst))
return atoms return n_atoms

210
python/sumk_dft_tools.py
View File

@ -22,7 +22,6 @@
from types import * from types import *
import numpy import numpy
import pytriqs.utility.dichotomy as dichotomy
from pytriqs.gf.local import * from pytriqs.gf.local import *
import pytriqs.utility.mpi as mpi import pytriqs.utility.mpi as mpi
from symmetry import * from symmetry import *
@ -36,7 +35,7 @@ class SumkDFTTools(SumkDFT):
parproj_data = 'dft_parproj_input', symmpar_data = 'dft_symmpar_input', bands_data = 'dft_bands_input', parproj_data = 'dft_parproj_input', symmpar_data = 'dft_symmpar_input', bands_data = 'dft_bands_input',
transp_data = 'dft_transp_input'): transp_data = 'dft_transp_input'):
self.G_upfold_refreq = None #self.G_latt_w = None # DEBUG -- remove later
SumkDFT.__init__(self, hdf_file=hdf_file, h_field=h_field, use_dft_blocks=use_dft_blocks, SumkDFT.__init__(self, hdf_file=hdf_file, h_field=h_field, use_dft_blocks=use_dft_blocks,
dft_data=dft_data, symmcorr_data=symmcorr_data, parproj_data=parproj_data, dft_data=dft_data, symmcorr_data=symmcorr_data, parproj_data=parproj_data,
symmpar_data=symmpar_data, bands_data=bands_data, transp_data=transp_data) symmpar_data=symmpar_data, bands_data=bands_data, transp_data=transp_data)
@ -60,7 +59,7 @@ class SumkDFTTools(SumkDFT):
"""Local <-> Global rotation of a GF block. """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'),"rotloc_all: Give direction 'toLocal' or 'toGlobal'."
gf_rotated = gf_to_rotate.copy() gf_rotated = gf_to_rotate.copy()
@ -82,69 +81,6 @@ class SumkDFTTools(SumkDFT):
return gf_rotated return gf_rotated
def lattice_gf_realfreq(self, ik, mu, broadening, mesh=None, with_Sigma=True):
"""Calculates the lattice Green function on the real frequency axis. If self energy is
present and with_Sigma=True, the mesh is taken from Sigma. Otherwise, the mesh has to be given."""
ntoi = self.spin_names_to_ind[self.SO]
spn = self.spin_block_names[self.SO]
if not hasattr(self,"Sigma_imp"): with_Sigma=False
if with_Sigma:
assert all(type(gf) == GfReFreq for bname,gf in self.Sigma_imp[0]), "Real frequency Sigma needed for lattice_gf_realfreq!"
stmp = self.add_dc()
else:
assert (not mesh is None),"Without Sigma, give the mesh=(om_min,om_max,n_points) for lattice_gf_realfreq!"
if self.G_upfold_refreq is None:
# first setting up of G_upfold_refreq
block_structure = [ range(self.n_orbitals[ik,ntoi[sp]]) for sp in spn ]
gf_struct = [ (spn[isp], block_structure[isp]) for isp in range(self.n_spin_blocks[self.SO]) ]
block_ind_list = [block for block,inner in gf_struct]
if with_Sigma:
glist = lambda : [ GfReFreq(indices = inner, mesh=self.Sigma_imp[0].mesh) for block,inner in gf_struct]
else:
glist = lambda : [ GfReFreq(indices = inner, window=(mesh[0],mesh[1]), n_points=mesh[2]) for block,inner in gf_struct]
self.G_upfold_refreq = BlockGf(name_list = block_ind_list, block_list = glist(),make_copies=False)
self.G_upfold_refreq.zero()
GFsize = [ gf.N1 for bname,gf in self.G_upfold_refreq]
unchangedsize = all( [ self.n_orbitals[ik,ntoi[spn[isp]]] == GFsize[isp]
for isp in range(self.n_spin_blocks[self.SO]) ] )
if not unchangedsize:
block_structure = [ range(self.n_orbitals[ik,ntoi[sp]]) for sp in spn ]
gf_struct = [ (spn[isp], block_structure[isp]) for isp in range(self.n_spin_blocks[self.SO]) ]
block_ind_list = [block for block,inner in gf_struct]
if with_Sigma:
glist = lambda : [ GfReFreq(indices = inner, mesh =self.Sigma_imp[0].mesh) for block,inner in gf_struct]
else:
glist = lambda : [ GfReFreq(indices = inner, window=(mesh[0],mesh[1]), n_points=mesh[2]) for block,inner in gf_struct]
self.G_upfold_refreq = BlockGf(name_list = block_ind_list, block_list = glist(),make_copies=False)
self.G_upfold_refreq.zero()
idmat = [numpy.identity(self.n_orbitals[ik,ntoi[sp]],numpy.complex_) for sp in spn]
self.G_upfold_refreq << Omega + 1j*broadening
M = copy.deepcopy(idmat)
for isp in range(self.n_spin_blocks[self.SO]):
ind = ntoi[spn[isp]]
n_orb = self.n_orbitals[ik,ind]
M[isp] = self.hopping[ik,ind,0:n_orb,0:n_orb] - (idmat[isp]*mu) - (idmat[isp] * self.h_field * (1-2*isp))
self.G_upfold_refreq -= M
if with_Sigma:
tmp = self.G_upfold_refreq.copy() # init temporary storage
for icrsh in range(self.n_corr_shells):
for bname,gf in tmp: tmp[bname] << self.upfold(ik,icrsh,bname,stmp[icrsh][bname],gf)
self.G_upfold_refreq -= tmp # adding to the upfolded GF
self.G_upfold_refreq.invert()
return self.G_upfold_refreq
def check_input_dos(self, om_min, om_max, n_om, beta=10, broadening=0.01): def check_input_dos(self, om_min, om_max, n_om, beta=10, broadening=0.01):
@ -153,16 +89,16 @@ class SumkDFTTools(SumkDFT):
for i in range(n_om): om_mesh[i] = om_min + delta_om * i for i in range(n_om): om_mesh[i] = om_min + delta_om * i
DOS = {} DOS = {}
for bn in self.spin_block_names[self.SO]: for sp in self.spin_block_names[self.SO]:
DOS[bn] = numpy.zeros([n_om],numpy.float_) DOS[sp] = numpy.zeros([n_om],numpy.float_)
DOSproj = [ {} for ish in range(self.n_inequiv_shells) ] DOSproj = [ {} for ish in range(self.n_inequiv_shells) ]
DOSproj_orb = [ {} for ish in range(self.n_inequiv_shells) ] DOSproj_orb = [ {} for ish in range(self.n_inequiv_shells) ]
for ish in range(self.n_inequiv_shells): for ish in range(self.n_inequiv_shells):
for bn in self.spin_block_names[self.corr_shells[self.inequiv_to_corr[ish]]['SO']]: for sp in self.spin_block_names[self.corr_shells[self.inequiv_to_corr[ish]]['SO']]:
dim = self.corr_shells[self.inequiv_to_corr[ish]]['dim'] dim = self.corr_shells[self.inequiv_to_corr[ish]]['dim']
DOSproj[ish][bn] = numpy.zeros([n_om],numpy.float_) DOSproj[ish][sp] = numpy.zeros([n_om],numpy.float_)
DOSproj_orb[ish][bn] = numpy.zeros([n_om,dim,dim],numpy.float_) DOSproj_orb[ish][sp] = numpy.zeros([n_om,dim,dim],numpy.float_)
# init: # init:
Gloc = [] Gloc = []
@ -174,18 +110,18 @@ class SumkDFTTools(SumkDFT):
for ik in range(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_latt_w=self.lattice_gf(ik=ik,mu=self.chemical_potential,iw_or_w="w",broadening=broadening,mesh=(om_min,om_max,n_om),with_Sigma=False)
G_upfold *= self.bz_weights[ik] G_latt_w *= self.bz_weights[ik]
# non-projected DOS # non-projected DOS
for iom in range(n_om): for iom in range(n_om):
for bname,gf in G_upfold: for bname,gf in G_latt_w:
asd = gf.data[iom,:,:].imag.trace()/(-3.1415926535) asd = gf.data[iom,:,:].imag.trace()/(-3.1415926535)
DOS[bname][iom] += asd DOS[bname][iom] += asd
for icrsh in range(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_latt_w[bname],gf) # downfolding G
Gloc[icrsh] += tmp Gloc[icrsh] += tmp
@ -205,21 +141,21 @@ class SumkDFTTools(SumkDFT):
# output: # output:
if mpi.is_master_node(): if mpi.is_master_node():
for bn in self.spin_block_names[self.SO]: for sp in self.spin_block_names[self.SO]:
f=open('DOS%s.dat'%bn, 'w') f=open('DOS%s.dat'%sp, '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[sp][i]))
f.close() f.close()
for ish in range(self.n_inequiv_shells): for ish in range(self.n_inequiv_shells):
f=open('DOS%s_proj%s.dat'%(bn,ish),'w') f=open('DOS%s_proj%s.dat'%(sp,ish),'w')
for i in range(n_om): f.write("%s %s\n"%(om_mesh[i],DOSproj[ish][bn][i])) for i in range(n_om): f.write("%s %s\n"%(om_mesh[i],DOSproj[ish][sp][i]))
f.close() f.close()
for i in range(self.corr_shells[self.inequiv_to_corr[ish]]['dim']): for i in range(self.corr_shells[self.inequiv_to_corr[ish]]['dim']):
for j in range(i,self.corr_shells[self.inequiv_to_corr[ish]]['dim']): for j in range(i,self.corr_shells[self.inequiv_to_corr[ish]]['dim']):
Fname = 'DOS'+bn+'_proj'+str(ish)+'_'+str(i)+'_'+str(j)+'.dat' Fname = 'DOS'+sp+'_proj'+str(ish)+'_'+str(i)+'_'+str(j)+'.dat'
f=open(Fname,'w') f=open(Fname,'w')
for iom in range(n_om): f.write("%s %s\n"%(om_mesh[iom],DOSproj_orb[ish][bn][iom,i,j])) for iom in range(n_om): f.write("%s %s\n"%(om_mesh[iom],DOSproj_orb[ish][sp][iom,i,j]))
f.close() f.close()
@ -239,7 +175,7 @@ class SumkDFTTools(SumkDFT):
def dos_partial(self,broadening=0.01): def dos_partial(self,broadening=0.01):
"""calculates the orbitally-resolved DOS""" """calculates the orbitally-resolved DOS"""
assert hasattr(self,"Sigma_imp"), "Set Sigma First!!" assert hasattr(self,"Sigma_imp_w"), "dos_partial: Set Sigma_imp_w first."
value_read = self.read_parproj_input_from_hdf() value_read = self.read_parproj_input_from_hdf()
if not value_read: return value_read if not value_read: return value_read
@ -248,41 +184,41 @@ class SumkDFTTools(SumkDFT):
mu = self.chemical_potential mu = self.chemical_potential
gf_struct_proj = [ [ (sp, range(self.shells[i]['dim'])) for sp in self.spin_block_names[self.SO] ] for i in range(self.n_shells) ] gf_struct_proj = [ [ (sp, range(self.shells[i]['dim'])) for sp in self.spin_block_names[self.SO] ] for i in range(self.n_shells) ]
Gproj = [BlockGf(name_block_generator = [ (block,GfReFreq(indices = inner, mesh = self.Sigma_imp[0].mesh)) for block,inner in gf_struct_proj[ish] ], make_copies = False ) Gproj = [BlockGf(name_block_generator = [ (block,GfReFreq(indices = inner, mesh = self.Sigma_imp_w[0].mesh))
for ish in range(self.n_shells)] for block,inner in gf_struct_proj[ish] ], make_copies = False ) for ish in range(self.n_shells)]
for ish in range(self.n_shells): Gproj[ish].zero() for ish in range(self.n_shells): Gproj[ish].zero()
Msh = [x.real for x in self.Sigma_imp[0].mesh] mesh = [x.real for x in self.Sigma_imp_w[0].mesh]
n_om = len(Msh) n_om = len(mesh)
DOS = {} DOS = {}
for bn in self.spin_block_names[self.SO]: for sp in self.spin_block_names[self.SO]:
DOS[bn] = numpy.zeros([n_om],numpy.float_) DOS[sp] = 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.spin_block_names[self.SO]: for sp in self.spin_block_names[self.SO]:
dim = self.shells[ish]['dim'] dim = self.shells[ish]['dim']
DOSproj[ish][bn] = numpy.zeros([n_om],numpy.float_) DOSproj[ish][sp] = numpy.zeros([n_om],numpy.float_)
DOSproj_orb[ish][bn] = numpy.zeros([n_om,dim,dim],numpy.float_) DOSproj_orb[ish][sp] = numpy.zeros([n_om,dim,dim],numpy.float_)
ikarray=numpy.array(range(self.n_k)) ikarray=numpy.array(range(self.n_k))
for ik in mpi.slice_array(ikarray): for ik in mpi.slice_array(ikarray):
S = self.lattice_gf_realfreq(ik=ik,mu=mu,broadening=broadening) G_latt_w = self.lattice_gf(ik=ik,mu=mu,iw_or_w="w",broadening=broadening)
S *= self.bz_weights[ik] G_latt_w *= self.bz_weights[ik]
# non-projected DOS # non-projected DOS
for iom in range(n_om): for iom in range(n_om):
for bname,gf in S: DOS[bname][iom] += gf.data[iom,:,:].imag.trace()/(-3.1415926535) for bname,gf in G_latt_w: DOS[bname][iom] += gf.data[iom,:,:].imag.trace()/(-3.1415926535)
#projected DOS: #projected DOS:
for ish in range(self.n_shells): for ish in range(self.n_shells):
tmp = Gproj[ish].copy() tmp = Gproj[ish].copy()
for ir in range(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,G_latt_w[bname],gf)
Gproj[ish] += tmp Gproj[ish] += tmp
# collect data from mpi: # collect data from mpi:
@ -307,22 +243,22 @@ class SumkDFTTools(SumkDFT):
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 sp in self.spin_block_names[self.SO]:
f=open('./DOScorr%s.dat'%bn, 'w') f=open('./DOScorr%s.dat'%sp, '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"%(mesh[i],DOS[sp][i]))
f.close() f.close()
# partial # partial
for ish in range(self.n_shells): for ish in range(self.n_shells):
f=open('DOScorr%s_proj%s.dat'%(bn,ish),'w') f=open('DOScorr%s_proj%s.dat'%(sp,ish),'w')
for i in range(n_om): f.write("%s %s\n"%(Msh[i],DOSproj[ish][bn][i])) for i in range(n_om): f.write("%s %s\n"%(mesh[i],DOSproj[ish][sp][i]))
f.close() f.close()
for i in range(self.shells[ish]['dim']): for i in range(self.shells[ish]['dim']):
for j in range(i,self.shells[ish]['dim']): for j in range(i,self.shells[ish]['dim']):
Fname = './DOScorr'+bn+'_proj'+str(ish)+'_'+str(i)+'_'+str(j)+'.dat' Fname = './DOScorr'+sp+'_proj'+str(ish)+'_'+str(i)+'_'+str(j)+'.dat'
f=open(Fname,'w') f=open(Fname,'w')
for iom in range(n_om): f.write("%s %s\n"%(Msh[iom],DOSproj_orb[ish][bn][iom,i,j])) for iom in range(n_om): f.write("%s %s\n"%(mesh[iom],DOSproj_orb[ish][sp][iom,i,j]))
f.close() f.close()
@ -332,7 +268,7 @@ class SumkDFTTools(SumkDFT):
""" Calculates the correlated band structure with a real-frequency self energy. """ Calculates the correlated band structure with a real-frequency self energy.
ATTENTION: Many things from the original input file are overwritten!!!""" ATTENTION: Many things from the original input file are overwritten!!!"""
assert hasattr(self,"Sigma_imp"), "Set Sigma First!!" assert hasattr(self,"Sigma_imp_w"), "spaghettis: Set Sigma_imp_w first."
things_to_read = ['n_k','n_orbitals','proj_mat','hopping','n_parproj','proj_mat_pc'] things_to_read = ['n_k','n_orbitals','proj_mat','hopping','n_parproj','proj_mat_pc']
value_read = self.read_input_from_hdf(subgrp=self.bands_data,things_to_read=things_to_read) value_read = self.read_input_from_hdf(subgrp=self.bands_data,things_to_read=things_to_read)
if not value_read: return value_read if not value_read: return value_read
@ -370,12 +306,12 @@ class SumkDFTTools(SumkDFT):
spn = self.spin_block_names[self.SO] spn = self.spin_block_names[self.SO]
# init DOS: # init DOS:
M = [x.real for x in self.Sigma_imp[0].mesh] mesh = [x.real for x in self.Sigma_imp_w[0].mesh]
n_om = len(M) n_om = len(mesh)
if plot_range is None: if plot_range is None:
om_minplot = M[0]-0.001 om_minplot = mesh[0]-0.001
om_maxplot = M[n_om-1] + 0.001 om_maxplot = mesh[n_om-1] + 0.001
else: else:
om_minplot = plot_range[0] om_minplot = plot_range[0]
om_maxplot = plot_range[1] om_maxplot = plot_range[1]
@ -395,20 +331,20 @@ class SumkDFTTools(SumkDFT):
if not ishell is None: if not ishell is None:
GFStruct_proj = [ (sp, range(self.shells[ishell]['dim'])) for sp in spn ] GFStruct_proj = [ (sp, range(self.shells[ishell]['dim'])) for sp in spn ]
Gproj = BlockGf(name_block_generator = [ (block,GfReFreq(indices = inner, mesh = self.Sigma_imp[0].mesh)) for block,inner in GFStruct_proj ], make_copies = False) Gproj = BlockGf(name_block_generator = [ (block,GfReFreq(indices = inner, mesh = self.Sigma_imp_w[0].mesh)) for block,inner in GFStruct_proj ], make_copies = False)
Gproj.zero() Gproj.zero()
for ik in range(self.n_k): for ik in range(self.n_k):
S = self.lattice_gf_realfreq(ik=ik,mu=mu,broadening=broadening) G_latt_w = self.lattice_gf(ik=ik,mu=mu,iw_or_w="w",broadening=broadening)
if ishell is None: 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 (mesh[iom] > om_minplot) and (mesh[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 G_latt_w: Akw[bname][ik,0] += gf.data[iom,:,:].imag.trace()/(-3.1415926535) * (mesh[1]-mesh[0])
else: else:
for bname,gf in S: Akw[bname][ik,iom] += gf.data[iom,:,:].imag.trace()/(-3.1415926535) for bname,gf in G_latt_w: Akw[bname][ik,iom] += gf.data[iom,:,:].imag.trace()/(-3.1415926535)
Akw[bname][ik,iom] += ik*shift # shift Akw for plotting in xmgrace -- REMOVE Akw[bname][ik,iom] += ik*shift # shift Akw for plotting in xmgrace -- REMOVE
@ -417,7 +353,7 @@ class SumkDFTTools(SumkDFT):
Gproj.zero() Gproj.zero()
tmp = Gproj.copy() tmp = Gproj.copy()
for ir in range(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,G_latt_w[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
@ -427,7 +363,7 @@ class SumkDFTTools(SumkDFT):
# 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 (mesh[iom] > om_minplot) and (mesh[iom] < om_maxplot):
for ish in range(self.shells[ishell]['dim']): for ish in range(self.shells[ishell]['dim']):
for ibn in spn: for ibn in spn:
Akw[ibn][ish,ik,iom] = Gproj[ibn].data[iom,ish,ish].imag/(-3.1415926535) Akw[ibn][ish,ik,iom] = Gproj[ibn].data[iom,ish,ish].imag/(-3.1415926535)
@ -458,13 +394,13 @@ class SumkDFTTools(SumkDFT):
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 (mesh[iom] > om_minplot) and (mesh[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'%(mesh[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,mesh[iom],Akw[ibn][ik,iom]))
f.write('\n') f.write('\n')
@ -482,13 +418,13 @@ class SumkDFTTools(SumkDFT):
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 (mesh[iom] > om_minplot) and (mesh[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'%(mesh[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,mesh[iom],Akw[ibn][ish,ik,iom]))
f.write('\n') f.write('\n')
@ -511,10 +447,10 @@ class SumkDFTTools(SumkDFT):
mu = self.chemical_potential mu = self.chemical_potential
GFStruct_proj = [ [ (sp, range(self.shells[i]['dim'])) for sp in spn ] for i in range(self.n_shells) ] GFStruct_proj = [ [ (sp, range(self.shells[i]['dim'])) for sp in spn ] for i in range(self.n_shells) ]
if hasattr(self,"Sigma_imp"): if hasattr(self,"Sigma_imp_iw"):
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_iw[0].mesh)) for block,inner in GFStruct_proj[ish] ], make_copies = False)
for ish in range(self.n_shells)] for ish in range(self.n_shells)]
beta = self.Sigma_imp[0].mesh.beta beta = self.Sigma_imp_iw[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 range(self.n_shells)] for ish in range(self.n_shells)]
@ -524,13 +460,13 @@ class SumkDFTTools(SumkDFT):
ikarray=numpy.array(range(self.n_k)) ikarray=numpy.array(range(self.n_k))
for ik in mpi.slice_array(ikarray): for ik in mpi.slice_array(ikarray):
S = self.lattice_gf_matsubara(ik=ik,mu=mu,beta=beta) G_latt_iw = self.lattice_gf(ik=ik,mu=mu,iw_or_w="iw",beta=beta)
S *= self.bz_weights[ik] G_latt_iw *= self.bz_weights[ik]
for ish in range(self.n_shells): for ish in range(self.n_shells):
tmp = Gproj[ish].copy() tmp = Gproj[ish].copy()
for ir in range(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,G_latt_iw[bname],gf)
Gproj[ish] += tmp Gproj[ish] += tmp
#collect data from mpi: #collect data from mpi:
@ -619,7 +555,7 @@ class SumkDFTTools(SumkDFT):
# Define mesh for Greens function and the used energy range # Define mesh for Greens function and the used energy range
if (with_Sigma == False): if (with_Sigma == False):
self.omega = numpy.array([round(x.real,12) for x in self.Sigma_imp[0].mesh]) self.omega = numpy.array([round(x.real,12) for x in self.Sigma_imp_w[0].mesh])
mu = self.chemical_potential mu = self.chemical_potential
n_om = len(self.omega) n_om = len(self.omega)
print "Using omega mesh provided by Sigma." print "Using omega mesh provided by Sigma."
@ -632,15 +568,15 @@ class SumkDFTTools(SumkDFT):
# Truncate Sigma to given omega window # Truncate Sigma to given omega window
for icrsh in range(self.n_corr_shells): for icrsh in range(self.n_corr_shells):
Sigma_save = self.Sigma_imp[icrsh].copy() Sigma_save = self.Sigma_imp_w[icrsh].copy()
spn = self.spin_block_names[self.corr_shells[icrsh]['SO']] spn = self.spin_block_names[self.corr_shells[icrsh]['SO']]
glist = lambda : [ GfReFreq(indices = inner, window=(self.omega[0], self.omega[-1]),n_points=n_om) for block, inner in self.gf_struct_sumk[icrsh]] glist = lambda : [ GfReFreq(indices = inner, window=(self.omega[0], self.omega[-1]),n_points=n_om) for block, inner in self.gf_struct_sumk[icrsh]]
self.Sigma_imp[icrsh] = BlockGf(name_list = spn, block_list = glist(),make_copies=False) self.Sigma_imp_w[icrsh] = BlockGf(name_list = spn, block_list = glist(),make_copies=False)
for i,g in self.Sigma_imp[icrsh]: for i,g in self.Sigma_imp_w[icrsh]:
for iL in g.indices: for iL in g.indices:
for iR in g.indices: for iR in g.indices:
for iom in xrange(n_om): for iom in xrange(n_om):
g.data[iom,iL,iR]= Sigma_save[i].data[ioffset+iom,iL,iR] g.data[iom,iL,iR] = Sigma_save[i].data[ioffset+iom,iL,iR] # FIXME IS THIS CLEAN? MANIPULATING data DIRECTLY?
else: else:
assert n_om is not None, "Number of omega points (n_om) needed!" assert n_om is not None, "Number of omega points (n_om) needed!"
@ -684,8 +620,8 @@ class SumkDFTTools(SumkDFT):
ikarray = numpy.array(range(self.n_k)) ikarray = numpy.array(range(self.n_k))
for ik in mpi.slice_array(ikarray): for ik in mpi.slice_array(ikarray):
unchangesize = all([ self.n_orbitals[ik][isp] == mupat[isp].shape[0] for isp in range(n_inequiv_spin_blocks)]) unchangedsize = all([ self.n_orbitals[ik][isp] == mupat[isp].shape[0] for isp in range(n_inequiv_spin_blocks)])
if (not unchangesize): if not unchangedsize:
# recontruct green functions. # recontruct green functions.
S = BlockGf(name_block_generator=[(bln[isp], GfReFreq(indices=range(self.n_orbitals[ik][isp]), window = (self.omega[0], self.omega[-1]), n_points = n_om)) S = BlockGf(name_block_generator=[(bln[isp], GfReFreq(indices=range(self.n_orbitals[ik][isp]), window = (self.omega[0], self.omega[-1]), n_points = n_om))
for isp in range(n_inequiv_spin_blocks) ], make_copies=False) for isp in range(n_inequiv_spin_blocks) ], make_copies=False)
@ -707,10 +643,10 @@ class SumkDFTTools(SumkDFT):
if (with_Sigma == False): if (with_Sigma == False):
tmp = S.copy() # init temporary storage tmp = S.copy() # init temporary storage
# form self energy from impurity self energy and double counting term. # form self energy from impurity self energy and double counting term.
stmp = self.add_dc() sigma_minus_dc = self.add_dc(iw_or_w="w")
## substract self energy # substract self energy
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 sig, gf in tmp: tmp[sig] << self.upfold(ik, icrsh, sig, sigma_minus_dc[icrsh][sig], gf)
S -= tmp S -= tmp
S.invert() S.invert()

6
python/symmetry.py
View File

@ -39,7 +39,7 @@ class Symmetry:
SP: Flag for spin polarisation. SP: Flag for spin polarisation.
""" """
assert type(hdf_file) == StringType, "hdf_file must be a filename" assert type(hdf_file) == StringType, "Symmetry: hdf_file must be a filename."
self.hdf_file = hdf_file self.hdf_file = hdf_file
things_to_read = ['n_symm','n_atoms','perm','orbits','SO','SP','time_inv','mat','mat_tinv'] things_to_read = ['n_symm','n_atoms','perm','orbits','SO','SP','time_inv','mat','mat_tinv']
for it in things_to_read: setattr(self,it,0) for it in things_to_read: setattr(self,it,0)
@ -73,8 +73,8 @@ class Symmetry:
def symmetrize(self,obj): def symmetrize(self,obj):
assert isinstance(obj,list), "symmetry: obj has to be a list of objects." assert isinstance(obj,list), "symmetrize: obj has to be a list of objects."
assert len(obj) == self.n_orbits, "symmetry: obj has to be a list of the same length as defined in the init." assert len(obj) == self.n_orbits, "symmetrize: obj has to be a list of the same length as defined in the init."
if isinstance(obj[0],BlockGf): if isinstance(obj[0],BlockGf):
symm_obj = [ obj[i].copy() for i in range(len(obj)) ] # here the result is stored, it is a BlockGf! symm_obj = [ obj[i].copy() for i in range(len(obj)) ] # here the result is stored, it is a BlockGf!