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

Current version of transport code (to be deleted later)

Browse Source
This commit is contained in:
Manuel Zingl 2014-10-28 10:22:44 +01:00
parent 8251d76553
commit 23096e0e76
2 changed files with 1373 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

942
python/SumK_LDA_Transport.py Normal file
View File

@ -0,0 +1,942 @@
################################################################################
#
# TRIQS: a Toolbox for Research in Interacting Quantum Systems
#
# Copyright (C) 2011 by M. Aichhorn, L. Pourovskii, V. Vildosola
#
# TRIQS is free software: you can redistribute it and/or modify it under the
# terms of the GNU General Public License as published by the Free Software
# Foundation, either version 3 of the License, or (at your option) any later
# version.
#
# TRIQS is distributed in the hope that it will be useful, but WITHOUT ANY
# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
# details.
#
# You should have received a copy of the GNU General Public License along with
# TRIQS. If not, see <http://www.gnu.org/licenses/>.
#
################################################################################
#=======================================================================================================================
# #################################################################
# Code for Transport/Optic calculations based on SumK_LDA... class
# by Xiaoyu Deng <xiaoyu.deng@gmail.com>
# #################################################################
#=======================================================================================================================
from types import *
import numpy
# NEW
import pytriqs.utility.dichotomy as Dichotomy
# OLD
#import pytriqs.Base.Utility.Dichotomy as Dichotomy
# NEW
#from pytriqs.gf.local.descriptors import A_Omega_Plus_B
from pytriqs.gf.local import *
# OLD
#from pytrigs.Base.GF_Local.GF import GF
#from pytriqs.Base.GF_local.GFBloc_ImFreq import GFBloc_ImFreq
#from pytriqs.Base.GF_Local.GFBloc_ReFreq import GFBloc_ReFreq
#from pytriqs.Base.GF_Local.GFBloc_ImTime import GFBloc_ImTime
#from pytriqs.Base.GF_Local import GF_Initializers
# NEW
#from pytriqs.operators.operators2 import *
# OLD
#from pytriqs.Solvers.Operators import *
# NEW
# NOT FOUND
# OLD
#from pytriqs.Base.Utility.myUtils import Sum
# NEW
import pytriqs.utility.mpi as myMPI
# OLD
#import pytriqs.Base.Utility.MPI as myMPI
from datetime import datetime
#from pytriqs.Wien2k.Symmetry import *
# NEW
from pytriqs.applications.dft.sumk_lda import *
from pytriqs.applications.dft.sumk_lda_tools import *
#k in xrange(self OLD
#from pytriqs.Wien2k.SumK_LDA import SumK_LDA
#from pytriqs.Wien2k.SumK_LDA_tools import SumK_LDA_tools
import string
import copy
import SumK_LDA_Transport_Wien2k_input as Wien
def Read_Fortran_File (filename):
""" Returns a generator that yields all numbers in the Fortran file as float, one by one"""
import os.path
if not(os.path.exists(filename)) : raise IOError, "File %s does not exists" % filename
for line in open(filename, 'r') :
for x in line.replace('D', 'E').split() :
yield string.atof(x)
def Read_Fortran_File2(filename):
""" Returns a generator that yields all numbers in the Fortran file as float, one by one"""
import os.path
if not(os.path.exists(filename)) : raise IOError, "File %s does not exists" % filename
for line in open(filename, 'r') :
for x in line.replace('D', 'E').split() :
try:
yield string.atof(x)
except GeneratorExit:
raise
except:
yield x.replace('E', 'D')
def fermidis(x):
return 1.0/(numpy.exp(x)+1)
# OLD
#class TransportEtOptic(SumK_LDA_tools):
# NEW
class TransportEtOptic(SumkLDATools):
"""transport and optic related functions
calculates distributions: Tr A(k,w) v(k) A(k, w) v(k) or Tr A(k,w) v(k) A(k, w+Omega) v(k)
.based on thi for ik in xrange(selfs distribution other properties could be obtained.
!!! worked for cases with spin orbital interaction.
!!! for non-SOI, need check! Be careful!.
"""
def initbyWien(self, wiencase):
""" read in necessary parameters from wien file.
symmetries:
volume:
velocities:
"""
# k_dep_projection is the general case!.
assert self.k_dep_projection == 1, "Not implemented!"
self.Nspinblocs = self.SP + 1 - self.SO
# suffix of wien2k output file
self.blocksuffix=[[""],["up","dn"]]
self.velocities=None
self.bandwin=None
self.Vol = None
self.symm = None
self.nsymm = None
if myMPI.is_master_node():
CASE = Wien.WienStruct(wiencase)
CASE.readSGsymm(wiencase)
if(self.Nspinblocs == 1): # paramagnetic , without SOI; or spin polarized with SOI
self.velocities = [Wien.Velocities(wiencase,bln) for bln in self.blocksuffix[0]]
else:
self.velocities = [Wien.Velocities(wiencase, bln) for bln in self.blocksuffix[1] ]
# read in band window for each k points. "CASE.oubwin".
bandwin=self.bandwinfromwiencase(wiencase=wiencase)
self.bandwin=bandwin
self.Vol = CASE.VolumePC
self.symm = CASE.symmcartesian
self.nsymm = CASE.nsymm
myMPI.barrier()
self.velocities=myMPI.bcast(self.velocities)
self.bandwin=myMPI.bcast(self.bandwin)
self.Vol=myMPI.bcast(self.Vol)
self.symm=myMPI.bcast(self.symm)
self.nsymm=myMPI.bcast(self.nsymm)
def Transportdistribution_Boltz(self, wiencase, mshape=None, broadening=0.01, energywindow=None, Nomega=1000, loadpw=False):
"""This is for transport calculation of Boltzmann theory
calculate \sum_k Tr vv delta(\omega-enk) which is transportdistribution in Boltzmann theory
Just use the LDA hamiltonian and velocity, DMFT self energy is not needed.
mshape defines the indices of directions. xx,yy,zz,xy,yz,zx.
mshape is 3x3 matrix, mshape[0,0]=1 --> calculate xx, mshape[1,1]=1 --> calculate yy, mshape[1,2]=1 --> calculate xy,
by default, xx is calculated.
"""
if mshape==None:
mshape=numpy.zeros(9).reshape(3,3)
mshape[0,0]=1
assert mshape.shape == (3, 3), "mshape should be 3x3"
velocities = self.velocities
mu = 0
assert self.k_dep_projection == 1, "k_dep_projection = 0 is NOT implemented!"
if energywindow is None:
omminplot = -1.0
ommaxplot = 1.0
else:
omminplot = energywindow[0]
ommaxplot = energywindow[1]
deltaomega = (ommaxplot - omminplot) / Nomega
# transport distribution :output P(\omega)_xy should has the same dimension as defined in mshape.
self.Pw = numpy.zeros((mshape.sum(), Nomega), dtype=numpy.float_)
mlist = []
for ir in xrange(3):
for ic in xrange(3):
if(mshape[ir][ic] == 1):
mlist.append((ir, ic))
if loadpw:
with open("TD_Boltz_mp.dat", "r") as pwin:
for iw in xrange(Nomega):
fstr = pwin.readline().split()
aomega = iw * deltaomega + omminplot
assert abs(float(fstr[0]) - aomega) <= 1e-8, "mesh not match when load transportdistribution"
for ipw in range(mshape.sum()): self.Pw[ipw, iw] = float(fstr[ipw + 1])
print "Blotz Pw loaded"
return
# for ik=0
ik = 0
# block_names for green function and self energy
bln = self.block_names[self.SO]
ntoi = self.names_to_ind[self.SO]
ikarray = numpy.array(range(self.n_k))
for isp in range(self.Nspinblocs):
for ik in myMPI.slice_array(ikarray):
n_orb = self.n_orbitals[ik][isp]
mupat = numpy.ones(n_orb, numpy.float_) * mu
Ms = copy.deepcopy(mupat)
ind = ntoi[bln[isp]]
Ms = numpy.diag(self.hopping[ik,ind,0:n_orb,0:n_orb].real) - mupat
if(ik%100==0):
print "ik,isp", ik, isp
kvel = velocities[isp].vks[ik]
# in general, bandwindows for Annkw and for velocities are not the same.
# one should make sure the same window are using before go further. otherwise the matrix size are not match.
Pwtem = numpy.zeros((mshape.sum(), Nomega), dtype=numpy.float_)
#symmetry loop
# this symmetrization could be done first to speed up... To be done.
for Rmat in self.symm:
# get new velocity.
Rkvel = copy.deepcopy(kvel.vel)
for vnb1 in xrange(kvel.bandwin[1] - kvel.bandwin[0] + 1):
Rkvel[vnb1][vnb1][:] = numpy.dot(Rmat, Rkvel[vnb1][vnb1][:])
ipw = 0
for (ir, ic) in mlist:
for iw in xrange(Nomega):
omega = deltaomega * iw + omminplot
bmin = max(self.bandwin[isp][ik, 0], kvel.bandwin[0])
bmax = min(self.bandwin[isp][ik, 1], kvel.bandwin[1])
for ib in xrange(bmax - bmin + 1):
ibb = ib + bmin - self.bandwin[isp][ik, 0]
ibv = ib + bmin - kvel.bandwin[0]
enk = Ms[ibb] - omega
Ag = 1.0 / numpy.sqrt(2 * numpy.pi) / broadening * numpy.exp(-enk ** 2 / 2.0 / broadening ** 2)
vkr = Rkvel[ibv][ibv][ir]
vkc = Rkvel[ibv][ibv][ic]
#Pwtem[ipw,iw]+=(vkr*vkc*Ag).real
Pwtem[ipw, iw] += vkr * vkc * Ag
ipw += 1
self.Pw += Pwtem * self.bz_weights[ik] / self.nsymm
myMPI.barrier()
self.Pw = myMPI.all_reduce(myMPI.world, self.Pw, lambda x, y : x + y)
# for non-magnetic case, the total weight is doubled because of spin degeneracy.
self.Pw *= (2 - self.SP)
# scattering is needed here.
self.Pw *= 1.0/numpy.pi/2.0/broadening
if myMPI.is_master_node():
with open("TD_Boltz.dat", "w") as pwout:
for iw in xrange(Nomega):
omega = deltaomega * iw + omminplot
pwout.write(str(omega) + " ")
for i in range(self.Pw.shape[0]):
pwout.write(str(self.Pw[i, iw]) + " ")
pwout.write("\n")
def Transportdistribution(self, wiencase, mshape=None, broadening=0.00, energywindow=None, loadpw=False):
"""calculate Tr A(k,w) v(k) A(k, w) v(k).
mshape defines the indices of directions. xx,yy,zz,xy,yz,zx.
mshape is 3x3 matrix, mshape[0,0]=1 --> calculate xx, mshape[1,1]=1 --> calculate yy, mshape[1,2]=1 --> calculate xy,
by default, xx is calculated.
"""
if mshape==None:
mshape=numpy.zeros(9).reshape(3,3)
mshape[0,0]=1
assert mshape.shape == (3, 3), "mshape should be 3x3"
assert hasattr(self, "Sigma_imp"), "Set Sigma First!"
velocities = self.velocities
mu = self.chemical_potential
if myMPI.is_master_node():
print "Chemical_Potential", mu
# use k-dependent-projections.
assert self.k_dep_projection == 1, "Not implemented!"
# form self energy from impurity self energy and double counting term.
stmp = self.add_dc()
#set mesh and energy range for spectrals functions
# one could construct self energy within only a small energy range when calculating transports.
M = [x for x in self.Sigma_imp[0].mesh]
N_om = len(M)
if energywindow is None:
omminplot = M[0].real - 0.001
ommaxplot = M[N_om - 1].real + 0.001
else:
omminplot = energywindow[0]
ommaxplot = energywindow[1]
# set mesh for Pw, only mesh in focused energyrange is needed. Mpw is just the index of mesh need in M
Mpw = [i for i in xrange(len(M)) if (M[i].real > omminplot and M[i].real < ommaxplot)]
# output P(\omega)_xy should has the same dimension as defined in mshape.
self.Pw = numpy.zeros((mshape.sum(), N_om), dtype=numpy.float)
mlist = []
for ir in xrange(3):
for ic in xrange(3):
if(mshape[ir][ic] == 1):
mlist.append((ir, ic))
ik = 0
# will construct G in the end; don't be mislead by
# nomenclature; S is sometimes sigma, sometimes sigma^-1, etc.
bln = self.block_names[self.SO]
ntoi = self.names_to_ind[self.SO]
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)
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)]
ikarray = numpy.array(range(self.n_k))
for ik in myMPI.slice_array(ikarray):
#for ik in xrange(self.n_k):
unchangesize = all([ self.n_orbitals[ik][isp] == mupat[isp].shape[0] for isp in range(self.Nspinblocs)])
if (not unchangesize):
# recontruct green functions.
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)
# change size of 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],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)]
# get lattice green functions.
# S <<= A_Omega_Plus_B(A=1, B=1j * broadening)
S <<= 1*Omega+1j*broadening
Ms = copy.deepcopy(mupat)
for ibl in range(self.Nspinblocs):
n_orb = self.n_orbitals[ik][ibl]
ind = ntoi[bln[ibl]]
Ms[ibl] = self.hopping[ik,ind,0:n_orb,0:n_orb].real - mupat[ibl]
S -= Ms
tmp = S.copy()
for icrsh in xrange(self.n_corr_shells):
for sig, gf in tmp: tmp[sig] <<= self.upfold(ik, icrsh, sig, stmp[icrsh][sig], gf)
S -= tmp
S.invert()
#hence we have A(k,\omega)_nn' for a special k points.
for isp in range(self.Nspinblocs):
# Annkw[isp].real = -copy.deepcopy(S[self.block_names[self.SO][isp]]._data.array).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
# A=-1/pi*Im G
# for different spin velocties might be different
for isp in range(self.Nspinblocs):
if(ik%100==0):
print "ik,isp", ik, isp
kvel = velocities[isp].vks[ik]
# in general, bandwindows for Annkw and for velocities are not the same.
# one should make sure the same window are using before go further. otherwise the matrix size are not match.
Pwtem = numpy.zeros((mshape.sum(), N_om), dtype=numpy.float_)
#symmetry loop
# how to symmetrize this part???
for Rmat in self.symm:
# get new velocity.
Rkvel = copy.deepcopy(kvel.vel)
for vnb1 in xrange(kvel.bandwin[1] - kvel.bandwin[0] + 1):
for vnb2 in xrange(kvel.bandwin[1] - kvel.bandwin[0] + 1):
Rkvel[vnb1][vnb2][:] = numpy.dot(Rmat, Rkvel[vnb1][vnb2][:])
ipw = 0
bmin = max(self.bandwin[isp][ik, 0], kvel.bandwin[0])
bmax = min(self.bandwin[isp][ik, 1], kvel.bandwin[1])
Astart = bmin - self.bandwin[isp][ik, 0]
Aend = bmax - self.bandwin[isp][ik, 0] + 1
vstart = bmin - kvel.bandwin[0]
vend = bmax - kvel.bandwin[0] + 1
for (ir, ic) in mlist:
#for iw in xrange(N_om):
for iw in Mpw:
#if (M[iw]>omminplot) and (M[iw]<ommaxplot):
# here use bandwin to construct match matrix for A and velocity.
Annkwt = Annkw[isp][Astart:Aend, Astart:Aend, iw]
Rkveltr = Rkvel[vstart:vend, vstart:vend, ir]
Rkveltc = Rkvel[vstart:vend, vstart:vend, ic]
#print Annkwt.shape,Rkvel[...,ir].shape
Pwtem[ipw, iw] += numpy.dot(numpy.dot(numpy.dot(Rkveltr, Annkwt), Rkveltc), Annkwt).trace().real
ipw += 1
# k sum and spin sum.
self.Pw += Pwtem * self.bz_weights[ik] / self.nsymm
self.Pw = myMPI.all_reduce(myMPI.world, self.Pw, lambda x, y : x + y)
# for non-magnetic case, the total weight is doubled because of spin degeneracy.
self.Pw *= (2 - self.SP)
if myMPI.is_master_node():
with open("TD_DMFT.dat", "w") as pwout:
for iw in xrange(N_om):
if (M[iw].real > omminplot) and (M[iw].real < ommaxplot):
pwout.write(str(M[iw].real) + " ")
for i in range(self.Pw.shape[0]):
pwout.write(str(self.Pw[i, iw]) + " ")
pwout.write("\n")
def OpticDistribution(self, wiencase, mshape=None, broadening=0.01, energywindow=None, Qmesh=[0.5], Beta=50, loadpw=False):
"""calculate Tr A(k,w) v(k) A(k, w+q) v(k) and optics.
energywindow is the regime for omega integral
Qmesh contains the frequencies of the optic conductivitity. I repin the Qmesh to the self-energy mesh,
so the exact value might not exactly the same as given in the list.
mshape defines the indices of directions. xx,yy,zz,xy,yz,zx.
mshape is 3x3 matrix, mshape[0,0]=1 --> calculate xx, mshape[1,1]=1 --> calculate yy, mshape[1,2]=1 --> calculate xy,
by default, xx is calculated.
"""
assert mshape.shape == (3, 3), "mshape should be 3x3"
assert hasattr(self, "Sigma_imp"), "Set Sigma First!"
assert ((self.SP == 0) and (self.SO == 0)), "For SP and SO implementation of spaghettis has to be changed!"
velocities = self.velocities
# calculate A(k,w):
mu = self.chemical_potential
#we need this somehow for k_dep_projections. So we have to face the problem that the size of A(k,\omega) will
#change, and more, the band index for the A(k,\omega) matrix is not known yet.
# use k-dependent-projections.
assert self.k_dep_projection == 1, "Not implemented!"
# form self energy from impurity self energy and double counting term.
stmp = self.add_dc()
#set mesh and energyrange.
M = [x for x in self.Sigma_imp[0].mesh]
deltaM = numpy.abs(M[0] - M[1])
N_om = len(M)
if energywindow is None:
omminplot = M[0].real - 0.001
ommaxplot = M[N_om - 1].real + 0.001
else:
omminplot = energywindow[0]
ommaxplot = energywindow[1]
# define exact mesh for optic conductivity
Qmesh_ex = [int(x / deltaM) for x in Qmesh]
if myMPI.is_master_node():
print "Qmesh ", Qmesh
print "mesh interval in self energy ", deltaM
print "Qmesh / mesh interval ", Qmesh_ex
# output P(\omega)_xy should has the same dimension as defined in mshape.
self.Pw_optic = numpy.zeros((mshape.sum(), len(Qmesh), N_om), dtype=numpy.float_)
mlist = []
for ir in xrange(3):
for ic in xrange(3):
if(mshape[ir][ic] == 1):
mlist.append((ir, ic))
ik = 0
bln = self.block_names[self.SO]
ntoi = self.names_to_ind[self.SO]
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)
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)]
ikarray = numpy.array(range(self.n_k))
for ik in myMPI.slice_array(ikarray):
unchangesize = all([ self.n_orbitals[ik][isp] == mupat[isp].shape[0] for isp in range(self.Nspinblocs)])
if (not unchangesize):
# recontruct green functions.
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)
# 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][isp], numpy.complex_) * mu for isp in range(self.Nspinblocs)] # construct mupat
# mupat *= mu
#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][isp], self.n_orbitals[ik][isp], N_om), dtype=numpy.complex_) for isp in range(self.Nspinblocs)]
# get lattice green functions.
# S <<= A_Omega_Plus_B(A=1, B=1j * broadening)
S <<= 1*Omega + 1j*broadening
Ms = copy.deepcopy(mupat)
for ibl in range(self.Nspinblocs):
ind = ntoi[bln[ibl]]
n_orb = self.n_orbitals[ik][ibl]
Ms[ibl] = self.hopping[ik,ind,0:n_orb,0:n_orb].real - mupat[ibl]
S -= Ms
tmp = S.copy() # init temporary storage
## substract self energy
for icrsh in xrange(self.n_corr_shells):
for sig, gf in tmp: tmp[sig] <<= self.upfold(ik, icrsh, sig, stmp[icrsh][sig], gf)
S -= tmp
S.invert()
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
for isp in range(self.Nspinblocs):
if(ik%100==0):
print "ik,isp", ik, isp
kvel = velocities[isp].vks[ik]
Pwtem = numpy.zeros((mshape.sum(), len(Qmesh_ex), N_om), dtype=numpy.float_)
#symmetry loop
for Rmat in self.symm:
# get new velocity.
Rkvel = copy.deepcopy(kvel.vel)
for vnb1 in xrange(kvel.bandwin[1] - kvel.bandwin[0] + 1):
for vnb2 in xrange(kvel.bandwin[1] - kvel.bandwin[0] + 1):
Rkvel[vnb1][vnb2][:] = numpy.dot(Rmat, Rkvel[vnb1][vnb2][:])
ipw = 0
for (ir, ic) in mlist:
for iw in xrange(N_om):
if(M[iw].real > 5.0 / Beta):
continue
for iq in range(len(Qmesh_ex)):
#if(Qmesh_ex[iq]==0 or iw+Qmesh_ex[iq]>=N_om ):
# here use fermi distribution to truncate self energy mesh.
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):
continue
if (M[iw].real > omminplot) and (M[iw].real < ommaxplot):
# here use bandwin to construct match matrix for A and velocity.
bmin = max(self.bandwin[isp][ik, 0], kvel.bandwin[0])
bmax = min(self.bandwin[isp][ik, 1], kvel.bandwin[1])
Astart = bmin - self.bandwin[isp][ik, 0]
Aend = bmax - self.bandwin[isp][ik, 0] + 1
vstart = bmin - kvel.bandwin[0]
vend = bmax - kvel.bandwin[0] + 1
Annkwl = Annkw[isp][Astart:Aend, Astart:Aend, iw]
Annkwr = Annkw[isp][Astart:Aend, Astart:Aend, iw + Qmesh_ex[iq]]
Rkveltr = Rkvel[vstart:vend, vstart:vend, ir]
Rkveltc = Rkvel[vstart:vend, vstart:vend, ic]
#print Annkwt.shape,Rkvel[...,ir].shape
Pwtem[ipw, iq, iw] += numpy.dot(numpy.dot(numpy.dot(Rkveltr, Annkwl), Rkveltc), Annkwr).trace().real
ipw += 1
# k sum and spin sum.
self.Pw_optic += Pwtem * self.bz_weights[ik] / self.nsymm
self.Pw_optic = myMPI.all_reduce(myMPI.world, self.Pw_optic, lambda x, y : x + y)
self.Pw_optic *= (2 - self.SP)
# just back up TD_optic data
if myMPI.is_master_node():
with open("TD_Optic_DMFT.dat", "w") as pwout:
#shape
L1,L2,L3=self.Pw_optic.shape
pwout.write("%s %s %s\n"%(L1,L2,L3))
#dump Qmesh
Qmeshr=[i*deltaM for i in Qmesh_ex]
#dump self energy mesh
#dump Pw_optic
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 iq in xrange(L2):
for iw in xrange(L3):
pwout.write(str(self.Pw_optic[i, iq, iw]) + " ")
pwout.write("\n")
# sum over omega to get optic conductivity for ik in xrange(self
if myMPI.is_master_node():
OpticConductivity = numpy.zeros((mshape.sum(), len(Qmesh)), dtype=numpy.float_)
for im in range(mshape.sum()):
for iq in range(len(Qmesh)):
for iw in xrange(N_om):
omegaT = M[iw].real * Beta
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 *= deltaM
OpticConductivity *= 10700 / self.Vol
with open("Optic_con.dat", "wOptic_con") as opt:
for iq in range(len(Qmesh_ex)):
opt.write(str(Qmesh_ex[iq] * deltaM) + " ")
for im in range(mshape.sum()):
opt.write(str(OpticConductivity[im, iq]) + " ")
opt.write("\n")
def loadOpticTD(self,OpticTDFile="TD_Optic_DMFT.dat",Beta=40):
""" load optic conductivity distribution and calculate Optical Conductivty
"""
if myMPI.is_master_node():
with open(OpticTDFile,"r") as pw:
L1,L2,L3=(int(i) for i in pw.readline().split())
#QMeshr=numpy.zeros(L2,dtype=numpy.float)
#M=numpy.zeros(L3,dtype=numpy.float)
#Pw_optic=numpy.zeros((L1,L2,L3), dtype=numpy.float)
QMeshr=numpy.array([float(i) for i in pw.readline().split()])
M=numpy.array([float(i) for i in pw.readline().split()])
Pw_optic=numpy.array([float(i) for i in pw.readline().split()]).reshape(L1,L2,L3)
OpticConductivity = numpy.zeros((L1, L2), dtype=numpy.float)
deltaM=M[1]-M[0]
for im in xrange(L1):
for iq in xrange(L2):
for iw in xrange(L3):
omegaT = M[iw] * Beta
omega_aug = QMeshr[iq]
OpticConductivity[im, iq] += Pw_optic[im, iq, iw] * (fermidis(omegaT) - fermidis(omegaT + omega_aug * Beta)) / omega_aug
OpticConductivity *= deltaM
## transform to standard unit as in resistivity
OpticConductivity *= 10700 / self.Vol
##
with open("Optic_con.dat", "w") as opt:
for iq in xrange(L2):
opt.write(str(QMeshr[iq]) + " ")
for im in xrange(L1):
opt.write(str(OpticConductivity[im, iq]) + " ")
opt.write("\n")
def OpticDistribution_LDA(self, wiencase, mshape=None, broadening=0.01, energywindow=None, Qmesh=[0.5], Beta=50, loadpw=False):
"""calculate Tr A(k,w) v(k) A(k, w+q) v(k) and optics. A constant self-energy is used to mimick noninteracting case.
It is not the best way to calculate optic for LDA. Just to compare.
energywindow is the regime for omega integral
Qmesh contains the frequencies of the optic conductivitity. I repin the Qmesh to the self-energy mesh,
so the exact value might not exactly the same as given in the list.
mshape defines the indices of directions. xx,yy,zz,xy,yz,zx.
mshape is 3x3 matrix, mshape[0,0]=1 --> calculate xx, mshape[1,1]=1 --> calculate yy, mshape[1,2]=1 --> calculate xy,
by default, xx is calculated.
"""
assert mshape.shape == (3, 3), "mshape should be 3x3"
assert hasattr(self, "Sigma_imp"), "Set Sigma First!"
assert ((self.SP == 0) and (self.SO == 0)), "For SP and SO implementation of spaghettis has to be changed!"
velocities = self.velocities
# calculate A(k,w):
mu = self.chemical_potential
#we need this somehow for k_dep_projections. So we have to face the problem that the size of A(k,\omega) will
#change, and more, the band index for the A(k,\omega) matrix is not known yet.
# use k-dependent-projections.
assert self.k_dep_projection == 1, "Not implemented!"
# form self energy from impurity self energy and double counting term.
stmp = self.add_dc()
#set mesh and energyrange.
M = [x for x in self.Sigma_imp[0].mesh]
deltaM = numpy.abs(M[0] - M[1])
N_om = len(M)
if energywindow is None:
omminplot = M[0] - 0.001
ommaxplot = M[N_om - 1] + 0.001
else:
omminplot = energywindow[0]
ommaxplot = energywindow[1]
# define exact mesh for optic conductivity
Qmesh_ex = [int(x / deltaM) for x in Qmesh]
if myMPI.is_master_node():
print "Qmesh ", Qmesh
print "mesh interval in self energy ", deltaM
print "Qmesh / mesh interval ", Qmesh_ex
# output P(\omega)_xy should has the same dimension as defined in mshape.
self.Pw_optic = numpy.zeros((mshape.sum(), len(Qmesh), N_om), dtype=numpy.float_)
mlist = []
for ir in xrange(3):
for ic in xrange(3):
if(mshape[ir][ic] == 1):
mlist.append((ir, ic))
ik = 0
bln = self.block_names[self.SO]
ntoi = self.names_to_ind[self.SO]
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)
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)]
ikarray = numpy.array(range(self.n_k))
for ik in myMPI.slice_array(ikarray):
unchangesize = all([ self.n_orbitals[ik][isp] == mupat[isp].shape[0] for isp in range(self.Nspinblocs)])
if (not unchangesize):
# recontruct green functions.
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)
# 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][isp], numpy.complex_) * mu for isp in range(self.Nspinblocs)] # construct mupat
# mupat *= mu
#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][isp], self.n_orbitals[ik][isp], N_om), dtype=numpy.complex_) for isp in range(self.Nspinblocs)]
# get lattice green functions.
# S <<= A_Omega_Plus_B(A=1, B=1j * broadening)
S <<= 1*Omega + 1j*broadening
Ms = copy.deepcopy(mupat)
for ibl in range(self.Nspinblocs):
ind = ntoi[bln[ibl]]
n_orb = self.n_orbitals[ik][ibl]
Ms[ibl] = self.hopping[ik,ind,0:n_orb,0:n_orb].real - mupat[ibl]
S -= Ms
# tmp = S.copy() # init temporary storage
# ## substract self energy
# for icrsh in xrange(self.n_corr_shells):
# for sig, gf in tmp: tmp[sig] <<= self.upfold(ik, icrsh, sig, stmp[icrsh][sig], gf)
# S -= tmp
S.invert()
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
for isp in range(self.Nspinblocs):
if(ik%100==0):
print "ik,isp", ik, isp
kvel = velocities[isp].vks[ik]
Pwtem = numpy.zeros((mshape.sum(), len(Qmesh_ex), N_om), dtype=numpy.float_)
#symmetry loop
for Rmat in self.symm:
# get new velocity.
Rkvel = copy.deepcopy(kvel.vel)
for vnb1 in xrange(kvel.bandwin[1] - kvel.bandwin[0] + 1):
for vnb2 in xrange(kvel.bandwin[1] - kvel.bandwin[0] + 1):
Rkvel[vnb1][vnb2][:] = numpy.dot(Rmat, Rkvel[vnb1][vnb2][:])
ipw = 0
for (ir, ic) in mlist:
for iw in xrange(N_om):
if(M[iw].real > 5.0 / Beta):
continue
for iq in range(len(Qmesh_ex)):
#if(Qmesh_ex[iq]==0 or iw+Qmesh_ex[iq]>=N_om ):
# here use fermi distribution to truncate self energy mesh.
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):
continue
if (M[iw].real > omminplot) and (M[iw].real < ommaxplot):
# here use bandwin to construct match matrix for A and velocity.
bmin = max(self.bandwin[isp][ik, 0], kvel.bandwin[0])
bmax = min(self.bandwin[isp][ik, 1], kvel.bandwin[1])
Astart = bmin - self.bandwin[isp][ik, 0]
Aend = bmax - self.bandwin[isp][ik, 0] + 1
vstart = bmin - kvel.bandwin[0]
vend = bmax - kvel.bandwin[0] + 1
Annkwl = Annkw[isp][Astart:Aend, Astart:Aend, iw]
Annkwr = Annkw[isp][Astart:Aend, Astart:Aend, iw + Qmesh_ex[iq]]
Rkveltr = Rkvel[vstart:vend, vstart:vend, ir]
Rkveltc = Rkvel[vstart:vend, vstart:vend, ic]
#print Annkwt.shape,Rkvel[...,ir].shape
Pwtem[ipw, iq, iw] += numpy.dot(numpy.dot(numpy.dot(Rkveltr, Annkwl), Rkveltc), Annkwr).trace().real
ipw += 1
# k sum and spin sum.
self.Pw_optic += Pwtem * self.bz_weights[ik] / self.nsymm
self.Pw_optic = myMPI.all_reduce(myMPI.world, self.Pw_optic, lambda x, y : x + y)
self.Pw_optic *= (2 - self.SP)
# just back up TD_optic data # just back up TD_optic data
if myMPI.is_master_node():
with open("TD_Optic_LDA.dat", "w") as pwout:
L1,L2,L3=self.Pw_optic.shape
pwout.write("%s %s %s\n"%(L1,L2,L3))
for i in xrange(L1):
for iq in xrange(L2):
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("\n")
pwout.write("\n")
# sum over omega to get optic conductivity
if myMPI.is_master_node():
OpticConductivity = numpy.zeros((mshape.sum(), len(Qmesh)), dtype=numpy.float_)
for im in range(mshape.sum()):
for iq in range(len(Qmesh)):
for iw in xrange(N_om):
omegaT = M[iw] * Beta
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 *= deltaM
OpticConductivity *= 10700 / self.Vol
with open("Optic_con_LDA.dat", "w") as opt:
for iq in range(len(Qmesh_ex)):
opt.write(str(Qmesh_ex[iq] * deltaM) + " ")
for im in range(mshape.sum()):
opt.write(str(OpticConductivity[im, iq]) + " ")
opt.write("\n")
# load transportdistribution Pw from file.
def loadTD(self, filename, fermishift=0.0):
""" load transport distribution from file. Assume energy mesh is uniform
the first column is energy mesh. The others are TD values.
fermishift is used to shift the TD mesh to mimick a rigid band shift.
"""
myMPI.barrier()
self.TD = numpy.loadtxt(filename)
self.TD[:, 0] -= fermishift
# seebeck is just intetral of omega \int Pw f(\omega)f(-\omega)(\beta w).
def Seebeck(self, Beta, index=0):
""" get -A1/A0, that is Seebeck in unit k_B/e. index is used to select the column in self.tdintegral.
Note: for nodiagonal element in Sxy this might not be right. so take care.
"""
if myMPI.is_master_node():
print "A0, A1 %.5e %.5e " % (self.tdintegral(Beta, 0)[index], self.tdintegral(Beta, 1)[index])
seb = -self.tdintegral(Beta, 1)[index] / self.tdintegral(Beta, 0)[index]
print "Seebeck%d %.4f k_B/e %.4f x 10^(-6)V/K" % (index, seb, seb * 86.17)
return seb
def Conductivity(self, Beta, index=0):
""" #return 1/T*A0, that is Conductivity in unit 1/V
return Conductivity
"""
if myMPI.is_master_node():
Cond = Beta * self.tdintegral(Beta, 0)[index]
#print "Beta*A0 ", Cond
print "V in bohr^3 ", self.Vol
Cond *= 10700.0 / self.Vol
print "Conductivity%d %.4f x 10^4 Ohm^-1 cm^-1" % (index, Cond)
print "Resistivity%d %.4f x 10^-4 Ohm cm" % (index, 1.0 / Cond)
return Cond
def tdintegral(self, Beta, pn=0):
"""calculate { \pi *\int Pw f(omega)f(-omega)(\beta\omega)^(pn)d\omega }
"""
M = self.TD[:, 0]
domega = abs(M[1] - M[0])
pwint = numpy.zeros(self.TD.shape[1] - 1)
for ipw in range(self.TD.shape[1] - 1):
for iw in xrange(len(M)):
x = M[iw] * Beta
pwint[ipw] += fermidis(x) * fermidis(-x) * self.TD[iw, ipw + 1] * numpy.float(x) ** pn * domega
return pwint
def tdintcore(self, Beta):
"""calculate { Pw f(omega)f(-omega) for check data)
"""
M = self.TD[:, 0]
domega = abs(M[1] - M[0])
pwint = numpy.zeros((self.TD.shape[1] - 1, M.size))
for ipw in range(self.TD.shape[1] - 1):
for iw in xrange(M.size):
x = M[iw] * Beta
pwint[ipw, iw] = self.TD[iw, ipw + 1] * fermidis(x) * fermidis(-x)
#if(pwint[ipw,iw]>=0.3):
# self.Pw[ipw,iw]=0.0
if myMPI.is_master_node():
with open("tdintcore.dat", "w") as pwout:
for iw in xrange(M.size):
pwout.write(str(M[iw]) + " ")
for i in range(pwint.shape[0]):
pwout.write(str(pwint[i, iw]) + " ")
pwout.write("\n")
return pwint
def bandwinfromwiencase(self, wiencase):
""" read in the band window from wiencase.outbwin file.
"""
bandwin = [numpy.zeros(self.n_k * 2, dtype=int).reshape(self.n_k, 2) for isp in range(self.SP + 1 - self.SO)]
for isp in range(self.SP + 1 - self.SO):
if(self.SP == 0 or self.SO == 1):
winfile = Read_Fortran_File2(wiencase + ".oubwin")
elif self.SP == 1 and isp == 0:
winfile = Read_Fortran_File2(wiencase + ".oubwinup")
elif self.SP == 1 and isp == 1:
winfile = Read_Fortran_File2(wiencase + ".oubwindn")
else:
assert 0, "Reading bandwin error! Check self.SP and self.SO!"
Nk = int(winfile.next())
assert Nk == self.n_k, "Number of K points is unconsistent in case.oubwin"
SO = int(winfile.next())
assert SO == self.SO, "SO is unconsistent in case.oubwin"
for i in xrange(self.n_k):
winfile.next()
bandwin[isp][i, 0] = winfile.next()
bandwin[isp][i, 1] = winfile.next()
winfile.next()
return bandwin

431
python/SumK_LDA_Transport_Wien2k_input.py Normal file
View File

@ -0,0 +1,431 @@
################################################################################
#
# TRIQS: a Toolbox for Research in Interacting Quantum Systems
#
# Copyright (C) 2011 by M. Aichhorn, L. Pourovskii, V. Vildosola
#
# TRIQS is free software: you can redistribute it and/or modify it under the
# terms of the GNU General Public License as published by the Free Software
# Foundation, either version 3 of the License, or (at your option) any later
# version.
#
# TRIQS is distributed in the hope that it will be useful, but WITHOUT ANY
# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
# details.
#
# You should have received a copy of the GNU General Public License along with
# TRIQS. If not, see <http://www.gnu.org/licenses/>.
#
################################################################################
#=======================================================================================================================
# #################################################################
# Code for Transport/Optic calculations based on SumK_LDA... class
# by Xiaoyu Deng <xiaoyu.deng@gmail.com>
# The code read in files needed for transport/Optic calculations from Wien outputs.
# including: symmetry, velocity, lattice constants.
# The HDF convention is not adopted here since momentum file from Wien output is usually quite large
# and it is not necessary to keep in in the HDF file.
# #################################################################
#=======================================================================================================================
import numpy, sys, os.path
def Read_Fortran_File (filename):
""" Returns a generator that yields all numbers in the Fortran file as float, one by one"""
import os.path
if not(os.path.exists(filename)) : raise IOError, "File %s does not exists" % filename
for line in open(filename, 'r') :
for x in line.replace('D', 'E').split() :
yield string.atof(x)
class Velocity_k:
"""momentum matrix for a single k points.
"""
def __init__(self):
self.kp = [0.0, 0.0, 0.0]
self.bandwin = [999, 0]
# here a matrix for 3 vels. use list since the size of vel is to be determined.
self.vel = []
class Velocities:
"""Class containig the velocities.
Provides container for velocities
as a function of k and method to read them from case.pmat Wien2k file
use as ClassInstance.vks[ik].vel[iband][jband][ix]
"""
def __init__(self, wiencase, spinbl=""):
if not(os.path.exists(wiencase + ".pmat" + spinbl)) : raise IOError, "File %s does not exists" % wiencase + ".pmat"
# expected format:
# k nu nu2 (k denotes k-point, nu1 starting band index, nu2 ending band index
# kpos ( read in if needed.
# (real, im) (of velocity for all nu <=nu_prime within nu1 and nu2)
# ... ...
# k nu nu2
# ...
f = open(wiencase + ".pmat" + spinbl)
self.vks = []
while 1:
try:
s = f.readline()
if (s == ""):
break
except:
break
try:
vec = Velocity_k()
[k, nu1, nu2] = [int (x) for x in s.strip().split()]
vec.bandwin[0] = nu1
vec.bandwin[1] = nu2
vec.kp = f.readline().strip().split()
dim = vec.bandwin[1] - vec.bandwin[0] + 1
shape = (dim, dim, 3)
vxyz = numpy.zeros(shape, dtype=complex)
for nu in xrange(dim):
for nu_prime in xrange(nu, dim):
for i in xrange(3):
s = f.readline().strip("\n ()").split(',')
vxyz[nu][nu_prime][i] = float(s[0]) + float(s[1]) * 1j
if(nu_prime != nu):
vxyz[nu_prime][nu][i] = vxyz[nu][nu_prime][i].conjugate()
vec.vel = vxyz
self.vks.append(vec)
except IOError:
print("Reading case.pmat error. Wrong format?\n ")
raise
f.close()
def getvel(self, k):
# should return array at a given k. use as [iband][jband][i]
return self.vks[k].vel
def plot(self):
f = open("velband.dat", "w")
bandid = numpy.array([vk.bandwin for vk in self.vks]).flatten()
minb = bandid.min()
maxb = bandid.max()
for ib in range(minb, maxb + 1):
ik = 0
for vk in self.vks:
if(vk.bandwin[0] <= ib and vk.bandwin[1] >= ib):
f.write(str(ik) + " ")
for i in range(3):
f.write(str(vk.vel[ib - vk.bandwin[0]][ib - vk.bandwin[0]][i].real) + " ")
f.write("\n")
ik = ik + 1
f.write("&\n")
f.close()
class SGsymmetry():
""" read in symmetry of space group from wiencase.outputs other than wiencase.struct since
in this file there are also symmetry operations in xyz coordinates..
"""
def __init__(self, wiencase):
structfile = wiencase + ".outputs"
self.nsymm = 1
self.symm = []
self.tau = []
self.bravaismatrix = numpy.zeros((3, 3), dtype=numpy.float_)
self.symmcartesian = []
self.taucartesian = []
with open(structfile, "r") as f:
f.readline()
f.readline()# bravais matrix
for i in range(3):
line = f.readline().strip().split()
self.bravaismatrix[i, :] = numpy.array([numpy.float(item) for item in line])[:]
print "bravais matrix", self.bravaismatrix
while 1:
s = f.readline().strip(" ").split()
try:
if(s[0] == "PGBSYM:"):
self.nsymm = int(s[-1])
break
except:
continue
f.readline()
f.readline()
for i in range(self.nsymm):
f.readline()
## read symmcartesian
symmt = numpy.zeros((3, 3), dtype=float)
taut = numpy.zeros((3), dtype=float)
for ir in range(3):
s = f.readline().strip().split()
for ic in range(3):
symmt[ir, ic] = float(s[ic])
s = f.readline().strip().split()
for ir in range(3):
taut[ir] = float(s[ir])
self.symmcartesian.append(symmt)
self.taucartesian.append(taut)
##read symm
symmt = numpy.zeros((3, 3), dtype=float)
taut = numpy.zeros((3), dtype=float)
for ir in range(3):
s = f.readline().strip().split()
for ic in range(3):
symmt[ir, ic] = float(s[ic])
taut[ir] = float(s[3])
self.symm.append(symmt)
self.tau.append(taut)
f.readline()
# end
f.close()
print "Read wiencase.outputs done!"
def checksymmxyz(self):
''' This is to check symm in cartesian coordinates and in primitive cell lattice
For details of symm, one should check wien/SRC_symmetry/, latsym.f, pglsym.f,pgbsym.f
In general, if a lattice is orthorhombic, then symmcartesian is the same as symmprimitive
(at most different by a transpose (or inversion)). If a lattice is not orthorhombic, these two symms could
be related by bravais matrix.
One special case is CXZ lattice. In Wien CXZ lattice contains two cases: orthorhombic, and
monoclinic (with only gamma not equal to 90). For CXZ lattice monoclinic, symmcartesian is also the same as
symmprimitive (or with transpose (or inversion)).
'''
for i in range(self.nsymm):
mat = self.symmcartesian[i].transpose() # according to wien2k, why?
bm = self.bravaismatrix
bminv = numpy.linalg.inv(bm)#.transpose()
res = numpy.dot(bminv, numpy.dot(mat, bm)) - self.symm[i]
print i, res
def size(self):
return self.nsymm
def cellvolume(latticetype, latticeconstants, latticeangle):
""" calculate cell volume: volumecc conventional cell, volumepc, primitive cell.
"""
for i in range(3):
latticeangle[i] *= 1.0 / 180 * numpy.pi
a = latticeconstants[0]
b = latticeconstants[1]
c = latticeconstants[2]
c_al = numpy.cos(latticeangle[0])
c_be = numpy.cos(latticeangle[1])
c_ga = numpy.cos(latticeangle[2])
volumecc = a * b * c * numpy.sqrt(1 + 2 * c_al * c_be * c_ga - c_al ** 2 - c_be * 82 - c_ga ** 2)
det = {"P":1,
"F":4,
"B":2,
"R":3,
"H":1,
"CXY":2,
"CYZ":2,
"CXZ":2
}
volumepc = volumecc / det[latticetype]
return volumecc, volumepc
class WienStruct():
""" parsing Wien Struct file
"""
def __init__(self, wiencase):
structfile = wiencase + ".struct"
with open(structfile, "r") as infile:
print "read in Wien case file %s" % structfile
infile.readline()#title
tem = infile.readline() #lattice
self.latticetype = tem[0:10].split()[0]
self.ineqvsite = int(tem[27:30])
try:
self.sgrnumber = int(tem[30:33])
except:
self.sgrnumber = None
try:
self.sgrlabel = tem[34:38]
except:
self.sgrlabel = None
print self.latticetype, self.ineqvsite, self.sgrnumber, self.sgrlabel
infile.readline()
tem = infile.readline() # lattice constants
self.latticeconstants = [float(tem[0:10]), float(tem[10:20]), float(tem[20:30])]
self.latticeangle = [float(tem[30:40]), float(tem[40:50]), float(tem[50:60])]
print "Cell"
print self.latticeconstants[:]
print self.latticeangle[:]
self.positions = []
self.atomsymbols = []
self.multi = []
self.atomnumbers = []
self.locrotmatrix = []
for isite in range(self.ineqvsite):
tem = infile.readline()
positions = []
positions.append([float(tem[12:22]), float(tem[25:35]), float(tem[38:48])])
tem = infile.readline()
multi = int(tem[15:17])
self.multi.append(multi)
for im in range(multi - 1):
tem = infile.readline()
positions.append([float(tem[12:22]), float(tem[25:35]), float(tem[38:48])])
#print positions
self.positions.append(positions)
#print self.positions
tem = infile.readline().strip(" ").split()
self.atomsymbols.append(tem[0])
self.atomnumbers.append(float(tem[-1]))
mat = numpy.zeros((3, 3), dtype=numpy.float_)
tem = infile.readline().strip(" ").split()
#print tem[-3:],mat[0,:]
mat[0, :] = tem[-3:]
tem = infile.readline().strip(" ").split()
mat[1, :] = tem[-3:]
tem = infile.readline().strip(" ").split()
mat[2, :] = tem[-3:]
self.locrotmatrix.append(mat)
for ia in range(len(self.atomsymbols)):
print "atom symbol : %s atom number: %d atom multi: %d" % (self.atomsymbols[ia], self.atomnumbers[ia], self.multi[ia])
print "positions:"
for im in self.positions[ia]:
print im[:]
# symmetry with lattice vector
tem = infile.readline().strip(" ").split()
self.symm = []
self.tau = []
self.Nsymm = int(tem[0])
for isymm in range(self.Nsymm):
symmt = numpy.zeros((3, 3), dtype=float)
taut = numpy.zeros((3), dtype=float)
for ir in range(3):
s = infile.readline()
for ic in range(3):
symmt[ir][ic] = float(s[ic * 2:ic * 2 + 2])
taut[ir] = float(s[7:17])
self.symm.append(symmt)
self.tau.append(taut)
infile.readline()
#############
print "Read in %s.struct done!" % wiencase
## convential Cell Volume and primitive Cell. In bohr^3 unit
self.VolumeCC, self.VolumePC = cellvolume(self.latticetype, self.latticeconstants, self.latticeangle)
def readSGsymm(self, wiencase):
""" read in symmetry of space group from wiencase.outputs other than wiencase.struct since
in this file there are also symmetry operations in xyz coordinates..
"""
structfile = wiencase + ".outputs"
self.nsymm = 1
self.symm = []
self.tau = []
# note bravaismatrix is not accurate enough in wiencase.outputs file. Just use it for test.
self.bravaismatrix = numpy.zeros((3, 3), dtype=numpy.float_)
self.symmcartesian = []
self.taucartesian = []
with open(structfile, "r") as f:
f.readline()
f.readline()# bravais matrix
for i in range(3):
line = f.readline().strip().split()
self.bravaismatrix[i, :] = numpy.array([numpy.float(item) for item in line])[:]
print "bravais matrix", self.bravaismatrix
while 1:
try:
s = f.readline().strip(" ").split()
if(s[0] == "PGBSYM:"):
self.nsymm = int(s[-1])
break
except:
assert "Error in read case.outputs"
#f.readline()
#f.readline()
for i in range(self.nsymm):
while 1:
s = f.readline().strip().split()
if s[0] == "Symmetry":
break
## read symmcartesian
symmt = numpy.zeros((3, 3), dtype=float)
taut = numpy.zeros((3), dtype=float)
for ir in range(3):
s = f.readline().strip().split()
for ic in range(3):
symmt[ir, ic] = float(s[ic])
s = f.readline().strip().split()
for ir in range(3):
taut[ir] = float(s[ir])
self.symmcartesian.append(symmt)
self.taucartesian.append(taut)
##read symm
symmt = numpy.zeros((3, 3), dtype=float)
taut = numpy.zeros((3), dtype=float)
for ir in range(3):
s = f.readline().strip().split()
for ic in range(3):
symmt[ir, ic] = numpy.float(s[ic])
taut[ir] = numpy.float(s[3])
self.symm.append(symmt)
self.tau.append(taut)
# end
f.close()
def checksymmxyz(self):
''' This is to check symm in cartesian coordinates and in primitive cell lattice
For details of symm, should check wien/SRC_symmetry/, latsym.f, pglsym.f,pgbsym.f
In general, if a lattice is orthorhombic, then symmcartesian is the same as symmprimitive
(at most different by a transpose (or inversion)). If a lattice is not orthorhombic, these two symms could
be related by bravais matrix.
One special case is CXZ lattice. In Wien CXZ lattice contains two cases: orthorhombic, and
monoclinic (with only gamma not equal to 90). For CXZ lattice monoclinic, symmcartesian is also the same as
symmprimitive (or with transpose (or inversion)).
'''
ortho = numpy.abs(numpy.array(self.latticeangle) - 90.0).sum() <= 1e-6
for i in range(self.nsymm):
mat = self.symmcartesian[i].transpose() # according to wien2k, why?
res = mat
if (not ortho) and (self.latticetype != "CXZ") :
bm = self.bravaismatrix
bminv = numpy.linalg.inv(bm)#.transpose()
res = numpy.dot(bminv, numpy.dot(mat, bm))
res -= self.symm[i]
print i, numpy.abs(res).sum()
## primitive cell vectors.+