LCPQ/dft_tools
3
0
Fork 0
mirror of https://github.com/triqs/dft_tools synced 2024-10-04 23:36:05 +02:00
Code Issues Releases Wiki Activity

Fix bugs in sumk_lda_tools.py

Browse Source 
This is again a merge of Leonid's and Markus' changes.
It must be checked! I had to make some choices!

  modified:   python/sumk_lda_tools.py
This commit is contained in:
Michel Ferrero 2014-01-21 16:47:45 +01:00
parent 2d5cf0abc7
commit 8db8a52eaa
1 changed files with 102 additions and 50 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

152
python/sumk_lda_tools.py
View File

@ -42,8 +42,8 @@ def read_fortran_file (filename):
import os.path import os.path
if not(os.path.exists(filename)) : raise IOError, "File %s does not exists"%filename if not(os.path.exists(filename)) : raise IOError, "File %s does not exists"%filename
for line in open(filename,'r') : for line in open(filename,'r') :
for x in line.replace('D','E').split() : for x in line.replace('D','E').split() :
yield string.atof(x) yield string.atof(x)
class SumkLDATools(SumkLDA): class SumkLDATools(SumkLDA):
@ -99,22 +99,21 @@ class SumkLDATools(SumkLDA):
return gf_rotated return gf_rotated
def lattice_gf_realfreq(self, ik, mu, broadening, mesh=None, beta=40, with_Sigma=True): 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 """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.names_to_ind[self.SO]
bln = self.blocnames[self.SO] bln = self.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 type(self.Sigma_imp[0]) == GfReFreq, "Real frequency Sigma needed for lattice_gf_realfreq!" assert self.Sigma_imp[0].note == 'ReFreq', "Real frequency Sigma needed for lattice_gf_realfreq!"
beta = self.Sigma_imp[0].mesh.beta
stmp = self.add_dc() stmp = self.add_dc()
else: else:
assert (not (mesh is None)),"Without Sigma, give the mesh 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.Gupf_refreq is None): if (self.Gupf_refreq is None):
# first setting up of Gupf_refreq # first setting up of Gupf_refreq
BS = [ range(self.n_orbitals[ik,ntoi[ib]]) for ib in bln ] BS = [ range(self.n_orbitals[ik,ntoi[ib]]) for ib in bln ]
gf_struct = [ (bln[ib], BS[ib]) for ib in range(self.n_spin_blocks_gf[self.SO]) ] gf_struct = [ (bln[ib], BS[ib]) for ib in range(self.n_spin_blocks_gf[self.SO]) ]
@ -122,7 +121,7 @@ class SumkLDATools(SumkLDA):
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 = al, mesh =self.Sigma_imp[0].mesh) for a,al in gf_struct]
else: else:
glist = lambda : [ GfReFreq(indices = al, beta = beta, mesh_array = mesh) for a,al in gf_struct] glist = lambda : [ GfReFreq(indices = al, window=(mesh[0],mesh[1]),n_points=mesh[2]) for a,al in gf_struct]
self.Gupf_refreq = BlockGf(name_list = a_list, block_list = glist(),make_copies=False) self.Gupf_refreq = BlockGf(name_list = a_list, block_list = glist(),make_copies=False)
self.Gupf_refreq.zero() self.Gupf_refreq.zero()
@ -137,7 +136,7 @@ class SumkLDATools(SumkLDA):
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 = al, mesh =self.Sigma_imp[0].mesh) for a,al in gf_struct]
else: else:
glist = lambda : [ GfReFreq(indices = al, beta = beta, mesh_array = mesh) for a,al in gf_struct] glist = lambda : [ GfReFreq(indices = al, window=(mesh[0],mesh[1]),n_points=mesh[2]) for a,al in gf_struct]
self.Gupf_refreq = BlockGf(name_list = a_list, block_list = glist(),make_copies=False) self.Gupf_refreq = BlockGf(name_list = a_list, block_list = glist(),make_copies=False)
self.Gupf_refreq.zero() self.Gupf_refreq.zero()
@ -165,9 +164,10 @@ class SumkLDATools(SumkLDA):
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):
delta_om = (om_max-om_min)/(n_om-1) delta_om = (om_max-om_min)/(n_om-1)
mesh = numpy.zeros([n_om],numpy.float_) om_mesh = numpy.zeros([n_om],numpy.float_)
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.block_names[self.SO]:
@ -179,22 +179,20 @@ class SumkLDATools(SumkLDA):
for bn in self.block_names[self.corr_shells[self.invshellmap[icrsh]][4]]: for bn in self.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([dl,dl,n_om],numpy.float_) DOSproj_orb[icrsh][bn] = numpy.zeros([n_om,dl,dl],numpy.float_)
for i in range(n_om): mesh[i] = om_min + delta_om * i
# 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 = [a for a,al 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 = al, beta = beta, mesh_array = mesh) for a,al 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]]
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
for ik in xrange(self.n_k): for ik in xrange(self.n_k):
Gupf=self.lattice_gf_realfreq(ik=ik,mu=self.chemical_potential,broadening=broadening,beta=beta,mesh=mesh,with_Sigma=False) Gupf=self.lattice_gf_realfreq(ik=ik,mu=self.chemical_potential,broadening=broadening,mesh=(om_min,om_max,n_om),with_Sigma=False)
Gupf *= self.bz_weights[ik] Gupf *= self.bz_weights[ik]
# non-projected DOS # non-projected DOS
@ -227,19 +225,19 @@ class SumkLDATools(SumkLDA):
if (mpi.is_master_node()): if (mpi.is_master_node()):
for bn in self.block_names[self.SO]: for bn in self.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"%(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()
for ish in range(self.n_inequiv_corr_shells): for ish in range(self.n_inequiv_corr_shells):
f=open('DOS%s_proj%s.dat'%(bn,ish),'w') f=open('DOS%s_proj%s.dat'%(bn,ish),'w')
for i in range(n_om): f.write("%s %s\n"%(mesh[i],DOSproj[ish][bn][i])) for i in range(n_om): f.write("%s %s\n"%(om_mesh[i],DOSproj[ish][bn][i]))
f.close() f.close()
for i in range(self.corr_shells[self.invshellmap[ish]][3]): for i in range(self.corr_shells[self.invshellmap[ish]][3]):
for j in range(i,self.corr_shells[self.invshellmap[ish]][3]): for j in range(i,self.corr_shells[self.invshellmap[ish]][3]):
Fname = 'DOS'+bn+'_proj'+str(ish)+'_'+str(i)+'_'+str(j)+'.dat' Fname = 'DOS'+bn+'_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"%(mesh[iom],DOSproj_orb[ish][bn][i,j,iom])) for iom in range(n_om): f.write("%s %s\n"%(om_mesh[iom],DOSproj_orb[ish][bn][iom,i,j]))
f.close() f.close()
@ -274,7 +272,7 @@ class SumkLDATools(SumkLDA):
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()
Msh = [x for x in self.Sigma_imp[0].mesh] Msh = [x.real for x in self.Sigma_imp[0].mesh]
n_om = len(Msh) n_om = len(Msh)
DOS = {} DOS = {}
@ -287,7 +285,7 @@ class SumkLDATools(SumkLDA):
for bn in self.block_names[self.SO]: for bn in self.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([dl,dl,n_om],numpy.float_) DOSproj_orb[ish][bn] = numpy.zeros([n_om,dl,dl],numpy.float_)
ikarray=numpy.array(range(self.n_k)) ikarray=numpy.array(range(self.n_k))
@ -325,7 +323,7 @@ class SumkLDATools(SumkLDA):
for sig,gf in Gproj[ish]: for sig,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][sig][iom] += gf.data[iom,:,:].imag.trace()/(-3.1415926535)
DOSproj_orb[ish][sig][:,:,:] += gf.data[:,:,:].imag / (-3.1415926535) DOSproj_orb[ish][sig][:,:,:] += gf.data[:,:,:].imag / (-3.1415926535)
if (mpi.is_master_node()): if (mpi.is_master_node()):
# output to files # output to files
@ -344,7 +342,7 @@ class SumkLDATools(SumkLDA):
for j in range(i,self.shells[ish][3]): for j in range(i,self.shells[ish][3]):
Fname = './DOScorr'+bn+'_proj'+str(ish)+'_'+str(i)+'_'+str(j)+'.dat' Fname = './DOScorr'+bn+'_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][i,j,iom])) for iom in range(n_om): f.write("%s %s\n"%(Msh[iom],DOSproj_orb[ish][bn][iom,i,j]))
f.close() f.close()
@ -514,38 +512,93 @@ class SumkLDATools(SumkLDA):
f.close() f.close()
def constr_Sigma_ME(self, filename, beta, n_om, orb = 0):
"""Uses Data from files to construct a GF object on the real axis."""
def constr_Sigma_real_axis(self, filename, hdf=True, hdf_dataset='SigmaReFreq',n_om=0,orb=0, tol_mesh=1e-6):
"""Uses Data from files to construct Sigma (or GF) on the real axis."""
if not hdf:
# read sigma from text files
#first get the mesh out of one of the files:
if (len(self.gf_struct_solver[orb][0][1])==1):
Fname = filename+'_'+self.gf_struct_solver[orb][0][0]+'.dat'
else:
Fname = filename+'_'+self.gf_struct_solver[orb][0][0]+'/'+str(self.gf_struct_solver[orb][0][1][0])+'_'+str(self.gf_struct_solver[orb][0][1][0])+'.dat'
R = read_fortran_file(Fname)
mesh = numpy.zeros([n_om],numpy.float_)
try:
for i in xrange(n_om):
mesh[i] = R.next()
sk = R.next()
sk = R.next()
except StopIteration : # a more explicit error if the file is corrupted.
raise "SumkLDA.read_Sigma_ME : reading mesh failed!"
R.close()
# check whether the mesh is uniform
bin = (mesh[n_om-1]-mesh[0])/(n_om-1)
for i in xrange(n_om):
assert abs(i*bin+mesh[0]-mesh[i]) < tol_mesh, 'constr_Sigma_ME: real-axis mesh is non-uniform!'
# construct Sigma
a_list = [a for a,al in self.gf_struct_solver[orb]]
glist = lambda : [ GfReFreq(indices = al, window=(mesh[0],mesh[n_om-1]),n_points=n_om) for a,al in self.gf_struct_solver[orb]]
SigmaME = BlockGf(name_list = a_list, block_list = glist(),make_copies=False)
#read Sigma
#first get the mesh out of one of the files: for i,g in SigmaME:
if (len(self.gf_struct_solver[orb][0][1])==1): mesh=[w for w in g.mesh]
Fname = filename+'_'+self.gf_struct_solver[orb][0][0]+'.dat' for iL in g.indices:
for iR in g.indices:
if (len(g.indices) == 1):
Fname = filename+'_%s'%(i)+'.dat'
else:
Fname = 'SigmaME_'+'%s'%(i)+'_%s'%(iL)+'_%s'%(iR)+'.dat'
R = read_fortran_file(Fname)
try:
for iom in xrange(n_om):
sk = R.next()
rsig = R.next()
isig = R.next()
g.data[iom,iL,iR]=rsig+1j*isig
except StopIteration : # a more explicit error if the file is corrupted.
raise "SumkLDA.read_Sigma_ME : reading Sigma from file failed!"
R.close()
else: else:
Fname = filename+'_'+self.gf_struct_solver[orb][0][0]+'/'+str(self.gf_struct_solver[orb][0][1][0])+'_'+str(self.gf_struct_solver[orb][0][1][0])+'.dat'
R = read_fortran_file(Fname) # read sigma from hdf
mesh = numpy.zeros([n_om],numpy.float_) omega_min=0.0
try: omega_max=0.0
for i in xrange(n_om): n_om=0
mesh[i] = R.next() if (mpi.is_master_node()):
sk = R.next() ar = HDFArchive(filename)
sk = R.next() SigmaME = ar[hdf_dataset]
del ar
except StopIteration : # a more explicit error if the file is corrupted. # we need some parameters to construct Sigma on other nodes
raise "SumkLDA.read_Sigma_ME : reading file failed!" omega_min=SigmaME.mesh.omega_min
R.close() omega_max=SigmaME.mesh.omega_max
n_om=len(SigmaME.mesh)
omega_min=mpi.bcast(omega_min)
omega_max=mpi.bcast(omega_max)
n_om=mpi.bcast(n_om)
mpi.barrier()
# construct Sigma on other nodes
if (not mpi.is_master_node()):
a_list = [a for a,al in self.gf_struct_solver[orb]]
glist = lambda : [ GfReFreq(indices = al, window=(omega_min,omega_max),n_points=n_om) for a,al in self.gf_struct_solver[orb]]
SigmaME = BlockGf(name_list = a_list, block_list = glist(),make_copies=False)
# pass SigmaME to other nodes
SigmaME = mpi.bcast(SigmaME)
mpi.barrier()
# now initialize the GF with the mesh SigmaME.note='ReFreq'
a_list = [a for a,al in self.gf_struct_solver[orb]]
glist = lambda : [ GfReFreq(indices = al, beta = beta, mesh_array = mesh) for a,al in self.gf_struct_solver[orb] ]
SigmaME = BlockGf(name_list = a_list, block_list = glist(),make_copies=False)
SigmaME.load(filename)
return SigmaME return SigmaME
def partial_charges(self,beta=40): def partial_charges(self,beta=40):
@ -622,4 +675,3 @@ class SumkLDATools(SumkLDA):
return dens_mat return dens_mat