from typing import Optional, Union
import numpy as np
from .models import Kmodel
from .operators import KspaceOpera, Velocity
from .formulas import KuboFormula
from .energy_kernels import FermiDirac, EnergyKernel, SecondOrderKuboKernel
from ._common import invert
__all__ = ["Kubo", "SecondOrderKubo"]
[docs]
class Kubo(KuboFormula):
r"""Computes the kubo conductivity
.. math::
\sigma_{\alpha,\beta} = i\sum_k \sum_{mn} \frac{f_m - f_n}{(E_n - E_m) (E_n - E_m + i \eta)} \langle m|\beta|n \rangle \langle n|\alpha|m \rangle.
:math:`R_{nm} = \langle n|\alpha|m \rangle`.
We note :math:`L_{mn}=\langle m|\beta|n \rangle`, :math:`R_{nm} = \langle n|\alpha|m \rangle`
We then write it as a product of an energy kernel :math:`K`
and operator matrix :math:`O`
.. math::
\sigma_{\alpha,\beta}=\sum_k \sum_{mn}K_{mn} O_{mn}
with :math:`O = L R^T`
and :math:`K_{mn} = (f_m - f_n) I_{mn}.`
Where :math:`I^{mn} = \frac{1}{(E_n - E_m) (E_n - E_m + i \eta)}.`
We enable the construction of any kubo like formula by providing a tailored :math:`K`
kernel can depend on :math:`(E, \mu, T, \eta).`
In the case of the habitual kubo formula, the computation can be further simplified as
.. math::
\sigma_{\alpha,\beta} = \sum_k\sum_{m} f_m V_m
with :math:`V_m = \sum_n (\Gamma_{mn} - \Gamma_{nm})`
and :math:`\Gamma_{ij} = I_{ij} O_{ij}`.
Note:
this symmetrized decomposition is more memory efficient than the kernel based approach,
as it holds a reduced array.
Attributes:
model: An instance of :class:`~py4mulas.models.Kmodel`
alpha: An instance of :class:`~py4mulas.operators.KspaceOpera`
beta: An instance of :class:`~py4mulas.operators.KspaceOpera`
kspace_options: a dictionary specifying chunk_size and precomp.
If chunk_size is not given, a default is used.
Note: if precomp is True, the energy and temperature scans are momentum space free.
If eta is to be changed, then set precomp = False, as the kernels are eta-dependent.
kernel: An instance or list of instances of :class:`~py4mulas.energy_kernels.EnergyKernel`,
the energy kernel of the Kubo formula.
If a list of kernels is given they should assume the same shape
as :func:`main_kernel`, when precomp is `True`.
If precomp is `False` then the shape uniformity is not important.
Exemple:
>>> G = Kubo(some_model, alpha='x', beta=Velocity('y'),
kspace_options=dict(chunk_size=1000, precomp=True), kernel=EnergyKernel('inter_band'))
>>> conductivities = []
>>> for mu_i in np.arange(-1, 1):
>>> conductivity = G(mu=mu_i, temperature=0, eta=0)
>>> conductivities.append(conductivity)
"""
def __init__(
self,
model: Kmodel,
kspace_options: Optional[dict] = None,
alpha: Union[str, KspaceOpera] = "x",
beta: Union[str, KspaceOpera] = "y",
kernel: Optional[Union[EnergyKernel, list[EnergyKernel]]] = None,
) -> None:
if kernel is None:
# native interband kubo formula
self.contractions = "ki,ki->"
else:
if isinstance(kernel, list):
# these contractions should assume the same main_kernel
# so any one of them is a good representative
self.contractions = kernel[0].contractions
else:
# a single energy_kernel => a single contractions anyway
self.contractions = kernel.contractions
super().__init__(
model, kspace_options=kspace_options, contractions=self.contractions
)
_inspect_opera([alpha, beta])
self.alpha = alpha
if isinstance(alpha, str):
self.alpha = Velocity(model, alpha)
self.beta = beta
if isinstance(beta, str):
self.beta = Velocity(model, beta)
self.kernel = kernel
[docs]
def energy_kernel(self):
if self.kernel is None:
kernel = FermiDirac()
kernel.contractions = "ki,ki->"
return [kernel]
kernel = self.kernel
return kernel if isinstance(kernel, (list, tuple)) else [kernel]
[docs]
def main_kernel(self, k_args, eta):
hk = self.H(*k_args)
en, psi = np.linalg.eigh(hk)
alpha = self.alpha(k_args, en, psi)
opera_kernel = self.beta(k_args, en, psi)
opera_kernel *= np.swapaxes(alpha, 1, 2)
if self.kernel is not None:
return opera_kernel, en
delta_e = en[:, None, :] - en[:, :, None]
denom = delta_e * (delta_e + 1j * eta)
opera_kernel *= 1j * invert(denom)
# reduce into a vector for symmetrized Kubo formula
reduced = opera_kernel.sum(axis=2) - opera_kernel.sum(axis=1)
return reduced, en
[docs]
class SecondOrderKubo(Kubo):
r"""Computes the kubo conductivity
.. math::
\sigma_{\alpha,\beta \gamma} = \frac{1}{2}\sum_k \sum_{mnl} O^{nml}(k) K_{nml}(k).
with :
.. math:: O^{nml} = \langle n|\alpha|m \rangle
(\langle m|\beta|l \rangle \langle l|\gamma|n \rangle + \beta \leftrightarrow \alpha)
.. math:: K_{nml} = \frac{1}{E_n - E_m + 2 i \eta}
(\frac{f_n - f_l}{E_n - E_l + i \eta} - \frac{f_l - f_m}{E_l - E_m + i \eta})
Terms corresponding to the decomposition of this formula into odd and even
under time-reversal symmetry are implemented. See :meth:`main_kernel`.
See. Yatsushiro et al [PRB 104, 054412 (2021)].
Note:
This class can in principle implement any energy kernel :math:`K` provided that it assumes
an operator kernel which is exactly equal to :math:`O` described above.
For formulas containing combinations of operator kernels, each product of the type :math:`K O`
can be implemented independently from formulas.
In this class we also enable `even` and `odd` energy kernels which involve different operator
kernels.
These can be accessed through the parameter ``contribution``.
Attributes:
model: :class:`~py4mulas.models.Kmodel`
alpha: :class:`~py4mulas.operators.KspaceOpera`
beta: :class:`~py4mulas.operators.KspaceOpera`
gamma: :class:`~py4mulas.operators.KspaceOpera`
kspace_options: A dictionary specifying chunk_size and precomp.
If chunk_size is not given, a default is used.
Note: if precomp is True, the energy and temperature scans are momentum space free.
If eta is to be changed, then set precomp = False, as the kernels are eta-dependent.
kernel: :class:`~py4mulas.energy_kernels.SecondOrderKuboKernel`,
the second order energy kernel of the Kubo formula.
If a list of kernels is given they should assume the same shape
as :func:`main_kernel`, when precomp is `True`.
If precomp is `False` then the shape uniformity is not important.
contribution: Specifies a particular kernel from ['even', 'odd']. When this is set to None
the full undecomposed kernel is assumed.
`contribution=None` is equivalent to kernel=SecondOrderKuboKernel().
Exemple:
>>> G = Kubo(some_model, alpha='x', beta=Velocity('y'), gamma='z',
kspace_options=dict(chunk_size=10000, precomp=True), kernel=SecondOrderKuboKernel())
>>> conductivities = []
>>> for mu_i in np.arange(-1, 1):
>>> conductivity = G(mu=mu_i, temperature=0, eta=0)
>>> conductivities.append(conductivity)
"""
def __init__(
self,
model: Kmodel,
kspace_options: Optional[dict] = None,
alpha: Union[str, KspaceOpera] = "x",
beta: Union[str, KspaceOpera] = "y",
gamma: Union[str, KspaceOpera] = "x",
kernel: Union[SecondOrderKuboKernel, list[SecondOrderKuboKernel]] = None,
contribution: str = None,
) -> None:
super().__init__(
model, alpha=alpha, beta=beta, kspace_options=kspace_options, kernel=kernel
)
_inspect_opera([gamma])
self.gamma = gamma
if isinstance(gamma, str):
self.gamma = Velocity(model, gamma)
if contribution is not None and kernel is not None:
raise ValueError("kernel and contribution are mutually exclusive")
if contribution is not None:
if contribution not in ["even", "odd"]:
raise ValueError("contribution should be either even or odd")
self.contribution = contribution
[docs]
def main_kernel(
self, k_args: list[np.ndarray], eta: float
) -> tuple[np.ndarray, np.ndarray]:
hk = self.H(*k_args)
energy, psi = np.linalg.eigh(hk)
alpha = self.alpha(k_args, energy, psi)
beta = self.beta(k_args, energy, psi)
gamma = self.gamma(k_args, energy, psi)
if self.kernel is not None:
# opera_kernel = Bi_nm * (Aj_ml * Ak_ln + Ak_ml * Aj_ln)
opera_kernel = np.einsum(
"inm,iml,iln->inml", alpha, beta, gamma, optimize=True
)
opera_kernel += np.einsum(
"inm,iml,iln->inml", alpha, gamma, beta, optimize=True
)
return opera_kernel, energy
enm = energy[:, :, None] - energy[:, None, :]
if self.contribution is None:
opera_kernel = self._full_kernel(enm, alpha, beta, gamma, eta)
elif self.contribution == "even":
opera_kernel = self._even_kernel(enm, alpha, beta, gamma, eta)
elif self.contribution == "odd":
opera_kernel = self._odd_kernel(enm, alpha, beta, gamma, eta)
return opera_kernel, energy
def _full_kernel(
self,
enm: np.ndarray,
alpha: np.ndarray,
beta: np.ndarray,
gamma: np.ndarray,
eta: float,
) -> np.ndarray:
S = invert(enm + 1j * eta)
alpha *= invert(enm + 2j * eta)
reduced = _kernel_reducer(S, alpha, beta, gamma)
reduced += _kernel_reducer(S, alpha, gamma, beta)
return 0.5 * reduced
def _even_kernel(
self,
enm: np.ndarray,
alpha: np.ndarray,
beta: np.ndarray,
gamma: np.ndarray,
eta: float,
) -> np.ndarray:
# ordered (k, n, m, l)
S = invert(enm**2 + eta**2)
alpha *= invert(enm**2 + 4 * eta**2)
reduced = _kernel_reducer(S, alpha, beta, gamma)
reduced *= -2 * eta**2
S *= enm
alpha *= enm
reduced += _kernel_reducer(S, alpha, beta, gamma)
return reduced
def _odd_kernel(
self,
enm: np.ndarray,
alpha: np.ndarray,
beta: np.ndarray,
gamma: np.ndarray,
eta: float,
) -> np.ndarray:
S2 = invert(enm**2 + eta**2)
S1 = enm * S2
alpha *= invert(enm**2 + 4 * eta**2)
reduced = 2 * eta * _kernel_reducer(S1, alpha, beta, gamma)
alpha *= enm
reduced += eta * _kernel_reducer(S2, alpha, beta, gamma)
return -1j * reduced
def _kernel_reducer(
S: np.ndarray, alpha: np.ndarray, beta: np.ndarray, gamma: np.ndarray
) -> np.ndarray:
r"""Computes the operator kernel
.. math:: \sum_{ml} (O^{nml} S_{nl} - O^{lmn} (S_{ln} - S_{nm}) + O^{mnl} S_{ln})
where:
.. math:: O_{nml} = \frac{\alpha^{nm}\beta_{ml}\gamma^{ln}}{E_n - E_m + 2 i \eta}.
.. math:: S_{nm} = \frac{1}{E_n - E_m + i \eta}
Returns:
The contracted array in shape (nk, norbs) to be stored.
"""
reduced = np.einsum(
"knm, kml, kln, knl -> kn", alpha, beta, gamma, S, optimize=True
)
reduced -= np.einsum(
"klm, kmn, knl, kln -> kn", alpha, beta, gamma, S, optimize=True
)
reduced -= np.einsum(
"klm, kmn, knl, knm -> kn", alpha, beta, gamma, S, optimize=True
)
reduced += np.einsum(
"kmn, knl, klm, kln -> kn", alpha, beta, gamma, S, optimize=True
)
return reduced
def _inspect_opera(operas: list) -> None:
allowed_opera = (str, KspaceOpera)
for opera in operas:
if not isinstance(opera, allowed_opera):
raise TypeError(f"Unsupported operator type: {type(input)}")
if isinstance(opera, str):
if opera not in "xyz":
raise ValueError("Velocity operator should be x, y or z")