Source code for wannierberri.calculators.dynamic
import numpy as np
import abc
import functools
from ..__utility import Gaussian, Lorentzian, FermiDirac, alpha_A, beta_A
from ..result import EnergyResult
from .calculator import Calculator
from ..formula.covariant import SpinVelocity
from ..formula import Formula
from copy import copy
from ..point_symmetry import transform_ident, transform_trans, transform_odd, transform_odd_trans_021
from scipy.constants import elementary_charge, hbar, electron_mass, physical_constants, angstrom
from .. import __factors as factors
bohr_magneton = elementary_charge * hbar / (2 * electron_mass)
bohr = physical_constants['Bohr radius'][0] / angstrom
eV_au = physical_constants['electron volt-hartree relationship'][0]
Ang_SI = angstrom
#######################################
# #
# integration (Efermi-omega) #
# #
#######################################
[docs]
class DynamicCalculator(Calculator, abc.ABC):
def __init__(self, Efermi=None, omega=None, kBT=0, smr_fixed_width=0.1, smr_type='Lorentzian',
kwargs_formula=None, dtype=complex, **kwargs):
if kwargs_formula is None:
kwargs_formula = {}
super().__init__(**kwargs)
self.Efermi = Efermi
self.omega = omega
self.kBT = kBT
self.smr_fixed_width = smr_fixed_width
self.smr_type = smr_type
self.kwargs_formula = copy(kwargs_formula)
self.Formula = None
self.constant_factor = 1.
self.dtype = dtype
self.EFmin = self.Efermi.min()
self.EFmax = self.Efermi.max()
self.omegamin = self.omega.min()
self.omegamax = self.omega.max()
self.eocc1max = self.EFmin - 30 * self.kBT
self.eocc0min = self.EFmax + 30 * self.kBT
if self.smr_type == 'Lorentzian':
self.smear = functools.partial(Lorentzian, width=self.smr_fixed_width)
elif self.smr_type == 'Gaussian':
self.smear = functools.partial(Gaussian, width=self.smr_fixed_width, adpt_smr=False)
else:
raise ValueError("Invalid smearing type {self.smr_type}")
self.FermiDirac = functools.partial(FermiDirac, mu=self.Efermi, kBT=self.kBT)
@abc.abstractmethod
def factor_omega(self, E1, E2):
"""determines a frequency-dependent factor for bands with energies E1 and E2"""
def factor_Efermi(self, E1, E2):
return self.FermiDirac(E2) - self.FermiDirac(E1)
def nonzero(self, E1, E2):
if (E1 < self.eocc1max and E2 < self.eocc1max) or (E1 > self.eocc0min and E2 > self.eocc0min):
return False
else:
return True
def __call__(self, data_K):
formula = self.Formula(data_K, **self.kwargs_formula)
restot_shape = (len(self.omega), len(self.Efermi)) + (3,) * formula.ndim
restot_shape_tmp = (
len(self.omega), len(self.Efermi) * 3 ** formula.ndim) # we will first get it in this shape, then transpose
restot = np.zeros(restot_shape_tmp, self.dtype)
for ik in range(data_K.nk):
degen_groups = data_K.get_bands_in_range_groups_ik(
ik, -np.Inf, np.Inf, degen_thresh=self.degen_thresh, degen_Kramers=self.degen_Kramers)
# now find needed pairs:
# as a dictionary {((ibm1,ibm2),(ibn1,ibn2)):(Em,En)}
degen_group_pairs = [
(ibm, ibn, Em, En) for ibm, Em in degen_groups.items() for ibn, En in degen_groups.items()
if self.nonzero(Em, En)
]
npair = len(degen_group_pairs)
if npair == 0:
continue
matrix_elements = np.array(
[formula.trace_ln(ik, np.arange(*pair[0]), np.arange(*pair[1])) for pair in degen_group_pairs])
factor_Efermi = np.array([self.factor_Efermi(pair[2], pair[3]) for pair in degen_group_pairs])
factor_omega = np.array([self.factor_omega(pair[2], pair[3]) for pair in degen_group_pairs]).T
restot += factor_omega @ (factor_Efermi[:, :, None] *
matrix_elements.reshape(npair, -1)[:, None, :]).reshape(npair, -1)
restot = restot.reshape(restot_shape).swapaxes(0, 1) # swap the axes to get EF,omega,a,b,...
restot *= self.constant_factor / (data_K.nk * data_K.cell_volume)
try:
transformTR = self.transformTR
except AttributeError:
transformTR = formula.transformTR
try:
transformInv = self.transformInv
except AttributeError:
transformInv = formula.transformInv
return EnergyResult(
[self.Efermi, self.omega], restot,
transformTR=transformTR, transformInv=transformInv,
save_mode=self.save_mode
)
[docs]
class MultitermCalculator(DynamicCalculator):
def __init__(self, **kwargs):
super().__init__(**kwargs)
self.terms = []
def __call__(self, data_K):
return sum(cal(data_K) for cal in self.terms)
###############################################
###############################################
###############################################
###############################################
#### ######
#### Implemented calculators ######
#### ######
###############################################
###############################################
###############################################
###############################################
###############################
# JDOS #
###############################
[docs]
class Formula_dyn_ident(Formula):
def __init__(self, data_K):
super().__init__(data_K)
self.transformTR = transform_ident
self.transformInv = transform_ident
self.ndim = 0
[docs]
class JDOS(DynamicCalculator):
def __init__(self, **kwargs):
super().__init__(**kwargs)
self.sigma = self.smr_fixed_width
self.Formula = Formula_dyn_ident
self.dtype = float
[docs]
def nonzero(self, E1, E2):
return (E1 < self.Efermi.max()) and (E2 > self.Efermi.min()) and (
self.omega.min() - 5 * self.smr_fixed_width < E2 - E1 < self.omega.max() + 5 * self.smr_fixed_width)
[docs]
def energy_factor(self, E1, E2):
res = np.zeros((len(self.Efermi), len(self.omega)))
gauss = self.smear(E2 - E1 - self.omega, self.smr_fixed_width)
res[(E1 < self.Efermi) * (self.Efermi < E2)] = gauss[None, :]
return res
################################
# Optical Conductivity #
################################
[docs]
class Formula_OptCond(Formula):
def __init__(self, data_K, **parameters):
super().__init__(data_K, **parameters)
if self.external_terms:
A = data_K.A_H
else:
A = data_K.A_H_internal
self.AA = 1j * A[:, :, :, :, None] * A.swapaxes(1, 2)[:, :, :, None, :]
self.ndim = 2
self.transformTR = transform_trans
self.transformInv = transform_ident
[docs]
class OpticalConductivity(DynamicCalculator):
def __init__(self, **kwargs):
super().__init__(**kwargs)
self.Formula = Formula_OptCond
self.constant_factor = factors.factor_opt
[docs]
def factor_omega(self, E1, E2):
delta_arg_12 = E2 - E1 - self.omega # argument of delta function [iw, n, m]
cfac = 1. / (delta_arg_12 - 1j * self.smr_fixed_width)
if self.smr_type != 'Lorentzian':
cfac.imag = np.pi * self.smear(delta_arg_12)
return (E2 - E1) * cfac
###############################
# SHC #
###############################
[docs]
class Formula_SHC(Formula):
def __init__(self, data_K, SHC_type='ryoo', shc_abc=None, **parameters):
super().__init__(data_K, **parameters)
A = SpinVelocity(data_K, SHC_type, external_terms=self.external_terms).matrix
if self.external_terms:
B = -1j * data_K.A_H
else:
B = -1j * data_K.A_H_internal
self.imAB = np.imag(A[:, :, :, :, None, :] * B.swapaxes(1, 2)[:, :, :, None, :, None])
self.ndim = 3
if shc_abc is not None:
assert len(shc_abc) == 3
a, b, c = (x - 1 for x in shc_abc)
self.imAB = self.imAB[:, :, :, a, b, c]
self.ndim = 0
self.transformTR = transform_ident
self.transformInv = transform_ident
[docs]
class SHC(DynamicCalculator):
def __init__(self, SHC_type="ryoo", shc_abc=None, **kwargs):
super().__init__(**kwargs)
self.kwargs_formula.update(dict(SHC_type=SHC_type, shc_abc=shc_abc))
self.Formula = Formula_SHC
self.constant_factor = factors.factor_shc
[docs]
def factor_omega(self, E1, E2):
delta_arg_12 = E1 - E2 - self.omega # argument of delta function [iw, n, m]
cfac = 1. / (delta_arg_12 - 1j * self.smr_fixed_width)
if self.smr_type != 'Lorentzian':
cfac.imag = np.pi * self.smear(delta_arg_12)
return cfac / 2
####################################################
# Spatially-dispersive conductivity tensor #
####################################################
# To keep e^2/hbar as the global factor of SDCT
electron_g_factor = physical_constants['electron g factor'][0]
m_spin_prefactor = electron_g_factor * hbar / electron_mass
# _____ Antisymmetric (time-even) spatially-dispersive conductivity tensor _____ #
[docs]
class SDCT_asym(MultitermCalculator):
def __init__(self, fermi_sea=True, fermi_surf=True, **kwargs):
super().__init__(**kwargs)
# Fermi sea terms
if fermi_sea:
self.terms.extend([SDCT_asym_sea_I(**kwargs), SDCT_asym_sea_II(**kwargs)])
# Fermi surface terms
if fermi_surf:
self.terms.extend([SDCT_asym_surf_I(**kwargs), SDCT_asym_surf_II(**kwargs)])
[docs]
class Formula_SDCT_asym_sea_I(Formula):
def __init__(self, data_K, M1_terms=True, E2_terms=True, V_terms=True, spin=False, **parameters):
super().__init__(data_K, **parameters)
# Intrinsic multipole moments
if self.external_terms:
A = -1. * data_K.E1
B_M1, B_E2, B = data_K.Bln
if spin:
S = data_K.Xbar('SS')
B_M1[:, :, :, alpha_A, beta_A] += -0.5 * m_spin_prefactor * S
B_M1[:, :, :, beta_A, alpha_A] -= -0.5 * m_spin_prefactor * S
else:
A = -1. * data_K.E1_internal
B_M1, B_E2, B = data_K.Bln_internal
# Other quantities
Vn = data_K.Vn
Vnm_plus = 0.5 * (Vn[:, :, None, :] + Vn[:, None, :, :])
# --- Formula --- #
summ = np.zeros((data_K.nk, data_K.num_wann, data_K.num_wann, 3, 3, 3), dtype=complex)
if M1_terms:
summ += -np.imag(A[:, :, :, :, None, None] * B_M1.swapaxes(1, 2)[:, :, :, None, :, :])
if E2_terms:
summ += -np.imag(A[:, :, :, :, None, None] * B_E2.swapaxes(1, 2)[:, :, :, None, :, :])
if V_terms:
summ += Vnm_plus[:, :, :, :, None, None] * np.imag(A[:, :, :, None, :, None] * A.swapaxes(1, 2)[:, :, :, None, None, :])
summ = summ - summ.swapaxes(3, 4)
self.summ = summ
self.ndim = 3
self.transformTR = transform_ident
self.transformInv = transform_odd
[docs]
class SDCT_asym_sea_I(DynamicCalculator):
def __init__(self, M1_terms=True, E2_terms=True, V_terms=True, spin=False, **kwargs):
super().__init__(**kwargs)
self.kwargs_formula.update(dict(M1_terms=M1_terms, E2_terms=E2_terms, V_terms=V_terms, spin=spin))
self.Formula = Formula_SDCT_asym_sea_I
self.constant_factor = factors.factor_SDCT
[docs]
def factor_omega(self, E1, E2):
omega = self.omega + 1.j * self.smr_fixed_width
Z_arg_12 = (E2 - E1)**2 - omega**2 # argument of Z_ln function [iw, n, m]
Zfac = 1. / Z_arg_12
return omega * Zfac
[docs]
class Formula_SDCT_asym_sea_II(Formula):
def __init__(self, data_K, M1_terms=True, E2_terms=True, V_terms=True, spin=False, **parameters):
super().__init__(data_K, **parameters)
# Intrinsic multipole moments
if self.external_terms:
A = -1. * data_K.E1
else:
A = -1. * data_K.E1_internal
# Other quantities
Vn = data_K.Vn
Vnm_plus = 0.5 * (Vn[:, :, None, :] + Vn[:, None, :, :])
# --- Formula --- #
summ = np.zeros((data_K.nk, data_K.num_wann, data_K.num_wann, 3, 3, 3), dtype=complex)
if V_terms:
summ += np.imag(A[:, :, :, :, None, None] * A.swapaxes(1, 2)[:, :, :, None, :, None]) * Vnm_plus[:, :, :, None, None, :]
self.summ = summ
self.ndim = 3
self.transformTR = transform_ident
self.transformInv = transform_odd
[docs]
class SDCT_asym_sea_II(DynamicCalculator):
def __init__(self, M1_terms=True, E2_terms=True, V_terms=True, spin=False, **kwargs):
super().__init__(**kwargs)
self.kwargs_formula.update(dict(M1_terms=M1_terms, E2_terms=E2_terms, V_terms=V_terms, spin=spin))
self.Formula = Formula_SDCT_asym_sea_II
self.constant_factor = factors.factor_SDCT
[docs]
def factor_omega(self, E1, E2):
omega = self.omega + 1.j * self.smr_fixed_width
Z_arg_12 = (E2 - E1)**2 - omega**2 # argument of Z_ln function [iw, n, m]
Zfac = 1. / Z_arg_12
return omega * (3.0 * (E2 - E1)**2 - omega**2) * Zfac**2
[docs]
class Formula_SDCT_asym_surf_I(Formula):
def __init__(self, data_K, M1_terms=True, E2_terms=True, V_terms=True, spin=False, **parameters):
super().__init__(data_K, **parameters)
# Intrinsic multipole moments
if self.external_terms:
A = -1. * data_K.E1
else:
A = -1. * data_K.E1_internal
# Other quantities
Vn = data_K.Vn
# --- Formula --- #
summ = np.zeros((data_K.nk, data_K.num_wann, data_K.num_wann, 3, 3, 3), dtype=complex)
if V_terms:
summ += -np.imag(A[:, :, :, :, None, None] * A.swapaxes(1, 2)[:, :, :, None, :, None]) * Vn[:, :, None, None, None, :]
self.summ = summ
self.ndim = 3
self.transformTR = transform_ident
self.transformInv = transform_odd
[docs]
class SDCT_asym_surf_I(DynamicCalculator):
def __init__(self, M1_terms=True, E2_terms=True, V_terms=True, spin=False, **kwargs):
super().__init__(**kwargs)
self.kwargs_formula.update(dict(M1_terms=M1_terms, E2_terms=E2_terms, V_terms=V_terms, spin=spin))
self.Formula = Formula_SDCT_asym_surf_I
self.constant_factor = factors.factor_SDCT
[docs]
def factor_omega(self, E1, E2):
omega = self.omega + 1.j * self.smr_fixed_width
Z_arg_12 = (E2 - E1)**2 - omega**2 # argument of Z_ln function [iw, n, m]
Zfac = 1. / Z_arg_12
return omega * (E2 - E1) * Zfac
[docs]
def factor_Efermi(self, E1, E2):
return -self.FermiDirac(E1)**2 * np.exp((E1 - self.Efermi) / self.kBT) / self.kBT
[docs]
class Formula_SDCT_asym_surf_II(Formula):
def __init__(self, data_K, M1_terms=True, E2_terms=True, V_terms=True, spin=False, **parameters):
super().__init__(data_K, **parameters)
# Intrinsic multipole moments
if self.external_terms:
B_M1, B_E2, B = data_K.Bln
if spin:
S = data_K.Xbar('SS')
B_M1[:, :, :, alpha_A, beta_A] += -0.5 * m_spin_prefactor * S
B_M1[:, :, :, beta_A, alpha_A] -= -0.5 * m_spin_prefactor * S
else:
B_M1, B_E2, B = data_K.Bln_internal
Bn_M1 = np.diagonal(B_M1, axis1=1, axis2=2).transpose(0, 3, 1, 2)
# Other quantities
Vn = data_K.Vn
# --- Formula --- #
summ = np.zeros((data_K.nk, data_K.num_wann, 3, 3, 3), dtype=complex)
if M1_terms:
summ += Vn[:, :, :, None, None] * Bn_M1[:, :, None, :, :]
summ = summ - summ.swapaxes(3, 4)
self.summ = summ
self.ndim = 3
self.transformTR = transform_ident
self.transformInv = transform_odd
[docs]
def trace_ln(self, ik, inn1, inn2): # There is no sum over l
return self.summ[ik, inn1].sum(axis=0)
[docs]
class SDCT_asym_surf_II(DynamicCalculator):
def __init__(self, M1_terms=True, E2_terms=True, V_terms=True, spin=False, **kwargs):
super().__init__(**kwargs)
self.kwargs_formula.update(dict(M1_terms=M1_terms, E2_terms=E2_terms, V_terms=V_terms, spin=spin))
self.Formula = Formula_SDCT_asym_surf_II
self.constant_factor = factors.factor_SDCT
[docs]
def factor_omega(self, E1, E2):
omega = self.omega + 1.j * self.smr_fixed_width
return 1. / omega
[docs]
def factor_Efermi(self, E1, E2):
return -self.FermiDirac(E1)**2 * np.exp((E1 - self.Efermi) / self.kBT) / self.kBT
# _____ Symmetric (time-odd) spatially-dispersive conductivity tensor _____ #
[docs]
class SDCT_sym(MultitermCalculator):
def __init__(self, fermi_sea=True, fermi_surf=True, **kwargs):
super().__init__(**kwargs)
# Fermi sea terms
if fermi_sea:
self.terms.extend([SDCT_sym_sea_I(**kwargs), SDCT_sym_sea_II(**kwargs)])
# Fermi surface terms
if fermi_surf:
self.terms.extend([SDCT_sym_surf_I(**kwargs), SDCT_sym_surf_II(**kwargs)])
[docs]
class Formula_SDCT_sym_sea_I(Formula):
def __init__(self, data_K, M1_terms=True, E2_terms=True, V_terms=True, spin=False, **parameters):
super().__init__(data_K, **parameters)
# Intrinsic multipole moments
if self.external_terms:
A = -1. * data_K.E1
B_M1, B_E2, B = data_K.Bln
if spin:
S = data_K.Xbar('SS')
B_M1[:, :, :, alpha_A, beta_A] += -0.5 * m_spin_prefactor * S
B_M1[:, :, :, beta_A, alpha_A] -= -0.5 * m_spin_prefactor * S
else:
A = -1. * data_K.E1_internal
B_M1, B_E2, B = data_K.Bln_internal
# Other quantities
Vn = data_K.Vn
Vnm_plus = 0.5 * (Vn[:, :, None, :] + Vn[:, None, :, :])
# Formula
summ = np.zeros((data_K.nk, data_K.num_wann, data_K.num_wann, 3, 3, 3), dtype=complex)
if M1_terms:
summ += np.real(A[:, :, :, :, None, None] * B_M1.swapaxes(1, 2)[:, :, :, None, :, :])
if E2_terms:
summ += np.real(A[:, :, :, :, None, None] * B_E2.swapaxes(1, 2)[:, :, :, None, :, :])
if V_terms:
summ += Vnm_plus[:, :, :, :, None, None] * np.real(A[:, :, :, None, :, None] * A.swapaxes(1, 2)[:, :, :, None, None, :])
summ = summ + summ.swapaxes(3, 4)
self.summ = summ
self.ndim = 3
self.transformTR = transform_odd
self.transformInv = transform_odd
[docs]
class SDCT_sym_sea_I(DynamicCalculator):
def __init__(self, M1_terms=True, E2_terms=True, V_terms=True, spin=False, **kwargs):
super().__init__(**kwargs)
self.kwargs_formula.update(dict(M1_terms=M1_terms, E2_terms=E2_terms, V_terms=V_terms, spin=spin))
self.Formula = Formula_SDCT_sym_sea_I
self.constant_factor = factors.factor_SDCT
[docs]
def factor_omega(self, E1, E2):
omega = self.omega + 1.j * self.smr_fixed_width
Z_arg_12 = (E2 - E1)**2 - omega**2 # argument of Z_ln function [iw, n, m]
Zfac = 1. / Z_arg_12
return 1.j * (E2 - E1) * Zfac
[docs]
class Formula_SDCT_sym_sea_II(Formula):
def __init__(self, data_K, M1_terms=True, E2_terms=True, V_terms=True, spin=False, **parameters):
super().__init__(data_K, **parameters)
# Intrinsic multipole moments
if self.external_terms:
A = -1. * data_K.E1
else:
A = -1. * data_K.E1_internal
# Other quantities
Vn = data_K.Vn
Vnm_plus = Vn[:, :, None, :] + Vn[:, None, :, :]
# Formula
summ = np.zeros((data_K.nk, data_K.num_wann, data_K.num_wann, 3, 3, 3), dtype=complex)
if V_terms:
summ -= np.real(A[:, :, :, :, None, None] * A.swapaxes(1, 2)[:, :, :, None, :, None]) * Vnm_plus[:, :, :, None, None, :]
self.summ = summ
self.ndim = 3
self.transformTR = transform_odd
self.transformInv = transform_odd
[docs]
class SDCT_sym_sea_II(DynamicCalculator):
def __init__(self, M1_terms=True, E2_terms=True, V_terms=True, spin=False, **kwargs):
super().__init__(**kwargs)
self.kwargs_formula.update(dict(M1_terms=M1_terms, E2_terms=E2_terms, V_terms=V_terms, spin=spin))
self.Formula = Formula_SDCT_sym_sea_II
self.constant_factor = factors.factor_SDCT
[docs]
def factor_omega(self, E1, E2):
omega = self.omega + 1.j * self.smr_fixed_width
Z_arg_12 = (E2 - E1)**2 - omega**2 # argument of Z_ln function [iw, n, m]
Zfac = 1. / Z_arg_12
return 1.j * (E2 - E1)**3 * Zfac**2
[docs]
class Formula_SDCT_sym_surf_I(Formula):
def __init__(self, data_K, M1_terms=True, E2_terms=True, V_terms=True, spin=False, **parameters):
super().__init__(data_K, **parameters)
# Intrinsic multipole moments
if self.external_terms:
A = -1. * data_K.E1
else:
A = -1. * data_K.E1_internal
# Other quantities
Vn = data_K.Vn
# Formula
summ = np.zeros((data_K.nk, data_K.num_wann, data_K.num_wann, 3, 3, 3), dtype=complex)
if V_terms:
summ += np.real(A[:, :, :, :, None, None] * A.swapaxes(1, 2)[:, :, :, None, :, None]) * Vn[:, :, None, None, None, :]
self.summ = summ
self.ndim = 3
self.transformTR = transform_odd
self.transformInv = transform_odd
[docs]
class SDCT_sym_surf_I(DynamicCalculator):
def __init__(self, M1_terms=True, E2_terms=True, V_terms=True, spin=False, **kwargs):
super().__init__(**kwargs)
self.kwargs_formula.update(dict(M1_terms=M1_terms, E2_terms=E2_terms, V_terms=V_terms, spin=spin))
self.Formula = Formula_SDCT_sym_surf_I
self.constant_factor = factors.factor_SDCT
[docs]
def factor_omega(self, E1, E2):
omega = self.omega + 1.j * self.smr_fixed_width
Z_arg_12 = (E2 - E1)**2 - omega**2 # argument of Z_ln function [iw, n, m]
Zfac = 1. / Z_arg_12
return 1.j * (E2 - E1)**2 * Zfac
[docs]
def factor_Efermi(self, E1, E2):
return -self.FermiDirac(E1)**2 * np.exp((E1 - self.Efermi) / self.kBT) / self.kBT
[docs]
class Formula_SDCT_sym_surf_II(Formula):
def __init__(self, data_K, M1_terms=True, E2_terms=True, V_terms=True, spin=False, **parameters):
super().__init__(data_K, **parameters)
Vn = data_K.Vn
# Formula
summ = np.zeros((data_K.nk, data_K.num_wann, data_K.num_wann, 3, 3, 3), dtype=complex)
if V_terms:
summ = Vn[:, :, :, None, None] * Vn[:, :, None, :, None] * Vn[:, :, None, None, :]
self.summ = summ
self.ndim = 3
self.transformTR = transform_odd
self.transformInv = transform_odd
[docs]
def trace_ln(self, ik, inn1, inn2): # There is no sum over l
return self.summ[ik, inn1].sum(axis=0)
[docs]
class SDCT_sym_surf_II(DynamicCalculator):
def __init__(self, M1_terms=True, E2_terms=True, V_terms=True, spin=False, **kwargs):
super().__init__(**kwargs)
self.kwargs_formula.update(dict(M1_terms=M1_terms, E2_terms=E2_terms, V_terms=V_terms, spin=spin))
self.Formula = Formula_SDCT_sym_surf_II
self.constant_factor = factors.factor_SDCT
[docs]
def factor_omega(self, E1, E2):
omega = self.omega + 1.j * self.smr_fixed_width
return -1.j / omega**2
[docs]
def factor_Efermi(self, E1, E2):
return -self.FermiDirac(E1)**2 * np.exp((E1 - self.Efermi) / self.kBT) / self.kBT
###############################################################################
# ===============
# Shift current
# ===============
[docs]
class ShiftCurrentFormula(Formula):
def __init__(self, data_K, sc_eta, **parameters):
super().__init__(data_K, **parameters)
D_H = data_K.D_H
V_H = data_K.Xbar('Ham', 1)
del2E_H = data_K.Xbar('Ham', 2)
dEig_inv = data_K.dEig_inv.swapaxes(2, 1)
# define D using broadening parameter
E_K = data_K.E_K
dEig = E_K[:, :, None] - E_K[:, None, :]
dEig_inv_Pval = dEig / (dEig ** 2 + sc_eta ** 2)
D_H_Pval = -V_H * dEig_inv_Pval[:, :, :, None]
# commutators
# ** the spatial index of D_H_Pval corresponds to generalized derivative direction
# ** --> stored in the fourth column of output variables
sum_HD = (np.einsum('knlc,klma->knmca', V_H, D_H_Pval) -
np.einsum('knnc,knma->knmca', V_H, D_H_Pval) -
np.einsum('knla,klmc->knmca', D_H_Pval, V_H) +
np.einsum('knma,kmmc->knmca', D_H_Pval, V_H))
# ** this one is invariant under a<-->c
DV_bit = (np.einsum('knmc,knna->knmca', D_H, V_H) -
np.einsum('knmc,kmma->knmca', D_H, V_H) +
np.einsum('knma,knnc->knmca', D_H, V_H) -
np.einsum('knma,kmmc->knmca', D_H, V_H))
# generalized derivative
A_gen_der = (+ 1j * (del2E_H + sum_HD + DV_bit) * dEig_inv[:, :, :, np.newaxis, np.newaxis])
if self.external_terms:
# ** the spatial index of A_Hbar with diagonal terms corresponds to generalized derivative direction
# ** --> stored in the fourth column of output variables
A_Hbar = data_K.Xbar('AA')
A_Hbar_der = data_K.Xbar('AA', 1)
sum_AD = (np.einsum('knlc,klma->knmca', A_Hbar, D_H_Pval) -
np.einsum('knnc,knma->knmca', A_Hbar, D_H_Pval) -
np.einsum('knla,klmc->knmca', D_H_Pval, A_Hbar) +
np.einsum('knma,kmmc->knmca', D_H_Pval, A_Hbar))
AD_bit = (np.einsum('knnc,knma->knmac', A_Hbar, D_H) -
np.einsum('kmmc,knma->knmac', A_Hbar, D_H) +
np.einsum('knna,knmc->knmac', A_Hbar, D_H) -
np.einsum('kmma,knmc->knmac', A_Hbar, D_H))
AA_bit = (np.einsum('knnb,knma->knmab', A_Hbar, A_Hbar) -
np.einsum('kmmb,knma->knmab', A_Hbar, A_Hbar))
A_gen_der += A_Hbar_der + AD_bit - 1j * AA_bit + sum_AD
# generalized derivative is fourth index of A, we put it into third index of Imn
if self.external_terms:
A_H = data_K.A_H
else:
A_H = data_K.A_H_internal
# here we take the -real part to eliminate the 1j factor in the final factor
Imn = - np.einsum('knmca,kmnb->knmabc', A_gen_der, A_H).imag
Imn += Imn.swapaxes(4, 5) # symmetrize b and c
self.Imn = Imn
self.ndim = 3
self.transformTR = transform_ident
self.transformInv = transform_odd
[docs]
class ShiftCurrent(DynamicCalculator):
def __init__(self, sc_eta, **kwargs):
super().__init__(dtype=float, **kwargs)
self.kwargs_formula.update(dict(sc_eta=sc_eta))
self.Formula = ShiftCurrentFormula
self.constant_factor = factors.factor_shift_current
[docs]
def factor_omega(self, E1, E2):
delta_arg_12 = E1 - E2 - self.omega # argument of delta function [iw, n, m]
delta_arg_21 = E2 - E1 - self.omega
return self.smear(delta_arg_12) + self.smear(delta_arg_21)
# ===================
# Injection current
# ===================
[docs]
class InjectionCurrentFormula(Formula):
"""
Eq. (10) of Lihm and Park, PRB 105, 045201 (2022)
Use v_mn = i * r_mn * (e_m - e_n) / hbar to replace v with r.
"""
def __init__(self, data_K, **parameters):
super().__init__(data_K, **parameters)
if self.external_terms:
A_H = data_K.A_H
else:
A_H = data_K.A_H_internal
V_H = data_K.Xbar('Ham', 1) # (k, m, n, a)
V_H_diag = np.diagonal(V_H, axis1=1, axis2=2).transpose(0, 2, 1) # (k, m, a)
# compute delta_V[k, m, n, a] = V_H[k, m, m, a] - V_H[k, n, n, a]
delta_V = V_H_diag[:, :, None, :] - V_H_diag[:, None, :, :] # (k, m, n, a)
Imn = np.einsum('kmna,kmnb,knmc->kmnabc', delta_V, A_H, A_H)
self.Imn = Imn
self.ndim = 3
[docs]
class InjectionCurrent(DynamicCalculator):
def __init__(self, **kwargs):
super().__init__(**kwargs)
self.Formula = InjectionCurrentFormula
self.transformTR = transform_odd_trans_021
self.transformInv = transform_odd
self.constant_factor = factors.factor_injection_current
[docs]
def factor_omega(self, E1, E2):
delta_arg_12 = E1 - E2 - self.omega # argument of delta function [iw, n, m]
return self.smear(delta_arg_12)