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

[transport] Case Sigma=0 included (LDA_only)

Browse Source 
Some other minor corrections.
This commit is contained in:
Manuel Zingl 2014-11-03 19:15:05 +01:00 committed by Priyanka Seth
parent 2ee744854e
commit cc82ba2d5a
1 changed files with 32 additions and 33 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

65
python/SumK_LDA_Transport.py
View File

@ -436,10 +436,7 @@ class TransportEtOptic(SumkLDATools):
bln = self.block_names[self.SO] bln = self.block_names[self.SO]
ntoi = self.names_to_ind[self.SO] ntoi = self.names_to_ind[self.SO]
S = BlockGf(name_block_generator=[(bln[isp], GfReFreq(indices=range(self.n_orbitals[ik][isp]), S = BlockGf(name_block_generator=[(bln[isp], GfReFreq(indices=range(self.n_orbitals[ik][isp]), mesh=self.Sigma_imp[0].mesh)) for isp in range(self.Nspinblocs) ], make_copies=False)
mesh=self.Sigma_imp[0].mesh))
for isp in range(self.Nspinblocs) ],
make_copies=False)
mupat = [numpy.identity(self.n_orbitals[ik][isp], numpy.complex_) * mu for isp in range(self.Nspinblocs)] # construct mupat mupat = [numpy.identity(self.n_orbitals[ik][isp], numpy.complex_) * mu for isp in range(self.Nspinblocs)] # construct mupat
Annkw = [numpy.zeros((self.n_orbitals[ik][isp], self.n_orbitals[ik][isp], N_om), dtype=numpy.complex_) for isp in range(self.Nspinblocs)] Annkw = [numpy.zeros((self.n_orbitals[ik][isp], self.n_orbitals[ik][isp], N_om), dtype=numpy.complex_) for isp in range(self.Nspinblocs)]
@ -448,43 +445,40 @@ class TransportEtOptic(SumkLDATools):
unchangesize = all([ self.n_orbitals[ik][isp] == mupat[isp].shape[0] for isp in range(self.Nspinblocs)]) unchangesize = all([ self.n_orbitals[ik][isp] == mupat[isp].shape[0] for isp in range(self.Nspinblocs)])
if (not unchangesize): if (not unchangesize):
# recontruct green functions. # recontruct green functions.
S = BlockGf(name_block_generator=[(bln[isp], GfReFreq(indices=range(self.n_orbitals[ik][isp]), S = BlockGf(name_block_generator=[(bln[isp], GfReFreq(indices=range(self.n_orbitals[ik][isp]), mesh=self.Sigma_imp[0].mesh)) for isp in range(self.Nspinblocs) ], make_copies=False)
mesh=self.Sigma_imp[0].mesh))
for isp in range(self.Nspinblocs) ],
make_copies=False)
# S = GF(name_block_generator=[ (s, GFBloc_ReFreq(Indices=BS, Mesh=self.Sigma_imp[0].mesh)) for s in ['up', 'down'] ], Copy=False) # S = GF(name_block_generator=[ (s, GFBloc_ReFreq(Indices=BS, Mesh=self.Sigma_imp[0].mesh)) for s in ['up', 'down'] ], Copy=False)
# mupat = numpy.identity(self.n_orbitals[ik], numpy.complex_) # change size of mupat # mupat = numpy.identity(self.n_orbitals[ik], numpy.complex_) # change size of mupat
mupat = [numpy.identity(self.n_orbitals[ik][isp], numpy.complex_) * mu for isp in range(self.Nspinblocs)] # construct mupat mupat = [numpy.identity(self.n_orbitals[ik][isp], numpy.complex_) * mu for isp in range(self.Nspinblocs)] # construct mupat
# mupat *= mu # mupat *= mu
#set a temporary array storing spectral functions with band index. Note, usually we should have spin index #set a temporary array storing spectral functions with band index. Note, usually we should have spin index
#Annkw=numpy.zeros((self.n_orbitals[ik],self.n_orbitals[ik],N_om),dtype=numpy.complex_) #Annkw=numpy.zeros((self.n_orbitals[ik],self.n_orbitals[ik],N_om),dtype=numpy.complex_)
Annkw = [numpy.zeros((self.n_orbitals[ik][isp], self.n_orbitals[ik][isp], N_om), dtype=numpy.complex_) for isp in range(self.Nspinblocs)] Annkw = [numpy.zeros((self.n_orbitals[ik][isp], self.n_orbitals[ik][isp], N_om), dtype=numpy.complex_) for isp in range(self.Nspinblocs)]
# get lattice green functions. # get lattice green functions.
# S <<= A_Omega_Plus_B(A=1, B=1j * broadening) # S <<= A_Omega_Plus_B(A=1, B=1j * broadening)
S <<= 1*Omega + 1j*broadening
S <<= 1*Omega + 1j*broadening
Ms = copy.deepcopy(mupat) Ms = copy.deepcopy(mupat)
for ibl in range(self.Nspinblocs): for ibl in range(self.Nspinblocs):
ind = ntoi[bln[ibl]] ind = ntoi[bln[ibl]]
n_orb = self.n_orbitals[ik][ibl] n_orb = self.n_orbitals[ik][ibl]
Ms[ibl] = self.hopping[ik,ind,0:n_orb,0:n_orb].real - mupat[ibl] Ms[ibl] = self.hopping[ik,ind,0:n_orb,0:n_orb].real - mupat[ibl]
S -= Ms S -= Ms
# print S[self.block_names[self.SO][0]].data
tmp = S.copy() # init temporary storage tmp = S.copy() # init temporary storage
## 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, stmp[icrsh][sig], gf)
S -= tmp S -= tmp
S.invert() S.invert()
for isp in range(self.Nspinblocs): for isp in range(self.Nspinblocs):
Annkw[isp].real = -copy.deepcopy(S[self.block_names[self.SO][isp]].data.swapaxes(0,1).swapaxes(1,2)).imag / numpy.pi Annkw[isp].real = -copy.deepcopy(S[self.block_names[self.SO][isp]].data.swapaxes(0,1).swapaxes(1,2)).imag / numpy.pi
for isp in range(self.Nspinblocs): for isp in range(self.Nspinblocs):
if(ik%100==0): if(ik%100==0):
print "ik,isp", ik, isp print "ik,isp", ik, isp
@ -508,9 +502,8 @@ class TransportEtOptic(SumkLDATools):
for iq in range(len(Qmesh_ex)): for iq in range(len(Qmesh_ex)):
#if(Qmesh_ex[iq]==0 or iw+Qmesh_ex[iq]>=N_om ): #if(Qmesh_ex[iq]==0 or iw+Qmesh_ex[iq]>=N_om ):
# here use fermi distribution to truncate self energy mesh. # here use fermi distribution to truncate self energy mesh.
if(Qmesh_ex[iq] == 0 or iw + Qmesh_ex[iq] >= N_om if(Qmesh_ex[iq] == 0 or iw + Qmesh_ex[iq] >= N_om or M[iw].real + Qmesh[iq] < -10.0 / Beta or M[iw].real >10.0 / Beta):
or M[iw].real + Qmesh[iq] < -10.0 / Beta or M[iw].real >10.0 / Beta): continue
continue
if (M[iw].real > omminplot) and (M[iw].real < ommaxplot): if (M[iw].real > omminplot) and (M[iw].real < ommaxplot):
# here use bandwin to construct match matrix for A and velocity. # here use bandwin to construct match matrix for A and velocity.
@ -558,7 +551,7 @@ class TransportEtOptic(SumkLDATools):
for i in xrange(L1): for i in xrange(L1):
for iq in xrange(L2): for iq in xrange(L2):
for iw in xrange(L3): for iw in xrange(L3):
pwout.write(str(self.Pw_optic[i, iq, iw]) + " ") pwout.write(str(self.Pw_optic[i, iq, iw]) + " ")
pwout.write("\n") pwout.write("\n")
# sum over omega to get optic conductivity for ik in xrange(self # sum over omega to get optic conductivity for ik in xrange(self
@ -701,14 +694,15 @@ class TransportEtOptic(SumkLDATools):
# get lattice green functions. # get lattice green functions.
# S <<= A_Omega_Plus_B(A=1, B=1j * broadening) # S <<= A_Omega_Plus_B(A=1, B=1j * broadening)
S <<= 1*Omega + 1j*broadening S <<= 1*Omega + 1j*broadening
Ms = copy.deepcopy(mupat) Ms = copy.deepcopy(mupat)
for ibl in range(self.Nspinblocs): for ibl in range(self.Nspinblocs):
ind = ntoi[bln[ibl]] ind = ntoi[bln[ibl]]
n_orb = self.n_orbitals[ik][ibl] n_orb = self.n_orbitals[ik][ibl]
Ms[ibl] = self.hopping[ik,ind,0:n_orb,0:n_orb].real - mupat[ibl] Ms[ibl] = self.hopping[ik,ind,0:n_orb,0:n_orb].real - mupat[ibl]
S -= Ms S -= Ms
# tmp = S.copy() # init temporary storage # tmp = S.copy() # init temporary storage
# ## substract self energy # ## substract self energy
@ -717,10 +711,10 @@ class TransportEtOptic(SumkLDATools):
# S -= tmp # S -= tmp
S.invert() S.invert()
for isp in range(self.Nspinblocs): for isp in range(self.Nspinblocs):
Annkw[isp].real = -copy.deepcopy(S[self.block_names[self.SO][isp]].data.swapaxes(0,1).swapaxes(1,2)).imag / numpy.pi Annkw[isp].real = -copy.deepcopy(S[self.block_names[self.SO][isp]].data.swapaxes(0,1).swapaxes(1,2)).imag / numpy.pi
for isp in range(self.Nspinblocs): for isp in range(self.Nspinblocs):
if(ik%100==0): if(ik%100==0):
print "ik,isp", ik, isp print "ik,isp", ik, isp
@ -777,21 +771,26 @@ class TransportEtOptic(SumkLDATools):
with open("TD_Optic_LDA.dat", "w") as pwout: with open("TD_Optic_LDA.dat", "w") as pwout:
L1,L2,L3=self.Pw_optic.shape L1,L2,L3=self.Pw_optic.shape
pwout.write("%s %s %s\n"%(L1,L2,L3)) pwout.write("%s %s %s\n"%(L1,L2,L3))
Qmeshr=[i*deltaM for i in Qmesh_ex]
for iq in xrange(L2):
pwout.write(str(Qmeshr[iq])+" ")
pwout.write("\n")
for iw in xrange(L3):
pwout.write(str(M[iw].real)+" ")
pwout.write("\n")
for i in xrange(L1): for i in xrange(L1):
for iq in xrange(L2): for iq in xrange(L2):
for iw in xrange(L3): for iw in xrange(L3):
pwout.write(str(i)+" "+str(Qmesh_ex[iq] * deltaM) + " " + str(M[iw]) + " ") pwout.write(str(self.Pw_optic[i, iq, iw]) + " ")
pwout.write(str(self.Pw_optic[i, iq, iw]) + " ") pwout.write("\n")
pwout.write("\n")
pwout.write("\n")
# sum over omega to get optic conductivity # sum over omega to get optic conductivity
if myMPI.is_master_node(): if myMPI.is_master_node():
OpticConductivity = numpy.zeros((mshape.sum(), len(Qmesh)), dtype=numpy.float_) OpticConductivity = numpy.zeros((mshape.sum(), len(Qmesh)), dtype=numpy.float_)
for im in range(mshape.sum()): for im in range(mshape.sum()):
for iq in range(len(Qmesh)): for iq in range(len(Qmesh)):
for iw in xrange(N_om): for iw in xrange(N_om):
omegaT = M[iw] * Beta omegaT = M[iw].real * Beta
omega_aug = Qmesh_ex[iq] * deltaM omega_aug = Qmesh_ex[iq] * deltaM
OpticConductivity[im, iq] += self.Pw_optic[im, iq, iw] * (fermidis(omegaT) - fermidis(omegaT + omega_aug * Beta)) / omega_aug OpticConductivity[im, iq] += self.Pw_optic[im, iq, iw] * (fermidis(omegaT) - fermidis(omegaT + omega_aug * Beta)) / omega_aug
OpticConductivity *= deltaM OpticConductivity *= deltaM