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

Restore everything which was lost in rebase.

Browse Source
This commit is contained in:
Manuel Zingl 2014-12-19 11:53:06 +01:00
parent fd7ab682a4
commit b3b199bf40
6 changed files with 224 additions and 300 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

191
python/converters/wien2k_converter.py
View File

@ -54,7 +54,7 @@ class Wien2kConverter(ConverterTools):
self.bands_subgrp = bands_subgrp self.bands_subgrp = bands_subgrp
self.transp_subgrp = transp_subgrp self.transp_subgrp = transp_subgrp
self.fortran_to_replace = {'D':'E'} self.fortran_to_replace = {'D':'E'}
self.vel_file = filename+'.pmat' self.pmat_file = filename+'.pmat'
self.outputs_file = filename+'.outputs' self.outputs_file = filename+'.outputs'
self.struct_file = filename+'.struct' self.struct_file = filename+'.struct'
self.oubwin_file = filename+'.oubwin' self.oubwin_file = filename+'.oubwin'
@ -373,166 +373,137 @@ class Wien2kConverter(ConverterTools):
Reads the input files necessary for transport calculations Reads the input files necessary for transport calculations
and stores the data in the HDFfile and stores the data in the HDFfile
""" """
#Read and write files only on the master node
if not (mpi.is_master_node()): return if not (mpi.is_master_node()): return
# Check if SP, SO and n_k are already in h5 # Check if SP, SO and n_k are already in h5
ar = HDFArchive(self.hdf_file, 'a') ar = HDFArchive(self.hdf_file, 'a')
if not (self.dft_subgrp in ar): raise IOError, "No %s subgroup in hdf file found! Call convert_dmft_input first." %self.dft_subgrp if not (self.dft_subgrp in ar): raise IOError, "convert_transport_input: No %s subgroup in hdf file found! Call convert_dmft_input first." %self.dft_subgrp
SP = ar[self.dft_subgrp]['SP'] SP = ar[self.dft_subgrp]['SP']
SO = ar[self.dft_subgrp]['SO'] SO = ar[self.dft_subgrp]['SO']
n_k = ar[self.dft_subgrp]['n_k'] n_k = ar[self.dft_subgrp]['n_k']
del ar del ar
# Read relevant data from .pmat file # Read relevant data from .pmat/up/dn files
############################################ ###########################################
# band_window_optics: Contains the index of the lowest and highest band within the
# band window (used by optics) for each k-point.
# velocities_k: velocity (momentum) matrix elements between all bands in band_window_optics
# and each k-point.
if (SP == 0 or SO == 1): if (SP == 0 or SO == 1):
files = [self.vel_file] files = [self.pmat_file]
elif SP == 1: elif SP == 1:
files = [self.vel_file+'up', self.vel_file+'dn'] files = [self.pmat_file+'up', self.pmat_file+'dn']
else: # SO and SP can't both be 1 else: # SO and SP can't both be 1
assert 0, "convert_transport_input: reading velocity file error! Check SP and SO!" assert 0, "convert_transport_input: Reading velocity file error! Check SP and SO!"
vk = [[] for f in files]
kp = [[] for f in files]
bandwin_opt = []
velocities_k = [[] for f in files]
band_window_optics = []
for isp, f in enumerate(files): for isp, f in enumerate(files):
bandwin_opt_isp = [] if not os.path.exists(f) : raise IOError, "convert_transport_input: File %s does not exist" %f
if not os.path.exists(f) : raise IOError, "File %s does not exist" %f print "Reading input from %s..."%f,
print "Reading input from %s..."%f
with open(f) as R:
while 1:
try:
s = R.readline()
if (s == ''):
break
except:
break
try:
[k, nu1, nu2] = [int(x) for x in s.strip().split()]
bandwin_opt_isp.append((nu1,nu2))
dim = nu2 - nu1 +1
v_xyz = numpy.zeros((dim,dim,3), dtype = complex)
temp = R.readline().strip().split()
kp[isp].append(numpy.array([float(t) for t in temp[0:3]]))
for nu_i in xrange(dim):
for nu_j in xrange(nu_i, dim):
for i in xrange(3):
s = R.readline().strip("\n ()").split(',')
v_xyz[nu_i][nu_j][i] = float(s[0]) + float(s[1])*1j
if (nu_i != nu_j):
v_xyz[nu_j][nu_i][i] = v_xyz[nu_i][nu_j][i].conjugate()
vk[isp].append(v_xyz)
except IOError:
raise "Wien2k_converter : reading file %s failed" %self.vel_file
bandwin_opt.append(numpy.array(bandwin_opt_isp))
print "Read in %s file done!" %self.vel_file
R = ConverterTools.read_fortran_file(self, f, {'D':'E','(':'',')':'',',':''})
band_window_optics_isp = []
for ik in xrange(n_k):
R.next()
nu1 = int(R.next())
nu2 = int(R.next())
band_window_optics_isp.append((nu1, nu2))
n_bands = nu2 - nu1 + 1
velocity_xyz = numpy.zeros((n_bands, n_bands, 3), dtype = complex)
for _ in range(4): R.next()
for nu_i in range(n_bands):
for nu_j in range(nu_i, n_bands):
for i in range(3):
velocity_xyz[nu_i][nu_j][i] = R.next() + R.next()*1j
if (nu_i != nu_j): velocity_xyz[nu_j][nu_i][i] = velocity_xyz[nu_i][nu_j][i].conjugate()
velocities_k[isp].append(velocity_xyz)
band_window_optics.append(numpy.array(band_window_optics_isp))
print "DONE!"
# Read relevant data from .struct file # Read relevant data from .struct file
############################################ ######################################
if not (os.path.exists(self.struct_file)) : raise IOError, "File %s does not exist" %self.struct_file # lattice_type: bravais lattice type as defined by Wien2k
print "Reading input from %s..."%self.struct_file # lattice_constants: unit cell parameters in a. u.
# lattice_angles: unit cell angles in rad
with open(self.struct_file) as f: if not (os.path.exists(self.struct_file)) : raise IOError, "convert_transport_input: File %s does not exist" %self.struct_file
print "Reading input from %s..."%self.struct_file,
with open(self.struct_file) as R:
try: try:
f.readline() #title R.readline()
temp = f.readline() #lattice lattice_type = R.readline().split()[0]
latticetype = temp.split()[0] R.readline()
f.readline() temp = R.readline().strip().split()
temp = f.readline().strip().split() # lattice constants lattice_constants = numpy.array([float(t) for t in temp[0:3]])
latticeconstants = numpy.array([float(t) for t in temp[0:3]]) lattice_angles = numpy.array([float(t) for t in temp[3:6]]) * numpy.pi / 180.0
latticeangles = numpy.array([float(t) for t in temp[3:6]])
latticeangles *= numpy.pi/180.0
print 'Lattice: ', latticetype
print 'Lattice constants: ', latticeconstants
print 'Lattice angles: ', latticeangles
except IOError: except IOError:
raise "Wien2k_converter : reading file %s failed" %self.struct_file raise "convert_transport_input: reading file %s failed" %self.struct_file
print "DONE!"
print "Read in %s file done!" %self.struct_file
# Read relevant data from .outputs file # Read relevant data from .outputs file
############################################ #######################################
if not (os.path.exists(self.outputs_file)) : raise IOError, "File %s does not exist" %self.outputs_file # rot_symmetries: matrix representation of all (space group) symmetry operations
print "Reading input from %s..."%self.outputs_file
symmcartesian = [] if not (os.path.exists(self.outputs_file)) : raise IOError, "convert_transport_input: File %s does not exist" %self.outputs_file
taucartesian = [] print "Reading input from %s..."%self.outputs_file,
with open(self.outputs_file) as f: rot_symmetries = []
with open(self.outputs_file) as R:
try: try:
while 1: while 1:
temp = f.readline().strip(' ').split() temp = R.readline().strip(' ').split()
if (temp[0] =='PGBSYM:'): if (temp[0] =='PGBSYM:'):
nsymm = int(temp[-1]) n_symmetries = int(temp[-1])
break break
for i in range(nsymm): for i in range(n_symmetries):
while 1: while 1:
temp = f.readline().strip().split() if (R.readline().strip().split()[0] == 'Symmetry'): break
if (temp[0] == 'Symmetry'): sym_i = numpy.zeros((3, 3), dtype = float)
break
# read cartesian symmetries
symmt = numpy.zeros((3, 3), dtype = float)
taut = numpy.zeros(3, dtype = float)
for ir in range(3): for ir in range(3):
temp = f.readline().strip().split() temp = R.readline().strip().split()
for ic in range(3): for ic in range(3):
symmt[ir, ic] = float(temp[ic]) sym_i[ir, ic] = float(temp[ic])
temp = f.readline().strip().split() R.readline()
for ir in range(3): rot_symmetries.append(sym_i)
taut[ir] = float(temp[ir])
symmcartesian.append(symmt)
taucartesian.append(taut)
except IOError: except IOError:
raise "Wien2k_converter: reading file %s failed" %self.outputs_file raise "convert_transport_input: reading file %s failed" %self.outputs_file
print "DONE!"
print "Read in %s file done!" %self.outputs_file # Read relevant data from .oubwin/up/dn files
###############################################
# band_window: Contains the index of the lowest and highest band within the
# Read relevant data from .oubwin/up/down files # projected subspace (used by dmftproj) for each k-point.
############################################
if (SP == 0 or SO == 1): if (SP == 0 or SO == 1):
files = [self.oubwin_file] files = [self.oubwin_file]
elif SP == 1: elif SP == 1:
files = [self.oubwin_file+'up', self.oubwin_file+'dn'] files = [self.oubwin_file+'up', self.oubwin_file+'dn']
else: # SO and SP can't both be 1 else: # SO and SP can't both be 1
assert 0, "Reding oubwin error! Check SP and SO!" assert 0, "convert_transport_input: Reding oubwin error! Check SP and SO!"
bandwin = [numpy.zeros((n_k, 2), dtype=int) for isp in range(SP + 1 - SO)]
band_window = [numpy.zeros((n_k, 2), dtype=int) for isp in range(SP + 1 - SO)]
for isp, f in enumerate(files): for isp, f in enumerate(files):
if not os.path.exists(f): raise IOError, "File %s does not exist" %f if not os.path.exists(f): raise IOError, "convert_transport_input: File %s does not exist" %f
print "Reading input from %s..."%f print "Reading input from %s..."%f,
R = ConverterTools.read_fortran_file(self, f, self.fortran_to_replace)
assert int(R.next()) == n_k, "Wien2k_converter: Number of k-points is inconsistent in oubwin file!"
assert int(R.next()) == SO, "Wien2k_converter: SO is inconsistent in oubwin file!"
R = ConverterTools.read_fortran_file(self, f, self.fortran_to_replace)
assert int(R.next()) == n_k, "convert_transport_input: Number of k-points is inconsistent in oubwin file!"
assert int(R.next()) == SO, "convert_transport_input: SO is inconsistent in oubwin file!"
for ik in xrange(n_k): for ik in xrange(n_k):
R.next() R.next()
bandwin[isp][ik,0] = R.next() band_window[isp][ik,0] = R.next()
bandwin[isp][ik,1] = R.next() band_window[isp][ik,1] = R.next()
R.next() R.next()
print "DONE!"
print "Read in %s files done!" %self.oubwin_file
# Put data to HDF5 file # Put data to HDF5 file
ar = HDFArchive(self.hdf_file, 'a') ar = HDFArchive(self.hdf_file, 'a')
if not (self.transp_subgrp in ar): ar.create_group(self.transp_subgrp) if not (self.transp_subgrp in ar): ar.create_group(self.transp_subgrp)
# The subgroup containing the data. If it does not exist, it is created. If it exists, the data is overwritten!!! # The subgroup containing the data. If it does not exist, it is created. If it exists, the data is overwritten!!!
things_to_save = ['bandwin', 'bandwin_opt', 'kp', 'vk', 'latticetype', 'latticeconstants', 'latticeangles', 'nsymm', 'symmcartesian', things_to_save = ['band_window', 'band_window_optics', 'velocities_k', 'lattice_type', 'lattice_constants', 'lattice_angles', 'n_symmetries', 'rot_symmetries']
'taucartesian']
for it in things_to_save: ar[self.transp_subgrp][it] = locals()[it] for it in things_to_save: ar[self.transp_subgrp][it] = locals()[it]
del ar del ar

303
python/sumk_dft_tools.py
View File

@ -1,4 +1,3 @@
################################################################################ ################################################################################
# #
# TRIQS: a Toolbox for Research in Interacting Quantum Systems # TRIQS: a Toolbox for Research in Interacting Quantum Systems
@ -501,72 +500,72 @@ class SumkDFTTools(SumkDFT):
""" """
Reads the data for transport calculations from the HDF file Reads the data for transport calculations from the HDF file
""" """
thingstoread = ['band_window','band_window_optics','velocities_k','lattice_angles','lattice_constants','lattice_type','n_symmetries','rot_symmetries']
thingstoread = ['bandwin','bandwin_opt','kp','latticeangles','latticeconstants','latticetype','nsymm','symmcartesian','vk']
retval = self.read_input_from_hdf(subgrp=self.transp_data,things_to_read = thingstoread) retval = self.read_input_from_hdf(subgrp=self.transp_data,things_to_read = thingstoread)
return retval return retval
def cellvolume(self, latticetype, latticeconstants, latticeangle): def cellvolume(self, lattice_type, lattice_constants, latticeangle):
""" """
Calculate cell volume: volumecc conventional cell, volumepc, primitive cell. Calculate cell volume: volumecc conventional cell, volumepc, primitive cell.
""" """
a = latticeconstants[0] a = lattice_constants[0]
b = latticeconstants[1] b = lattice_constants[1]
c = latticeconstants[2] c = lattice_constants[2]
c_al = numpy.cos(latticeangle[0]) c_al = numpy.cos(latticeangle[0])
c_be = numpy.cos(latticeangle[1]) c_be = numpy.cos(latticeangle[1])
c_ga = numpy.cos(latticeangle[2]) 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) 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} det = {"P":1, "F":4, "B":2, "R":3, "H":1, "CXY":2, "CYZ":2, "CXZ":2}
volumepc = volumecc / det[latticetype] volumepc = volumecc / det[lattice_type]
return volumecc, volumepc return volumecc, volumepc
def transport_distribution(self, dir_list=[(0,0)], broadening=0.01, energywindow=None, Om_mesh=[0.0], beta=40, with_Sigma=False, n_om=None): def transport_distribution(self, directions=['xx'], energy_window=None, Om_mesh=[0.0], beta=40, with_Sigma=False, n_om=None, broadening=0.01):
"""calculate Tr A(k,w) v(k) A(k, w+q) v(k) and optics. """
energywindow: regime for omega integral calculate Tr A(k,w) v(k) A(k, w+Om) v(k).
Om_mesh: contains the frequencies of the optic conductivitity. Om_mesh is repinned to the self-energy mesh energy_window: regime for omega integral
(hence exact values might be different from those given in Om_mesh) Om_mesh: mesh for optic conductivitity. Om_mesh is repinned to the self-energy mesh!
dir_list: list to defines the indices of directions. xx,yy,zz,xy,yz,zx. directions: list of directions: xx,yy,zz,xy,yz,zx.
((0, 0) --> xx, (1, 1) --> yy, (0, 2) --> xz, default: (0, 0)) with_Sigma: Use Sigma_w = 0 if False (In this case it is necessary to specifiy the energywindow (energy_window),
with_Sigma: Use Sigma = 0 if False the number of omega points (n_om) in the window and the broadening (broadening)).
""" """
# Check if wien converter was called # Check if wien converter was called and read transport subgroup form hdf file
if mpi.is_master_node(): if mpi.is_master_node():
ar = HDFArchive(self.hdf_file, 'a') ar = HDFArchive(self.hdf_file, 'a')
if not (self.transp_data in ar): raise IOError, "No %s subgroup in hdf file found! Call convert_transp_input first." %self.transp_data if not (self.transp_data in ar): raise IOError, "transport_distribution: No %s subgroup in hdf file found! Call convert_transp_input first." %self.transp_data
self.dir_list = dir_list
self.read_transport_input_from_hdf() self.read_transport_input_from_hdf()
velocities = self.vk
n_inequiv_spin_blocks = self.SP + 1 - self.SO # up and down are equivalent if SP = 0
n_inequiv_spin_blocks = self.SP + 1 - self.SO # up and down are equivalent if SP = 0
self.directions = directions
dir_to_int = {'x':0, 'y':1, 'z':2}
# k-dependent-projections.
assert self.k_dep_projection == 1, "transport_distribution: k dependent projection is not implemented!"
# calculate A(k,w) # calculate A(k,w)
####################################### #######################################
# use k-dependent-projections. # Define mesh for Greens function and in the specified energy window
assert self.k_dep_projection == 1, "Not implemented!"
# Define mesh for Greens function and the used energy range
if (with_Sigma == True): if (with_Sigma == True):
self.omega = numpy.array([round(x.real,12) for x in self.Sigma_imp_w[0].mesh]) self.omega = numpy.array([round(x.real,12) for x in self.Sigma_imp_w[0].mesh])
mesh = None
mu = self.chemical_potential mu = self.chemical_potential
n_om = len(self.omega) n_om = len(self.omega)
print "Using omega mesh provided by Sigma." print "Using omega mesh provided by Sigma!"
if energywindow is not None: if energy_window is not None:
# Find according window in Sigma mesh # Find according window in Sigma mesh
ioffset = numpy.sum(self.omega < energywindow[0]) ioffset = numpy.sum(self.omega < energy_window[0])
self.omega = self.omega[numpy.logical_and(self.omega >= energywindow[0], self.omega <= energywindow[1])] self.omega = self.omega[numpy.logical_and(self.omega >= energy_window[0], self.omega <= energy_window[1])]
n_om = len(self.omega) n_om = len(self.omega)
# Truncate Sigma to given omega window # Truncate Sigma to given omega window
# In the future there should be an option in gf to manipulate the mesh (e.g. truncate) directly.
# For we stick with this:
for icrsh in range(self.n_corr_shells): for icrsh in range(self.n_corr_shells):
Sigma_save = self.Sigma_imp_w[icrsh].copy() Sigma_save = self.Sigma_imp_w[icrsh].copy()
spn = self.spin_block_names[self.corr_shells[icrsh]['SO']] spn = self.spin_block_names[self.corr_shells[icrsh]['SO']]
@ -576,179 +575,133 @@ class SumkDFTTools(SumkDFT):
for iL in g.indices: for iL in g.indices:
for iR in g.indices: for iR in g.indices:
for iom in xrange(n_om): for iom in xrange(n_om):
g.data[iom,iL,iR] = Sigma_save[i].data[ioffset+iom,iL,iR] # FIXME IS THIS CLEAN? MANIPULATING data DIRECTLY? g.data[iom,iL,iR] = Sigma_save[i].data[ioffset+iom,iL,iR]
else: else:
assert n_om is not None, "Number of omega points (n_om) needed!" assert n_om is not None, "transport_distribution: Number of omega points (n_om) needed to calculate transport distribution!"
assert energywindow is not None, "Energy window needed!" assert energy_window is not None, "transport_distribution: Energy window needed to calculate transport distribution!"
self.omega = numpy.linspace(energywindow[0],energywindow[1],n_om) assert broadening != 0.0 and broadening is not None, "transport_distribution: Broadening necessary to calculate transport distribution!"
mu = 0.0 self.omega = numpy.linspace(energy_window[0],energy_window[1],n_om)
mesh = [energy_window[0], energy_window[1], n_om]
mu = 0.0
# Check if energy_window is sufficiently large
if (abs(self.fermi_dis(self.omega[0]*beta)*self.fermi_dis(-self.omega[0]*beta)) > 1e-5 if (abs(self.fermi_dis(self.omega[0]*beta)*self.fermi_dis(-self.omega[0]*beta)) > 1e-5
or abs(self.fermi_dis(self.omega[-1]*beta)*self.fermi_dis(-self.omega[-1]*beta)) > 1e-5): or abs(self.fermi_dis(self.omega[-1]*beta)*self.fermi_dis(-self.omega[-1]*beta)) > 1e-5):
print "\n##########################################" print "\n####################################################################"
print "WARNING: Energywindow might be too narrow!" print "transport_distribution: WARNING - energy window might be too narrow!"
print "##########################################\n" print "####################################################################\n"
# Define mesh for optic conductivity
d_omega = round(numpy.abs(self.omega[0] - self.omega[1]), 12) d_omega = round(numpy.abs(self.omega[0] - self.omega[1]), 12)
iOm_mesh = numpy.array([int(Om / d_omega) for Om in Om_mesh])
# define exact mesh for optic conductivity self.Om_mesh = iOm_mesh * d_omega
Om_mesh_ex = numpy.array([int(x / d_omega) for x in Om_mesh])
self.Om_meshr= Om_mesh_ex*d_omega
if mpi.is_master_node(): if mpi.is_master_node():
print "Chemical potential: ", mu print "Chemical potential: ", mu
print "Using n_om = %s points in the energywindow [%s,%s]"%(n_om, self.omega[0], self.omega[-1]) print "Using n_om = %s points in the energy_window [%s,%s]"%(n_om, self.omega[0], self.omega[-1]),
print "omega vector is:" print "where the omega vector is:"
print self.omega print self.omega
print "Omega mesh interval ", d_omega print "Calculation requested for Omega mesh: ", numpy.array(Om_mesh)
print "Provided Om_mesh ", numpy.array(Om_mesh) print "Omega mesh automatically repinned to: ", self.Om_mesh
print "Pinnend Om_mesh to ", self.Om_meshr
# output P(\omega)_xy should have the same dimension as defined in mshape. self.Gamma_w = {direction: numpy.zeros((len(self.Om_mesh), n_om), dtype=numpy.float_) for direction in self.directions}
self.Pw_optic = numpy.zeros((len(dir_list), len(Om_mesh_ex), n_om), dtype=numpy.float_)
ik = 0
spn = self.spin_block_names[self.SO]
ntoi = self.spin_names_to_ind[self.SO]
G_w = BlockGf(name_block_generator=[(spn[isp], GfReFreq(indices=range(self.n_orbitals[ik][isp]), window=(self.omega[0], self.omega[-1]), n_points = n_om))
for isp in range(n_inequiv_spin_blocks) ], make_copies=False)
mupat = [numpy.identity(self.n_orbitals[ik][isp], numpy.complex_) * mu for isp in range(n_inequiv_spin_blocks)] # construct mupat
Annkw = [numpy.zeros((self.n_orbitals[ik][isp], self.n_orbitals[ik][isp], n_om), dtype=numpy.complex_) for isp in range(n_inequiv_spin_blocks)]
# Sum over all k-points
ikarray = numpy.array(range(self.n_k)) ikarray = numpy.array(range(self.n_k))
for ik in mpi.slice_array(ikarray): for ik in mpi.slice_array(ikarray):
unchangedsize = all([ self.n_orbitals[ik][isp] == mupat[isp].shape[0] for isp in range(n_inequiv_spin_blocks)]) # Calculate G_w for ik and initialize A_kw
if not unchangedsize: G_w = self.lattice_gf(ik, mu, iw_or_w="w", beta=beta, broadening=broadening, mesh=mesh, with_Sigma=with_Sigma)
# recontruct green functions. A_kw = [numpy.zeros((self.n_orbitals[ik][isp], self.n_orbitals[ik][isp], n_om), dtype=numpy.complex_)
G_w = BlockGf(name_block_generator=[(spn[isp], GfReFreq(indices=range(self.n_orbitals[ik][isp]), window = (self.omega[0], self.omega[-1]), n_points = n_om)) for isp in range(n_inequiv_spin_blocks)]
for isp in range(n_inequiv_spin_blocks) ], make_copies=False)
# construct mupat
mupat = [numpy.identity(self.n_orbitals[ik][isp], numpy.complex_) * mu for isp in range(n_inequiv_spin_blocks)]
#set a temporary array storing spectral functions with band index. Note, usually we should have spin index
Annkw = [numpy.zeros((self.n_orbitals[ik][isp], self.n_orbitals[ik][isp], n_om), dtype=numpy.complex_) for isp in range(n_inequiv_spin_blocks)]
# get lattice green function
G_w << 1*Omega + 1j*broadening for isp in range(n_inequiv_spin_blocks):
# Obtain A_kw from G_w (swapaxes is used to have omega in the 3rd dimension)
A_kw[isp].real = -copy.deepcopy(G_w[self.spin_block_names[self.SO][isp]].data.swapaxes(0,1).swapaxes(1,2)).imag / numpy.pi
b_min = max(self.band_window[isp][ik, 0], self.band_window_optics[isp][ik, 0])
b_max = min(self.band_window[isp][ik, 1], self.band_window_optics[isp][ik, 1])
A_i = slice(b_min - self.band_window[isp][ik, 0], b_max - self.band_window[isp][ik, 0] + 1)
v_i = slice(b_min - self.band_window_optics[isp][ik, 0], b_max - self.band_window_optics[isp][ik, 0] + 1)
MS = copy.deepcopy(mupat) # loop over all symmetries
for ibl in range(n_inequiv_spin_blocks): for R in self.rot_symmetries:
ind = ntoi[spn[ibl]] # get transformed velocity under symmetry R
n_orb = self.n_orbitals[ik][ibl] vel_R = copy.deepcopy(self.velocities_k[isp][ik])
MS[ibl] = self.hopping[ik,ind,0:n_orb,0:n_orb].real - mupat[ibl] for nu1 in range(self.band_window_optics[isp][ik, 1] - self.band_window_optics[isp][ik, 0] + 1):
G_w -= MS for nu2 in range(self.band_window_optics[isp][ik, 1] - self.band_window_optics[isp][ik, 0] + 1):
vel_R[nu1][nu2][:] = numpy.dot(R, vel_R[nu1][nu2][:])
if (with_Sigma == True): # calculate Gamma_w for each direction from the velocities vel_R and the spectral function A_kw
tmp = G_w.copy() # init temporary storage for direction in self.directions:
# form self energy from impurity self energy and double counting term.
sigma_minus_dc = self.add_dc(iw_or_w="w")
# substract self energy
for icrsh in range(self.n_corr_shells):
for sig, gf in tmp: tmp[sig] << self.upfold(ik, icrsh, sig, sigma_minus_dc[icrsh][sig], gf)
G_w -= tmp
G_w.invert()
for isp in range(n_inequiv_spin_blocks):
Annkw[isp].real = -copy.deepcopy(G_w[self.spin_block_names[self.SO][isp]].data.swapaxes(0,1).swapaxes(1,2)).imag / numpy.pi
for isp in range(n_inequiv_spin_blocks):
kvel = velocities[isp][ik]
Pwtem = numpy.zeros((len(dir_list), len(Om_mesh_ex), n_om), dtype=numpy.float_)
bmin = max(self.bandwin[isp][ik, 0], self.bandwin_opt[isp][ik, 0])
bmax = min(self.bandwin[isp][ik, 1], self.bandwin_opt[isp][ik, 1])
Astart = bmin - self.bandwin[isp][ik, 0]
Aend = bmax - self.bandwin[isp][ik, 0] + 1
vstart = bmin - self.bandwin_opt[isp][ik, 0]
vend = bmax - self.bandwin_opt[isp][ik, 0] + 1
#symmetry loop
for Rmat in self.symmcartesian:
# get new velocity.
Rkvel = copy.deepcopy(kvel)
for vnb1 in range(self.bandwin_opt[isp][ik, 1] - self.bandwin_opt[isp][ik, 0] + 1):
for vnb2 in range(self.bandwin_opt[isp][ik, 1] - self.bandwin_opt[isp][ik, 0] + 1):
Rkvel[vnb1][vnb2][:] = numpy.dot(Rmat, Rkvel[vnb1][vnb2][:])
ipw = 0
for (ir, ic) in dir_list:
for iw in xrange(n_om): for iw in xrange(n_om):
for iq in range(len(Om_mesh_ex)): for iq in range(len(self.Om_mesh)):
if(iw + Om_mesh_ex[iq] >= n_om): if(iw + iOm_mesh[iq] >= n_om): continue
continue self.Gamma_w[direction][iq, iw] += (numpy.dot(numpy.dot(numpy.dot(vel_R[v_i, v_i, dir_to_int[direction[0]]],
A_kw[isp][A_i, A_i, iw]), vel_R[v_i, v_i, dir_to_int[direction[1]]]),
A_kw[isp][A_i, A_i, iw + iOm_mesh[iq]]).trace().real * self.bz_weights[ik])
# construct matrix for A and velocity. for direction in self.directions:
Annkwl = Annkw[isp][Astart:Aend, Astart:Aend, iw] self.Gamma_w[direction] = (mpi.all_reduce(mpi.world, self.Gamma_w[direction], lambda x, y : x + y)
Annkwr = Annkw[isp][Astart:Aend, Astart:Aend, iw + Om_mesh_ex[iq]] / self.cellvolume(self.lattice_type, self.lattice_constants, self.lattice_angles)[1]) / self.n_symmetries
Rkveltr = Rkvel[vstart:vend, vstart:vend, ir]
Rkveltc = Rkvel[vstart:vend, vstart:vend, ic]
# print Annkwl.shape, Annkwr.shape, Rkveltr.shape, Rkveltc.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 = mpi.all_reduce(mpi.world, self.Pw_optic, lambda x, y : x + y) def transport_coefficient(self, direction, iq=0, n=0, beta=40):
self.Pw_optic *= (2 - self.SP) """
calculates the transport coefficients A_n in a given direction and for a given Omega. (see documentation)
A_1 is set to nan if requested for Omega != 0.0
iq: index of Omega point in Om_mesh
direction: 'xx','yy','zz','xy','xz','yz'
"""
if not (mpi.is_master_node()): return
assert hasattr(self,'Gamma_w'), "transport_coefficient: Run transport_distribution first or load data from h5!"
A = 0.0
omegaT = self.omega * beta
d_omega = self.omega[1] - self.omega[0]
if (self.Om_mesh[iq] == 0.0):
for iw in xrange(self.Gamma_w[direction].shape[1]):
A += self.Gamma_w[direction][iq, iw] * self.fermi_dis(omegaT[iw]) * self.fermi_dis(-omegaT[iw]) * numpy.float(omegaT[iw])**n * d_omega
elif (n == 0.0):
for iw in xrange(self.Gamma_w[direction].shape[1]):
A += (self.Gamma_w[direction][iq, iw] * (self.fermi_dis(omegaT[iw]) - self.fermi_dis(omegaT[iw] + self.Om_mesh[iq] * beta))
/ (self.Om_mesh[iq] * beta) * d_omega)
else:
A = numpy.nan
return A * numpy.pi * (2.0-self.SP)
def conductivity_and_seebeck(self, beta=40): def conductivity_and_seebeck(self, beta=40):
""" #return 1/T*A0, that is Conductivity in unit 1/V
calculate, save and return Conductivity
""" """
Calculates the Seebeck coefficient and the conductivity for a given Gamma_w
"""
if not (mpi.is_master_node()): return
if mpi.is_master_node(): assert hasattr(self,'Gamma_w'), "conductivity_and_seebeck: Run transport_distribution first or load data from h5!"
assert hasattr(self,'Pw_optic'), "Run transport_distribution first or load data from h5!" n_q = self.Gamma_w[self.directions[0]].shape[0]
assert hasattr(self,'latticetype'), "Run convert_transp_input first or load data from h5!"
volcc, volpc = self.cellvolume(self.latticetype, self.latticeconstants, self.latticeangles) A0 = {direction: numpy.full((n_q,),numpy.nan) for direction in self.directions}
A1 = {direction: numpy.full((n_q,),numpy.nan) for direction in self.directions}
self.seebeck = {direction: numpy.nan for direction in self.directions}
self.optic_cond = {direction: numpy.full((n_q,),numpy.nan) for direction in self.directions}
n_direction, n_q, n_w= self.Pw_optic.shape for direction in self.directions:
omegaT = self.omega * beta for iq in xrange(n_q):
A0 = numpy.zeros((n_direction,n_q), dtype=numpy.float_) A0[direction][iq] = self.transport_coefficient(direction, iq=iq, n=0, beta=beta)
q_0 = False A1[direction][iq] = self.transport_coefficient(direction, iq=iq, n=1, beta=beta)
self.seebeck = numpy.zeros((n_direction,), dtype=numpy.float_) print "A_0 in direction %s for Omega = %.2f %e a.u." % (direction, self.Om_mesh[iq], A0[direction][iq])
self.seebeck[:] = numpy.NAN print "A_1 in direction %s for Omega = %.2f %e a.u." % (direction, self.Om_mesh[iq], A1[direction][iq])
if ~numpy.isnan(A1[direction][iq]):
d_omega = self.omega[1] - self.omega[0] # Seebeck is overwritten if there is more than one Omega = 0 in Om_mesh
for iq in xrange(n_q): self.seebeck[direction] = - A1[direction][iq] / A0[direction][iq] * 86.17
# treat q = 0, caclulate conductivity and seebeck self.optic_cond[direction] = A0[direction] * 10700.0
if (self.Om_meshr[iq] == 0.0): for iq in xrange(n_q):
# if Om_meshr contains multiple entries with w=0, A1 is overwritten! print "Conductivity in direction %s for Omega = %.2f %f x 10^4 Ohm^-1 cm^-1" % (direction, self.Om_mesh[iq], self.optic_cond[direction][iq])
q_0 = True if not (numpy.isnan(A1[direction][iq])):
A1 = numpy.zeros((n_direction,), dtype=numpy.float_) print "Seebeck in direction %s for Omega = 0.00 %f x 10^(-6) V/K" % (direction, self.seebeck[direction])
for idir in range(n_direction):
for iw in xrange(n_w):
A0[idir, iq] += beta * self.Pw_optic[idir, iq, iw] * self.fermi_dis(omegaT[iw]) * self.fermi_dis(-omegaT[iw])
A1[idir] += beta * self.Pw_optic[idir, iq, iw] * self.fermi_dis(omegaT[iw]) * self.fermi_dis(-omegaT[iw]) * numpy.float(omegaT[iw])
self.seebeck[idir] = -A1[idir] / A0[idir, iq]
print "A0", A0[idir, iq] *d_omega/beta
print "A1", A1[idir] *d_omega/beta
# treat q ~= 0, calculate optical conductivity
else:
for idir in range(n_direction):
for iw in xrange(n_w):
A0[idir, iq] += self.Pw_optic[idir, iq, iw] * (self.fermi_dis(omegaT[iw]) - self.fermi_dis(omegaT[iw] + self.Om_meshr[iq] * beta)) / self.Om_meshr[iq]
A0 *= d_omega
#cond = beta * self.tdintegral(beta, 0)[index]
print "V in bohr^3 ", volpc
# transform to standard unit as in resistivity
self.optic_cond = A0 * 10700.0 / volpc
self.seebeck *= 86.17
# print
for im in range(n_direction):
for iq in xrange(n_q):
print "Conductivity in direction %s for Om_mesh %d %.4f x 10^4 Ohm^-1 cm^-1" % (self.dir_list[im], iq, self.optic_cond[im, iq])
print "Resistivity in direction %s for Om_mesh %d %.4f x 10^-4 Ohm cm" % (self.dir_list[im], iq, 1.0 / self.optic_cond[im, iq])
if q_0:
print "Seebeck in direction %s for q = 0 %.4f x 10^(-6) V/K" % (self.dir_list[im], self.seebeck[im])
def fermi_dis(self, x): def fermi_dis(self, x):
"""
fermi distribution at x = omega * beta
"""
return 1.0/(numpy.exp(x)+1) return 1.0/(numpy.exp(x)+1)

BIN
test/SrVO3.h5
View File

Binary file not shown.

BIN
test/SrVO3_Sigma.h5
View File

Binary file not shown.

BIN
test/srvo3_transp.output.h5
View File

Binary file not shown.

8
test/srvo3_transp.py
View File

@ -33,13 +33,13 @@ Converter.convert_transport_input()
SK = SumkDFTTools(hdf_file='SrVO3.h5', use_dft_blocks=True) SK = SumkDFTTools(hdf_file='SrVO3.h5', use_dft_blocks=True)
ar = HDFArchive('SrVO3_Sigma.h5', 'a') ar = HDFArchive('SrVO3_Sigma.h5', 'a')
Sigma = ar['dmft_transp_output']['Sigma'] Sigma = ar['dmft_transp_output']['Sigma_w']
SK.put_Sigma(Sigma_imp = [Sigma]) SK.put_Sigma(Sigma_imp = [Sigma])
del ar del ar
SK.transport_distribution(dir_list=[(0,0)], broadening=0.0, energywindow=[-0.3,0.3], Om_mesh=[0.00, 0.02] , beta=beta, with_Sigma=True) SK.transport_distribution(directions=['xx'], broadening=0.0, energy_window=[-0.3,0.3], Om_mesh=[0.00, 0.02] , beta=beta, with_Sigma=True)
#SK.save(['Pw_optic','Om_meshr','omega','dir_list']) #SK.save(['Gamma_w','Om_meshr','omega','directions'])
#SK.load(['Pw_optic','Om_meshr','omega','dir_list']) #SK.load(['Gamma_w','Om_meshr','omega','directions'])
SK.conductivity_and_seebeck(beta=beta) SK.conductivity_and_seebeck(beta=beta)
SK.hdf_file = 'srvo3_transp.output.h5' SK.hdf_file = 'srvo3_transp.output.h5'
SK.save(['seebeck','optic_cond']) SK.save(['seebeck','optic_cond'])