import numpy as np
import abc
import functools
from ..utility import Gaussian, Lorentzian, FermiDirac
from ..result import EnergyResult
from .calculator import Calculator
from ..formula.covariant import SpinVelocity
from ..formula import Formula
from copy import copy
from ..symmetry.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
)
###############################################
###############################################
###############################################
###############################################
#### ######
#### Implemented calculators ######
#### ######
###############################################
###############################################
###############################################
###############################################
###############################
# JDOS #
###############################
[docs]
class JDOS(DynamicCalculator):
r"""Joint Density of States
:math:`\rho(\omega) = \sum_{\mathbf{k}} \sum_{m,n} \delta(E_{m\mathbf{k}} - E_{n\mathbf{k}} - \omega) \times \left(f(E_{n\mathbf{k}}) - f(E_{m\mathbf{k}})\right)`
"""
def __init__(self, **kwargs):
super().__init__(**kwargs)
self.sigma = self.smr_fixed_width
self.Formula = Formula_dyn_ident
self.dtype = float
[docs]
def factor_omega(self, E1, E2):
return self.smear(E2 - E1 - self.omega)
[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)
################################
# Optical Conductivity #
################################
[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 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
# ===============
# Shift current
# ===============
[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 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)