from abc import ABC, abstractmethod
from typing import Optional
import numpy as np
from ._common import invert, fermi_dirac, fermi_prime, primitive_fermi_dirac
__all__ = [
"EnergyKernel",
"KuboKernel",
"KuboBastinKernel",
"SecondOrderKuboKernel",
"FermiDirac",
"PrimitiveFermiDirac",
]
[docs]
class EnergyKernel(ABC):
def __init__(self):
self._contractions = "kij,kij->"
@property
def contractions(self):
return self._contractions
@contractions.setter
def contractions(self, value):
self._contractions = value
[docs]
@abstractmethod
def __call__(self, *args):
pass
[docs]
class FermiDirac(EnergyKernel):
def __init__(self):
super().__init__()
self.contractions = "ki,ki->"
[docs]
def __call__(
self,
energies: np.ndarray,
mu: float = 0.0,
temperature: float = 0.0,
eta: Optional[float] = None,
) -> np.ndarray:
return fermi_dirac(energies, mu, temperature)
[docs]
class PrimitiveFermiDirac(EnergyKernel):
def __init__(self):
super().__init__()
self.contractions = "ki,ki->"
[docs]
def __call__(
self,
energies: np.ndarray,
mu: float = 0.0,
temperature: float = 0.0,
eta: Optional[float] = None,
) -> np.ndarray:
return primitive_fermi_dirac(energies, mu, temperature)
class BerryKernel(EnergyKernel):
def __init__(self) -> None:
super().__init__()
self.contractions = 'kij,kij->kj'
def __call__(self, energies, mu=0, temperature=0, eta=0):
delta_e = energies[:, None, :] - energies[:, :, None]
return 2 * 1j * invert(delta_e * delta_e)
[docs]
class KuboKernel(EnergyKernel):
"""Interband, intraband, anomalous, even or odd kernels.
Attributes:
name: Specifies the response to be computed.
"""
def __init__(self, name: str = "intra_band") -> None:
super().__init__()
self.name = name
if name not in ["inter_band", "intra_band", "anomalous", "even", "odd"]:
raise ValueError(
"name should be either inter_band, intra_band, even, odd or anomalous"
)
if self.name == "intra_band":
self.contractions = "ki,kii->"
[docs]
def interband_kubo_kernel(
self,
energies: np.ndarray,
mu: float = 0.0,
temperature: float = 0.0,
eta: float = 0.0,
) -> np.ndarray:
r"""Computes the kubo kernel in Eq(A5) of [1]
.. math:
i\frac{f_m - f_n}{(E_n - E_m)(E_n - E_m + i\eta)}
[1] Crepieux et Bruno, [PRB 64, 014416 (2001)]
Args:
energies: Energy eigenvalues
mu: Chemical potential
temperature: Temperature
eta: Broadening
"""
# order (m, n)
enm = energies[:, None, :] - energies[:, :, None]
distribution = fermi_dirac(energies, mu, temperature)
fn = distribution[:, None, :]
fm = distribution[:, :, None]
kernel = 1j * (fm - fn) * invert(enm * (enm + 1j * eta))
return kernel
[docs]
def anomalous_kernel(
self,
energies: np.ndarray,
mu: float = 0.0,
temperature: float = 0.0,
eta: float = 0.0,
) -> np.ndarray:
r"""Computes the anomalous kernel.
.. math:
-2i\frac{f_n}{(E_n - E_m)(E_n - E_m + i\eta)}
Args:
energies: Energy eigenvalues
mu: Chemical potential
temperature: Temperature
eta: Broadening
"""
enm = energies[:, None, :] - energies[:, :, None]
fn = fermi_dirac(energies, mu, temperature)[:, None, :]
kernel = -2 * 1j * fn * invert(enm * (enm + 1j * eta))
return kernel
[docs]
def intraband_kubo_kernel(
self,
energies: np.ndarray,
mu: float = 0.0,
temperature: float = 0.0,
eta: float = 0.0,
) -> np.ndarray:
r"""Computes the intra-band kubo kernel
.. math:
\frac{i}{\eta}\partial_E f|_{E=e_n}
See for instance Eq(5), Huhtinen et al, [PRB 108, 155108 (2023)]
Args:
energies: Energy eigenvalues
mu: Chemical potential
temperature: Temperature
eta: Broadening
"""
if eta == 0:
raise ValueError("for the intra band contribution, eta cannot be 0")
return 1j * fermi_prime(energies, mu, temperature) / eta
[docs]
def oddKernel(
self,
energies: np.ndarray,
mu: float = 0.0,
temperature: float = 0.0,
eta: float = 0.0,
) -> np.ndarray:
r"""Eq (5) of [1]. The effect of temperature is included in the constant broadening :math:`\eta`.
Args:
energies: Energy eigenvalues
mu: Chemical potential
temperature: Temperature
eta: Broadening
[1] Freimuth et al, [PRB 90, 174423 (2014)]
"""
denom1 = (mu - energies[:, None, :]) ** 2 + eta**2
denom2 = (mu - energies[:, :, None]) ** 2 + eta**2
return eta**2 * invert(denom1 * denom2) / np.pi
[docs]
def evenKernel(
self,
energies: np.ndarray,
mu: float = 0.0,
temperature: float = 0.0,
eta: float = 0.0,
) -> np.ndarray:
r"""Eq (4) of [1]. The effect of temperature is included in the constant broadening :math:`\eta`.
Args:
energies: Energy eigenvalues
mu: Chemical potential
temperature: Temperature
eta: Broadening
[1] Freimuth et al, [PRB 90, 174423 (2014)]
"""
ediff = energies[:, :, None] - energies[..., None, :] # em - en
muEm = mu - energies[:, :, None] # mu - em
muEn = mu - energies[:, None, :] # mu - en
denom1 = muEn**2 + eta**2 # 1/[(mu-em)^2 + \Gamma^2]
denom2 = muEm**2 + eta**2 # 1/[(mu-en)^2 + \Gamma^2]
term1 = invert(denom1 * denom2)
term1 *= eta * ediff
term2 = invert(ediff * denom1)
term2 *= 2 * eta
term3 = np.log((-muEm - 1j * eta) * invert(-muEn - 1j * eta)).imag
term3 *= 2 * invert(ediff**2)
term1 += term2 + term3
return (1j / (2 * np.pi)) * term1
[docs]
def __call__(
self,
energies: np.ndarray,
mu: float = 0.0,
temperature: float = 0.0,
eta: float = 0.0,
) -> np.ndarray:
params = dict(energies=energies, mu=mu, temperature=temperature, eta=eta)
if self.name == "inter_band":
kernel = self.interband_kubo_kernel
if self.name == "anomalous":
kernel = self.anomalous_kernel
elif self.name == "intra_band":
kernel = self.intraband_kubo_kernel
elif self.name == "odd":
kernel = self.oddKernel
elif self.name == "even":
kernel = self.evenKernel
return kernel(**params)
[docs]
class KuboBastinKernel(EnergyKernel):
"""A class for Kubo-Bastin formula kernels.
Attributes:
name: Specifies the response to be computed.
"""
def __init__(self, name: str = None) -> None:
super().__init__()
self.name = name
if self.name is not None:
if name not in ["BMII", "BMI", "Overlap", "StredaI", "StredaII"]:
raise ValueError(
"name should be among [BMII, BMI, Overlap, StredaI, StredaII]"
)
[docs]
def bmII_kernel(
self,
energies: np.ndarray,
mu: float = 0.0,
temperature: float = 0.0,
eta: float = 0.0,
) -> np.ndarray:
r"""Sea contribution of Bonbian-Manchon decomposition
[1] Bonbian and Manchon [PRB 102, 085113 (2020)]
Args:
energies: Energy eigenvalues
mu: Chemical potential
temperature: Temperature
eta: Broadening
"""
# order (m, n)
emn = energies[:, :, None] - energies[:, None, :]
distribution = fermi_dirac(energies, mu, temperature)
fm = distribution[:, :, None]
i1 = invert((emn + 1j * eta) ** 2)
i2 = invert((emn - 1j * eta) ** 2)
kernel = 1j * fm * (i1 + i2)
return kernel
[docs]
def full_kernel(
self,
energies: np.ndarray,
mu: float = 0.0,
temperature: float = 0.0,
eta: float = 0.0,
) -> np.ndarray:
r"""Computes the Kubo Bastin kernel in Eq(A9) of [1]
[1] Crepieux et Bruno, PRB 64, 014416 (2001)
Args:
energies: Energy eigenvalues
mu: Chemical potential
temperature: Temperature
eta: Broadening
"""
# order (m, n)
enm = energies[:, :, None] - energies[:, None, :]
distribution = fermi_dirac(energies, mu, temperature)
fn = distribution[:, None, :]
kernel = -2 * 1j * (fn * invert((enm + 1j * eta) ** 2))
return kernel
[docs]
def stredaI_kernel(
self,
energies: np.ndarray,
mu: float = 0.0,
temperature: float = 0.0,
eta: float = 0.0,
) -> np.ndarray:
r"""First term of Strada-Smrka decomposition
Args:
energies: Energy eigenvalues
mu: Chemical potential
temperature: Temperature
eta: Broadening
"""
# order (m, n)
enm = energies[:, None, :] - energies[:, :, None]
distribution = fermi_prime(energies, mu, temperature)
fnprime = distribution[:, None, :]
kernel = -1j * fnprime * invert((enm + 1j * eta))
return kernel
[docs]
def overlap_kernel(
self,
energies: np.ndarray,
mu: float = 0.0,
temperature: float = 0.0,
eta: float = 0.0,
) -> np.ndarray:
r"""Overlap between surface and sea contributions
Args:
energies: Energy eigenvalues
mu: Chemical potential
temperature: Temperature
eta: Broadening
"""
# order (m, n)
enm = energies[:, None, :] - energies[:, :, None]
distribution = fermi_prime(energies, mu, temperature)
fnprime = distribution[:, None, :]
i1 = invert((enm + 1j * eta))
i2 = invert((enm - 1j * eta))
kernel = -1j * 0.5 * fnprime * (i1 + i2)
return kernel
[docs]
def __call__(
self,
energies: np.ndarray,
mu: float = 0.0,
temperature: float = 0.0,
eta: float = 0.0,
) -> np.ndarray:
params = dict(energies=energies, mu=mu, temperature=temperature, eta=eta)
if self.name is None:
return self.full_kernel(**params)
else:
if self.name == "StredaII":
return self.bmII_kernel(**params) - self.overlap_kernel(**params)
elif self.name == "BMI":
return self.stredaI_kernel(**params) - self.overlap_kernel(**params)
if self.name == "BMII":
kernel = self.bmII_kernel
elif self.name == "StredaI":
kernel = self.stredaI_kernel
elif self.name == "Overlap":
kernel = self.overlap_kernel
return kernel(**params)
[docs]
class SecondOrderKuboKernel(EnergyKernel):
def __init__(self) -> None:
super().__init__()
self.contractions = "knlm,knlm->"
[docs]
def __call__(
self,
energies: np.ndarray,
mu: float = 0.0,
temperature: float = 0.0,
eta: float = 0.0,
) -> np.ndarray:
# computes the non linear energy kernel:
# 0.5 * (1 / (enm + 2 * 1j * eta)) * (fnl / (enl + 1j * eta) - flm / (elm + 1j * eta))
distribution = fermi_dirac(energies, mu, temperature)
# order (k, n, m, l)
en = energies[:, :, None, None]
em = energies[:, None, :, None]
el = energies[:, None, None, :]
enm = invert(en - em + 2 * 1j * eta)
enl = invert(en - el + 1j * eta)
elm = invert(el - em + 1j * eta)
fn = distribution[:, :, None, None]
fm = distribution[:, None, :, None]
fl = distribution[:, None, None, :]
flm = fl - fm
flm = flm * elm
fnl = fn - fl
fnl = fnl * enl
fnl = fnl - flm
return 0.5 * fnl * enm