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

remove any transport from sumk_dft_tools.py

Browse Source
This commit is contained in:
phibeck 2022-11-02 10:48:44 -04:00 committed by Alexander Hampel
parent 345fd14a39
commit d5e6d60258
2 changed files with 18 additions and 561 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

545
python/triqs_dft_tools/sumk_dft_tools.py
View File

@ -31,8 +31,6 @@ from .symmetry import *
from .sumk_dft import SumkDFT
from scipy.integrate import *
from scipy.interpolate import *
from scipy import constants as constants
from itertools import product
if not hasattr(numpy, 'full'):
# polyfill full for older numpy:
@ -1060,546 +1058,3 @@ class SumkDFTTools(SumkDFT):
f.close()
# ----------------- transport -----------------------
def read_transport_input_from_hdf(self):
r"""
Reads the data for transport calculations from the hdf5 archive.
"""
assert self.dft_code in ('wien2k','elk'), "Transport has only been implemented for wien2k and elk inputs"
thingstoread = ['band_window_optics', 'velocities_k']
self.read_input_from_hdf(
subgrp=self.transp_data, things_to_read=thingstoread)
if(self.dft_code=="wien2k"):
thingstoread = ['band_window', 'lattice_angles', 'lattice_constants',
'lattice_type', 'n_symmetries', 'rot_symmetries']
elif(self.dft_code=="elk"):
thingstoread = ['band_window', 'n_symmetries',
'rot_symmetries','cell_vol']
self.read_input_from_hdf(
subgrp=self.misc_data, things_to_read=thingstoread)
if(self.dft_code=="wien2k"):
self.cell_vol = self.cellvolume(self.lattice_type, self.lattice_constants, self.lattice_angles)[1]
def read_transport_input_from_hdf_wannier90(self):
r"""
Reads the data for transport calculations from the hdf5 archive.
"""
thingstoread = ['band_window_optics', 'nk_optics']
self.read_input_from_hdf(
subgrp=self.transp_data, things_to_read=thingstoread)
thingstoread = ['band_window', 'n_symmetries', 'rot_symmetries']
self.read_input_from_hdf(
subgrp=self.misc_data, things_to_read=thingstoread)
def cellvolume(self, lattice_type, lattice_constants, latticeangle):
r"""
Determines the conventional und primitive unit cell volumes.
Parameters
----------
lattice_type : string
Lattice type according to the Wien2k convention (P, F, B, R, H, CXY, CYZ, CXZ).
lattice_constants : list of double
Lattice constants (a, b, c).
lattice angles : list of double
Lattice angles (:math:`\alpha, \beta, \gamma`).
Returns
-------
vol_c : double
Conventional unit cell volume.
vol_p : double
Primitive unit cell volume.
"""
a = lattice_constants[0]
b = lattice_constants[1]
c = lattice_constants[2]
c_al = numpy.cos(latticeangle[0])
c_be = numpy.cos(latticeangle[1])
c_ga = numpy.cos(latticeangle[2])
vol_c = a * b * c * \
numpy.sqrt(1 + 2 * c_al * c_be * c_ga -
c_al ** 2 - c_be ** 2 - c_ga ** 2)
det = {"P": 1, "F": 4, "B": 2, "R": 3,
"H": 1, "CXY": 2, "CYZ": 2, "CXZ": 2}
vol_p = vol_c / det[lattice_type]
return vol_c, vol_p
# Uses .data of only GfReFreq objects.
def transport_distribution(self, beta, directions=['xx'], energy_window=None, Om_mesh=[0.0], with_Sigma=False, n_om=None, broadening=0.0, code='wien2k'):
r"""
Calculates the transport distribution
.. math::
\Gamma_{\alpha\beta}\left(\omega+\Omega/2, \omega-\Omega/2\right) = \frac{1}{V} \sum_k Tr\left(v_{k,\alpha}A_{k}(\omega+\Omega/2)v_{k,\beta}A_{k}\left(\omega-\Omega/2\right)\right)
in the direction :math:`\alpha\beta`. The velocities :math:`v_{k}` are read from the transport subgroup of the hdf5 archive.
Parameters
----------
beta : double
Inverse temperature :math:`\beta`.
directions : list of double, optional
:math:`\alpha\beta` e.g.: ['xx','yy','zz','xy','xz','yz'].
energy_window : list of double, optional
Specifies the upper and lower limit of the frequency integration for :math:`\Omega=0.0`. The window is automatically enlarged by the largest :math:`\Omega` value,
hence the integration is performed in the interval [energy_window[0]-max(Om_mesh), energy_window[1]+max(Om_mesh)].
Om_mesh : list of double, optional
:math:`\Omega` frequency mesh of the optical conductivity. For the conductivity and the Seebeck coefficient :math:`\Omega=0.0` has to be
part of the mesh. In the current version Om_mesh is repined to the mesh provided by the self-energy! The actual mesh is printed on the screen and stored as
member Om_mesh.
with_Sigma : boolean, optional
Determines whether the calculation is performed with or without self energy. If this parameter is set to False the self energy is set to zero (i.e. the DFT band
structure :math:`A(k,\omega)` is used). Note: For with_Sigma=False it is necessary to specify the parameters energy_window, n_om and broadening.
n_om : integer, optional
Number of equidistant frequency points in the interval [energy_window[0]-max(Om_mesh), energy_window[1]+max(Om_mesh)]. This parameters is only used if
with_Sigma = False.
broadening : double, optional
Lorentzian broadening. It is necessary to specify the boradening if with_Sigma = False, otherwise this parameter can be set to 0.0.
code : string
DFT code from which velocities are being read. Options: 'wien2k', 'wannier90'
"""
BOHRTOANG = constants.physical_constants['Bohr radius'][0]/constants.angstrom
HARTREETOEV = constants.physical_constants['Hartree energy'][0]/constants.eV
n_inequiv_spin_blocks = self.SP + 1 - self.SO
if code in ('wien2k'):
# Check if wien converter was called and read transport subgroup form
# hdf file
if mpi.is_master_node():
ar = HDFArchive(self.hdf_file, 'r')
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)
# check if outputs file was converted
if not ('n_symmetries' in ar['dft_misc_input']):
raise IOError("transport_distribution: n_symmetries missing. Check if case.outputs file is present and call convert_misc_input() or convert_dft_input().")
self.read_transport_input_from_hdf()
cell_volume = self.cellvolume(self.lattice_type, self.lattice_constants, self.lattice_angles)[1]
n_symmetries = self.n_symmetries
elif code in ('wannier90'):
# check if spin-unpolarized
assert n_inequiv_spin_blocks == 1, "Spin-polarized optical conductivity calculations not implemented with Wannier90"
# read in transport input
self.read_transport_input_from_hdf_wannier90()
# checks for right formatting of self.nk_optics
assert len(self.nk_optics) in [1,3], '"nk_optics" must be given as three integers or one float'
if len(self.nk_optics) == 1: assert np.array(list(self.nk_optics)).dtype in (int, float), '"nk_optics" single value must be float or integer'
if len(self.nk_optics) == 3: assert np.array(list(self.nk_optics)).dtype == int, '"nk_optics" mesh must be integers'
n_symmetries = 1
# calculate velocity
#wberri = wb.System_w90('/mnt/home/sbeck/Dropbox/ccqlin030/sro/I4_mmm_prim/wan_conv_12_v/sro', berry=True)
pathname = './'
seedname = 'sro'
fermi = 0.
if len(self.nk_optics) == 1:
interpolate_factor = self.nk_optics[0]
nk_x, nk_y, nk_z = list(map(lambda i: int(numpy.ceil(interpolate_factor * len(set(self.kpts[:,i])))), range(3)))
else:
nk_x, nk_y, nk_z = self.nk_optics
#nk_x, nk_y, nk_z = 10
#nk_x, nk_y, nk_z = 12
#nk_x, nk_y, nk_z = 34
n_orb = numpy.max([self.n_orbitals[ik][0] for ik in range(self.n_k)])
shift_gamma = [0.0,0.0,0.0]
#shift_gamma = [0.015,0.015,0.015]
things_to_modify = {'bz_weights': None, 'hopping': None, 'kpt_weights': None, 'kpts': None,
'n_k': None, 'n_orbitals': None, 'proj_mat': None, 'band_window': None, 'band_window_optics': None}
things_to_store = dict.fromkeys(things_to_modify, None)
# initialize variables
n_kpts = nk_x * nk_y * nk_z
kpts = numpy.zeros((n_kpts, 3))
hopping = numpy.zeros((n_kpts, 1, n_orb, n_orb), dtype=complex)
proj_mat = numpy.zeros(numpy.shape(hopping[:,0,0,0]) + numpy.shape(self.proj_mat[0,:]), dtype=complex)
if mpi.is_master_node():
print(hopping.shape, self.proj_mat.shape, numpy.shape(hopping[:,0,0,0]) + numpy.shape(self.proj_mat[0,:]))
# simple modifications
things_to_modify['n_k'] = n_kpts
things_to_modify['n_orbitals'] = numpy.full((n_kpts, 1), n_orb)
for key in ['bz_weights', 'kpt_weights']:
things_to_modify[key] = numpy.full(n_kpts, 1/n_kpts)
n_inequiv_spin_blocks = self.SP + 1 - self.SO
for key in ['band_window', 'band_window_optics']:
things_to_modify[key] = [numpy.full((n_kpts, 2), self.band_window[isp][0]) for isp in range(n_inequiv_spin_blocks)]
velocities_k = None
cell_volume = None
kpts = None
if mpi.is_master_node():
# try wannierberri import
try:
import wannierberri as wb
except ImportError:
print('ImportError: WannierBerri needs to be installed to run test "Py_w90_optics_Sr2RuO4"')
try:
mpi.MPI.COMM_WORLD.Abort(1)
except:
sys.exit()
# initialize WannierBerri system
wberri = wb.System_w90(pathname + seedname, berry=True)
grid = wb.Grid(wberri, NKdiv=1, NKFFT=[nk_x, nk_y, nk_z])
dataK = Data_K(wberri, dK=shift_gamma, grid=grid)
# construct velocities from dataK
V_H_diag = numpy.zeros(numpy.shape(dataK.V_H), dtype=complex)
V_H_diag[:, range(V_H_diag.shape[1]), range(V_H_diag.shape[1]), :] = numpy.diagonal(dataK.V_H[:,:,:,:],axis1=1, axis2=2).transpose(0,2,1).copy()
velocities_k = ( V_H_diag - dataK.A_Hbar * 1j*( dataK.E_K[:,None,:,None] - dataK.E_K[:,:,None,None] ) ) / HARTREETOEV / BOHRTOANG
#velocities_k = V_H_diag / HARTREETOEV / BOHRTOANG
# read in hoppings and proj_mat
hopping[:,0,range(hopping.shape[2]),range(hopping.shape[3])] = dataK.E_K
for isp in range(n_inequiv_spin_blocks):
iorb = 0
for icrsh in range(self.n_corr_shells):
dim = self.corr_shells[icrsh]['dim']
proj_mat[:,isp,icrsh,0:dim,:] = dataK.UU_K[:,iorb:iorb+dim,:]
iorb += dim
# read in rest from dataK
cell_volume = dataK.cell_volume / BOHRTOANG ** 3
kpts = dataK.kpoints_all
# broadcast everything
velocities_k = mpi.bcast(velocities_k)
cell_volume = mpi.bcast(cell_volume)
kpts = mpi.bcast(kpts)
hopping = mpi.bcast(hopping)
proj_mat = mpi.bcast(proj_mat)
# upgrade sumk quantities for interpolation
things_to_modify['kpts'] = kpts
things_to_modify['hopping'] = hopping
things_to_modify['proj_mat'] = proj_mat
mpi.barrier()
if mpi.is_master_node():
print(self.n_k, nk_x, nk_y, nk_z)
for key in things_to_modify:
things_to_store[key] = getattr(self, key)
setattr(self, key, things_to_modify[key])
#if mpi.is_master_node():
#print(key, things_to_store[key] )
#print(getattr(self, key))
# write velocities to file
if mpi.is_master_node():
ar = HDFArchive(self.hdf_file, 'a')
ar['dft_transp_input']['velocities_k'] = velocities_k
if mpi.is_master_node():
# k-dependent-projections.
# to be checked. But this should be obsolete atm, works for both cases
# k_dep_projection is nowhere used
# assert sum_k.k_dep_projection == 0, "transport_distribution: k dependent projection is not implemented!"
# positive Om_mesh
assert all(
Om >= 0.0 for Om in Om_mesh), "transport_distribution: Om_mesh should not contain negative values!"
# Check if energy_window is sufficiently large and correct
if (energy_window[0] >= energy_window[1] or energy_window[0] >= 0 or energy_window[1] <= 0):
assert 0, "transport_distribution: energy_window wrong!"
if (abs(self.fermi_dis(energy_window[0], beta) * self.fermi_dis(-energy_window[0], beta)) > 1e-5
or abs(self.fermi_dis(energy_window[1], beta) * self.fermi_dis(-energy_window[1], beta)) > 1e-5):
mpi.report(
"\n####################################################################")
mpi.report(
"transport_distribution: WARNING - energy window might be too narrow!")
mpi.report(
"####################################################################\n")
# up and down are equivalent if SP = 0
self.directions = directions
dir_to_int = {'x': 0, 'y': 1, 'z': 2}
# calculate A(k,w)
#######################################
# Define mesh for Green's function and in the specified energy window
if (with_Sigma == True):
self.omega = numpy.array([round(x.real, 12)
for x in self.mesh])
mesh = None
mu = self.chemical_potential
n_om = len(self.omega)
mpi.report("Using omega mesh provided by Sigma!")
if energy_window:
# Find according window in Sigma mesh
ioffset = numpy.sum(
self.omega < energy_window[0] - max(Om_mesh))
self.omega = self.omega[numpy.logical_and(self.omega >= energy_window[
0] - max(Om_mesh), self.omega <= energy_window[1] + max(Om_mesh))]
n_om = len(self.omega)
# 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 now we stick with this:
for icrsh in range(self.n_corr_shells):
Sigma_save = self.Sigma_imp[icrsh].copy()
spn = self.spin_block_names[self.corr_shells[icrsh]['SO']]
glist = lambda: [GfReFreq(target_shape=(block_dim, block_dim), window=(self.omega[
0], self.omega[-1]), n_points=n_om) for block, block_dim in self.gf_struct_sumk[icrsh]]
self.Sigma_imp[icrsh] = BlockGf(
name_list=spn, block_list=glist(), make_copies=False)
for i, g in self.Sigma_imp[icrsh]:
for iL in g.indices[0]:
for iR in g.indices[0]:
for iom in range(n_om):
g.data[iom, int(iL), int(iR)] = Sigma_save[
i].data[ioffset + iom, int(iL), int(iR)]
else:
assert n_om is not None, "transport_distribution: Number of omega points (n_om) needed to calculate transport distribution!"
assert energy_window is not None, "transport_distribution: Energy window needed to calculate transport distribution!"
assert broadening != 0.0 and broadening is not None, "transport_distribution: Broadening necessary to calculate transport distribution!"
self.omega = numpy.linspace(
energy_window[0] - max(Om_mesh), energy_window[1] + max(Om_mesh), n_om)
mesh = MeshReFreq(energy_window[0] -
max(Om_mesh), energy_window[1] + max(Om_mesh), n_om)
mu = 0.0
# Define mesh for optic conductivity
d_omega = round(numpy.abs(self.omega[0] - self.omega[1]), 12)
iOm_mesh = numpy.array([round((Om / d_omega), 0) for Om in Om_mesh])
self.Om_mesh = iOm_mesh * d_omega
if mpi.is_master_node():
print("Chemical potential: ", mu)
print("Using n_om = %s points in the energy_window [%s,%s]" % (n_om, self.omega[0], self.omega[-1]), end=' ')
print("where the omega vector is:")
print(self.omega)
print("Calculation requested for Omega mesh: ", numpy.array(Om_mesh))
print("Omega mesh automatically repined to: ", self.Om_mesh)
self.Gamma_w = {direction: numpy.zeros(
(len(self.Om_mesh), n_om), dtype=float) for direction in self.directions}
max_orb = numpy.max([self.n_orbitals[ik][0] for ik in range(self.n_k)])
#Akw_write = numpy.zeros((self.n_k, max_orb, max_orb, n_om), dtype=numpy.complex_)
# Sum over all k-points
ikarray = numpy.array(list(range(self.n_k)))
for ik in mpi.slice_array(ikarray):
# Calculate G_w for ik and initialize A_kw
G_w = self.lattice_gf(ik, mu, broadening=broadening, mesh=mesh, with_Sigma=with_Sigma)
A_kw = [numpy.zeros((self.n_orbitals[ik][isp], self.n_orbitals[ik][isp], n_om), dtype=complex)
for isp in range(n_inequiv_spin_blocks)]
for isp in range(n_inequiv_spin_blocks):
# copy data from G_w (swapaxes is used to have omega in the 3rd
# dimension)
A_kw[isp] = copy.deepcopy(G_w[self.spin_block_names[self.SO][
isp]].data.swapaxes(0, 1).swapaxes(1, 2))
# calculate A(k,w) for each frequency
for iw in range(n_om):
A_kw[isp][:, :, iw] = -1.0 / (2.0 * numpy.pi * 1j) * (
A_kw[isp][:, :, iw] - numpy.conjugate(numpy.transpose(A_kw[isp][:, :, iw])))
#Akw_write[ik] = A_kw[isp].copy() * self.bz_weights[ik]
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)
# loop over all symmetries
for R in self.rot_symmetries:
# get transformed velocity under symmetry R
if code in ('wien2k'):
vel_R = copy.deepcopy(self.velocities_k[isp][ik])
elif code in ('wannier90'):
vel_R = copy.deepcopy(velocities_k[ik])
for nu1 in range(self.band_window_optics[isp][ik, 1] - self.band_window_optics[isp][ik, 0] + 1):
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][:])
# calculate Gamma_w for each direction from the velocities
# vel_R and the spectral function A_kw
for direction in self.directions:
for iw in range(n_om):
for iq in range(len(self.Om_mesh)):
if(iw + iOm_mesh[iq] >= n_om or self.omega[iw] < -self.Om_mesh[iq] + energy_window[0] or self.omega[iw] > self.Om_mesh[iq] + energy_window[1]):
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, int(iw + iOm_mesh[iq])]), vel_R[v_i, v_i, dir_to_int[direction[1]]]),
A_kw[isp][A_i, A_i, iw]).trace().real * self.bz_weights[ik])
#Akw_write = mpi.all_reduce(mpi.world, Akw_write, lambda x, y: x + y)
#mpi.barrier()
#if mpi.is_master_node():
# ar = HDFArchive(self.hdf_file, 'a')
# ar.create_group('Akw')
# ar['Akw'] = numpy.sum(numpy.trace(Akw_write, axis1=1, axis2=2), axis=0)
for direction in self.directions:
self.Gamma_w[direction] = (mpi.all_reduce(mpi.world, self.Gamma_w[direction], lambda x, y: x + y) / self.cell_vol / self.n_symmetries)
def transport_coefficient(self, direction, iq, n, beta, method=None):
r"""
Calculates the transport coefficient A_n in a given direction for a given :math:`\Omega`. The required members (Gamma_w, directions, Om_mesh) have to be obtained first
by calling the function :meth:`transport_distribution <dft.sumk_dft_tools.SumkDFTTools.transport_distribution>`. For n>0 A is set to NaN if :math:`\Omega` is not 0.0.
Parameters
----------
direction : string
:math:`\alpha\beta` e.g.: 'xx','yy','zz','xy','xz','yz'.
iq : integer
Index of :math:`\Omega` point in the member Om_mesh.
n : integer
Number of the desired moment of the transport distribution.
beta : double
Inverse temperature :math:`\beta`.
method : string
Integration method: cubic spline and scipy.integrate.quad ('quad'), simpson rule ('simps'), trapezoidal rule ('trapz'), rectangular integration (otherwise)
Note that the sampling points of the the self-energy are used!
Returns
-------
A : double
Transport coefficient.
"""
if not (mpi.is_master_node()):
return
assert hasattr(
self, 'Gamma_w'), "transport_coefficient: Run transport_distribution first or load data from h5!"
if (self.Om_mesh[iq] == 0.0 or n == 0.0):
A = 0.0
# setup the integrand
if (self.Om_mesh[iq] == 0.0):
A_int = self.Gamma_w[direction][iq] * (self.fermi_dis(
self.omega, beta) * self.fermi_dis(-self.omega, beta)) * (self.omega * beta)**n
elif (n == 0.0):
A_int = self.Gamma_w[direction][iq] * (self.fermi_dis(self.omega, beta) - self.fermi_dis(
self.omega + self.Om_mesh[iq], beta)) / (self.Om_mesh[iq] * beta)
# w-integration
if method == 'quad':
# quad on interpolated w-points with cubic spline
A_int_interp = interp1d(self.omega, A_int, kind='cubic')
A = quad(A_int_interp, min(self.omega), max(self.omega),
epsabs=1.0e-12, epsrel=1.0e-12, limit=500)
A = A[0]
elif method == 'simps':
# simpson rule for w-grid
A = simps(A_int, self.omega)
elif method == 'trapz':
# trapezoidal rule for w-grid
A = numpy.trapz(A_int, self.omega)
else:
# rectangular integration for w-grid (orignal implementation)
d_w = self.omega[1] - self.omega[0]
for iw in range(self.Gamma_w[direction].shape[1]):
A += A_int[iw] * d_w
A = A * numpy.pi * (2.0 - self.SP)
else:
A = numpy.nan
return A
def conductivity_and_seebeck(self, beta, method=None):
r"""
Calculates the Seebeck coefficient and the optical conductivity by calling
:meth:`transport_coefficient <dft.sumk_dft_tools.SumkDFTTools.transport_coefficient>`.
The required members (Gamma_w, directions, Om_mesh) have to be obtained first by calling the function
:meth:`transport_distribution <dft.sumk_dft_tools.SumkDFTTools.transport_distribution>`.
Parameters
----------
beta : double
Inverse temperature :math:`\beta`.
Returns
-------
optic_cond : dictionary of double vectors
Optical conductivity in each direction and frequency given by Om_mesh.
seebeck : dictionary of double
Seebeck coefficient in each direction. If zero is not present in Om_mesh the Seebeck coefficient is set to NaN.
kappa : dictionary of double.
thermal conductivity in each direction. If zero is not present in Om_mesh the thermal conductivity is set to NaN
"""
if not (mpi.is_master_node()):
return
assert hasattr(
self, 'Gamma_w'), "conductivity_and_seebeck: Run transport_distribution first or load data from h5!"
n_q = self.Gamma_w[self.directions[0]].shape[0]
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}
A2 = {direction: numpy.full((n_q,), numpy.nan)
for direction in self.directions}
self.seebeck = {direction: numpy.nan for direction in self.directions}
self.kappa = {direction: numpy.nan for direction in self.directions}
self.optic_cond = {direction: numpy.full(
(n_q,), numpy.nan) for direction in self.directions}
for direction in self.directions:
for iq in range(n_q):
A0[direction][iq] = self.transport_coefficient(
direction, iq=iq, n=0, beta=beta, method=method)
A1[direction][iq] = self.transport_coefficient(
direction, iq=iq, n=1, beta=beta, method=method)
A2[direction][iq] = self.transport_coefficient(
direction, iq=iq, n=2, beta=beta, method=method)
print("A_0 in direction %s for Omega = %.2f %e a.u." % (direction, self.Om_mesh[iq], A0[direction][iq]))
print("A_1 in direction %s for Omega = %.2f %e a.u." % (direction, self.Om_mesh[iq], A1[direction][iq]))
print("A_2 in direction %s for Omega = %.2f %e a.u." % (direction, self.Om_mesh[iq], A2[direction][iq]))
if ~numpy.isnan(A1[direction][iq]):
# Seebeck and kappa are overwritten if there is more than one Omega =
# 0 in Om_mesh
self.seebeck[direction] = - \
A1[direction][iq] / A0[direction][iq] * 86.17
self.kappa[direction] = A2[direction][iq] - A1[direction][iq]*A1[direction][iq]/A0[direction][iq]
self.kappa[direction] *= 293178.0
self.optic_cond[direction] = beta * \
A0[direction] * 10700.0 / numpy.pi
for iq in range(n_q):
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]))
if not (numpy.isnan(A1[direction][iq])):
print("Seebeck in direction %s for Omega = 0.00 %f x 10^(-6) V/K" % (direction, self.seebeck[direction]))
print("kappa in direction %s for Omega = 0.00 %f W/(m * K)" % (direction, self.kappa[direction]))
return self.optic_cond, self.seebeck, self.kappa
def fermi_dis(self, w, beta):
r"""
Fermi distribution.
.. math::
f(x) = 1/(e^x+1).
Parameters
----------
w : double
frequency
beta : double
inverse temperature
Returns
-------
f : double
"""
return 1.0 / (numpy.exp(w * beta) + 1)

30
python/triqs_dft_tools/sumk_dft_transport.py
View File

@ -47,12 +47,19 @@ def read_transport_input_from_hdf(sum_k):
sum_k : sum_k object
triqs SumkDFT object
"""
assert sum_k.dft_code in ('wien2k','elk'), "read_transport_input_from_hdf() is only implemented for wien2k and elk inputs"
thingstoread = ['band_window_optics', 'velocities_k']
sum_k.read_input_from_hdf(
subgrp=sum_k.transp_data, things_to_read=thingstoread)
if(sum_k.dft_code=="wien2k"):
thingstoread = ['band_window', 'lattice_angles', 'lattice_constants',
'lattice_type', 'n_symmetries', 'rot_symmetries']
elif(sum_k.dft_code=="elk"):
thingstoread = ['band_window', 'n_symmetries',
'rot_symmetries','cell_vol']
sum_k.read_input_from_hdf(subgrp=sum_k.misc_data, things_to_read=thingstoread)
if(self.dft_code=="wien2k"):
self.cell_vol = self.cellvolume(self.lattice_type, self.lattice_constants, self.lattice_angles)[1]
return sum_k
@ -278,7 +285,7 @@ def recompute_w90_input_on_different_mesh(sum_k, seedname, nk_optics, pathname='
kpts = dataK.kpoints_all
# broadcast everything
cell_volume = mpi.bcast(cell_volume)
sum_k.cell_vol = mpi.bcast(cell_volume)
kpts = mpi.bcast(kpts)
hopping = mpi.bcast(hopping)
proj_mat = mpi.bcast(proj_mat)
@ -319,7 +326,7 @@ def recompute_w90_input_on_different_mesh(sum_k, seedname, nk_optics, pathname='
ar = HDFArchive(sum_k.hdf_file, 'a')
ar['dft_transp_input']['inverse_mass'] = inverse_mass
return sum_k, cell_volume, things_to_store
return sum_k, things_to_store
# ----------------- transport -----------------------
@ -341,15 +348,13 @@ def init_spectroscopy(sum_k, code='wien2k', w90_params={}):
-------
sum_k : sum_k object
triqs SumkDFT object, interpolated
cell_volume : double
primitive unit cell volume
"""
n_inequiv_spin_blocks = sum_k.SP + 1 - sum_k.SO
# up and down are equivalent if SP = 0
# ----------------- set-up input from DFT -----------------------
if code in ('wien2k'):
if code in ('wien2k', 'elk'):
# Check if wien converter was called and read transport subgroup form
# hdf file
if mpi.is_master_node():
@ -361,7 +366,6 @@ def init_spectroscopy(sum_k, code='wien2k', w90_params={}):
raise IOError("transport_distribution: n_symmetries missing. Check if case.outputs file is present and call convert_misc_input() or convert_dft_input().")
sum_k = read_transport_input_from_hdf(sum_k)
_, cell_volume = cellvolume(sum_k.lattice_type, sum_k.lattice_constants, sum_k.lattice_angles)
elif code in ('wannier90'):
required_entries = ['seedname', 'nk_optics']
@ -380,7 +384,7 @@ def init_spectroscopy(sum_k, code='wien2k', w90_params={}):
assert all(isinstance(name, bool) for name in [calc_velocity, calc_inverse_mass]), f'Parameter {calc_velocity} or {calc_inverse_mass} not bool!'
# recompute sum_k instances on denser grid
sum_k, cell_volume, _ = recompute_w90_input_on_different_mesh(sum_k, w90_params['seedname'], nk_optics=w90_params['nk_optics'], pathname=pathname,
sum_k, _ = recompute_w90_input_on_different_mesh(sum_k, w90_params['seedname'], nk_optics=w90_params['nk_optics'], pathname=pathname,
calc_velocity=calc_velocity, calc_inverse_mass=calc_inverse_mass)
# k-dependent-projections.
@ -388,10 +392,10 @@ def init_spectroscopy(sum_k, code='wien2k', w90_params={}):
# k_dep_projection is nowhere used
# assert sum_k.k_dep_projection == 0, "transport_distribution: k dependent projection is not implemented!"
return sum_k, cell_volume
return sum_k
# Uses .data of only GfReFreq objects.
def transport_distribution(sum_k, beta, cell_volume, directions=['xx'], energy_window=None, Om_mesh=[0.0], with_Sigma=False, n_om=None, broadening=0.0, code='wien2k'):
def transport_distribution(sum_k, beta, directions=['xx'], energy_window=None, Om_mesh=[0.0], with_Sigma=False, n_om=None, broadening=0.0, code='wien2k'):
r"""
Calculates the transport distribution
@ -406,8 +410,6 @@ def transport_distribution(sum_k, beta, cell_volume, directions=['xx'], energy_w
triqs SumkDFT object
beta : double
Inverse temperature :math:`\beta`.
cell_volume : double
primitive unit cell volume
directions : list of string, optional
:math:`\alpha\beta` e.g.: ['xx','yy','zz','xy','xz','yz'].
energy_window : list of double, optional
@ -495,8 +497,8 @@ def transport_distribution(sum_k, beta, cell_volume, directions=['xx'], energy_w
assert broadening != 0.0 and broadening is not None, "transport_distribution: Broadening necessary to calculate transport distribution!"
omega = numpy.linspace(
energy_window[0] - max(Om_mesh), energy_window[1] + max(Om_mesh), n_om)
mesh = [energy_window[0] -
max(Om_mesh), energy_window[1] + max(Om_mesh), n_om]
mesh = MeshReFreq(energy_window[0] -
max(Om_mesh), energy_window[1] + max(Om_mesh), n_om)
mu = 0.0
dir_to_int = {'x': 0, 'y': 1, 'z': 2}
@ -569,7 +571,7 @@ def transport_distribution(sum_k, beta, cell_volume, directions=['xx'], energy_w
for direction in directions:
Gamma_w[direction] = (mpi.all_reduce(mpi.world, Gamma_w[direction], lambda x, y: x + y) / cell_volume / sum_k.n_symmetries)
Gamma_w[direction] = (mpi.all_reduce(mpi.world, Gamma_w[direction], lambda x, y: x + y) / sum_k.cell_vol / sum_k.n_symmetries)
return Gamma_w, omega, temp_Om_mesh