mirror of
https://github.com/triqs/dft_tools
synced 2024-12-22 12:23:41 +01:00
[transport] Delete obsolete SumK_Transport files
This commit is contained in:
Michel Ferrero
Priyanka Seth
parent
6f6c8d1c56
commit
287c44116b
@ -1,902 +0,0 @@
|
|||||||
################################################################################
|
|
||||||
#
|
|
||||||
# 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
|
|
||||||
|
|
||||||
import pytriqs.utility.dichotomy as Dichotomy
|
|
||||||
from pytriqs.gf.local import *
|
|
||||||
import pytriqs.utility.mpi as myMPI
|
|
||||||
from datetime import datetime
|
|
||||||
from pytriqs.applications.dft.sumk_lda import *
|
|
||||||
from pytriqs.applications.dft.sumk_lda_tools import *
|
|
||||||
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
|
|
||||||
|
|
||||||
# print S[self.block_names[self.SO][0]].data
|
|
||||||
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))
|
|
||||||
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 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
|
|
||||||
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_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) 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
|
|
||||||
|
|
||||||
@ -1,431 +0,0 @@
|
|||||||
|
|
||||||
################################################################################
|
|
||||||
#
|
|
||||||
# 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.+
|
|
||||||
|
|
||||||
Loading…
Reference in New Issue
Block a user