Build Your Own Formula
Py4mulas implements a number of built-in responses, which are the common formulas we use for computing
transport responses (See. Available Responses). However, one may sometimes need to define our own formula, which may not belong to the most
popular response formulas used in the field.
No worries, Py4mulas is designed to enable this construction, in a flexible manner.
For this advanced construction of formulae, the user should be familiar with basic numpy operations, as all formulae are required to be implemented in a vectorized manner.
Here, we go. To construct a formula, we mainly need to define two kernels. A kernel which maybe stored, that is, which does not depend on transport parameters such as \(\mu\), \(T\) and \(\eta\). In this way changing any of these parameters does not require re-diagonalizing the model’s Hamiltonian. This makes the computation faster. Here we present an example of formula to be constructed.
We want to compute the conductivity
We denote \(L_{mn}=\langle m|\beta|n \rangle\), \(R_{nm} = \langle n|\alpha|m \rangle\)
We then write it as a product of a \(\mu\) dependent kernel (mu_kernel) \(K\)
and another kernel \(O\) which involves Py4mulas operators. This is what we call an opera_kernel
with \(O = L R^T\)
and
Using np.einsum, this formula can be simply written as :
sigma = np.einsum('kij,kij->', O, K)
To implement this formula, we shall define the O kernel and the K kernel, with the appropriate contractions. Here, is how to proceed.
In py4mulas the O kernel is called opera_kernel and the K kernel is in principle (\(\mu\), \(T\), \(\eta\)) dependent, but
for short notation it is called mu_kernel.
Defining an opera_kernel
First, we define, the O kernel, which is required to be a callable with signature (k_args, energy, psi, eta).
It can be either a class or a function.
from py4mulas.operators import Velocity
class OperaKernel:
def __init__(self, model, alpha, beta):
self.alpha = alpha
self.beta = beta
if isinstance(alpha, str):
self.alpha = Velocity(model, alpha)
if isinstance(beta, str):
self.beta = Velocity(model, beta)
def __call__(self, k_args, energy, psi, eta):
alpha = self.alpha(k_args, energy, psi)
beta = self.beta(k_args, energy, psi)
kernel = beta * np.swapaxes(alpha, 1, 2)
return kernel
opera_kernel = OperaKernel(some_model, alpha='x', beta='y')
This kernel could be also defined as a function
def opera_kernel_fun(model, alpha, beta):
if isinstance(alpha, str):
alpha = Velocity(model, alpha)
if isinstance(beta, str):
beta = Velocity(model, beta)
def fun(k_args, energy, psi, eta):
alpha = alpha(k_args, energy, psi)
beta = beta(k_args, energy, psi)
kernel = beta * np.swapaxes(alpha, 1, 2)
return kernel
return fun
opera_kernel = opera_kernel_fun(some_model, alpha='x', beta='y')
Here, we have used alpha and beta as str s, but this kernel could take any py4mulas operator.
Note that operators always have signature (k_args, energy, psi).
Defining a mu_kernel
Now we shall define the second kernel which may depend on transport parameters.
class KuboKernel:
def __init__(self) -> None:
self.contractions = "kij,kij->"
def __call__(self, energy, mu, temperature, eta):
# order (m, n)
enm = energy[:, None, :] - energy[:, :, None]
distribution = fermi_dirac(energy, mu, temperature)
fn = distribution[:, None, :]
fm = distribution[:, :, None]
kernel = 1j * (fm - fn) * invert(enm * (enm + 1j * eta))
return kernel
mu_kernel = KuboKernel()
mu_kernel s should be callable with signature (energy, mu, temperature, eta).
Optionally a mu_kernel may implement an attribute contractions, as above. If this attribute is not implemented or if the kernel is
defined as a function, contractions should be passed as an argument to the formula.
This above mu_kernel is implemented also in KuboKernel.
So it can be also accessed as
from py4mulas.mu_kernels import KuboKernel
mu_kernel = KuboKernel("inter_band")
Here, the simple case of kernels is considered, where we have one mu_kernel and one opera_kernel.
In practice, one may need to implement a formula made of a sum of kernel products with different contractions.
In such a case both mu_kernel and opera_kernel can be passed as lists of kernels.
opera_kernel = 2 * [OperaKernel(some_model, alpha='x', beta='y')]
mu_kernel = 2 * [KuboKernel("inter_band")]
In this case each opera_kernel is multiplied by the corresponding mu_kernel.
Taking these kernels together, the formula is :
from py4mulas.formulas import KuboFormula
formula = KuboFormula(some_model, kspace_options=some_opteions,
contractions=contractions, mu_kernel=mu_kernel, opera_kernel=opera_kernel)
In the case above, the response would simply double.
See. advancedKernels.py for a complete example.
Note
It is important to note that contractions, opera_kernel and mu_kernel should have the same length.
And should be given in the same order.
We recommand the definition of a formula in this manner, by passing kernels as arguments.
Alternatively, one can write a class which inherit from KuboFormula. However,
in this case, one should implement _mu_kernel() and _opera_kernel()
methods.
Note that a mu_kernel should necessarily implement the attribute contractions
which is used to compute the response using np.einsum, as illustrated above.
contractions corresponds to np.einsum subscripts.
These should always start with k, which is reserved to the momentum summation symbole.
Further, contractions should be defined in a way consistent to the shape of the kernels.
Actually, Py4mulas intercepts
these contractions to deduce the assumed shapes of the kernels, entring the formula.
Note also that the mu_kernel should be callable with signature (energy, mu, temperature, eta).
Defining kernels as functions
Alternatively, one can define both mu_kernel and opera_kernel as functions, which return callables
or lists of callables. Here, is an example
def kernel1(energy, mu, temperature, eta):
pass
def kernel2(energy, mu, temperature, eta):
pass
mu_kernel = [kernel1, kernel2]
def fun_kernel(*args, **kwargs):
def fun1(k_args, en, psi, eta):
pass
def fun2(k_args, en, psi, eta):
pass
return [fun1, fun2]
opera_kernel = fun_kernel
Onece the formula is defined, all computation strategies can be applied including prallelization (see. Summary Of Simulation Steps and Complete Example), precomputation and memeory use optimization (see. Setting Kspace Options).
Defining Contractions
An important step to implement a formula is the definition of appropriate contractions.
In general for a model with n orbitals, kernel matrices should have the shape (nk, n, n)
with nk the number of k points. Often we are left with a product of 2 arrays each of this shape,
the corresponding contractions is therefore kij,kij-> which means a summation over k, i and j, resulting in a scalar.
The first symbole of contractions is necessary to be denoted as k for momentum summation, the other symbols are free.
One could use kmn,kmn->. This contraction obviously leads to a scalar. However, one might sometimes
need to have access to k_resolved quantities calling our foumula with k_resolved=True.
In this case Py4mulas automatically updates the contractions to kmn,kmn->k, so it refreins from
summing over k.
If the contractions is given with symbols after the arrow, Py4mulas assumes
that contractions is appropriate for the underlying kernels, so it does not update the subscripts
even if k_resolved is set to True. In this case we should be careful
and ensure that contractions is compatible with kernel shapes.
See for instance the BerryKernel where
a contractions of the form kij,kij->ki is defined for momentum and band resoved Berry curvature.
Having imposed the momentum space symbole to always be first in contractions, the latter should be always
given in the explicit mode, that is it should always contain the header '->' and specify the output if some indices would remain after summation.
The implicit convention is not allowed.
For more on defining np.einsum subscripts,
see Numpy documentation.
On Kernel Parameters
Above, we have discussed how to construct a formula in a general manner.
Particularly, when we are interested in varying the parameter \(\eta\),
the way we construct the formula matters for performance.
In fact, if the opera_kernel is \(\eta\) dependent
the use of the precomputation strategy would not help,
as Py4mulas would recompute the kernel (which is supposed to be frozen) whenever \(\eta\) is varied.
See. Setting Kspace Options. In this case, it is more interesting to define the kernels such that the kernel to be stored
(the opera_kernel) is \(\eta\) independent. This means that mu_kernel will then depend on \(\eta\).
One could in fact, make \(\eta\) always belong to the mu_kernel. However, in most cases treating
the \(\eta\) dependence in opera_kernel comes with an important reduction of the kernel matrix dimensions, in a manner
that drastically improves performance when \(\mu\) and/or \(T\) are to be varied.
To illustrate this in depth, let us analyse the Kubo formula written above. One can notice that the formula can be simplified as
where
and
Now the kernel to store is dimension reduced with the cost of transferring the \(\eta\) dependence
into opera_kernel. Therefore, this kernel cannot be stored if \(\eta\) is being varied.
Nonetheless, varying
\(\mu\) and/or \(T\) while using this decomposition is much more efficient than using the initial Kubo formula.
In our implementation of the Kubo formula this is appropriately handled as we enable to either provide an external \(\eta\)
dependent mu_kernel, which makes the computation efficient for \(\eta\) varied responses.
At the same time we keep
the benefit of the dimension reduction of the opera_kernel matrix when no external kernel is provided. In the latter case
the mu_kernel would simply correspond to the Fermi-Dirac distribution FermiDirac which is naturally \(\eta\) independent.
Here is how we would proceed using Py4mulas built-in Kubo response with the reduced kernel.
from py4mulas.responses import Kubo
options = dict(precomp=True, chunk_size=10000)
G = Kubo(some_model, alpha="x", beta="y", kspace_options=options, mu_kernel=None)
This is the most efficitent setting if \(\eta\) is fixed in the computation.
Note taht setting mu_kernel=None is equivalent to setting
from py4mulas.mu_kernels import FermiDirac
mu_kernel = FermiDirac()
Now, if we aim at varying \(\eta\)
it is more efficient to use an \(\eta\) dependent mu_kernel
from py4mulas.mu_kernels import KuboKernel
kernel = KuboKernel('inter_band')
G = Kubo(some_model, alpha="x", beta="y", kspace_options=options, mu_kernel=kernel)
In the latter case the opera_kernel to be stored has a shape (nk, n, n)
instead of (nk, n) for the reduced kernel.
Note that, using the reduced (or symmetrized) Kubo formula, contractions becomes simply km,km-> instead of kmn,kmn->
corresponding to a reduction of an axis of the opera_kernel array.
Notes On Vectorization
It is important to notice that Py4mulas requires the kernels to accept
arrays for energy and psi.
These arrays are results of hamiltonian diagonalization which is performed internally. The argument energy would have
the shape (nk, n) and psi will be of shape (nk, n, n).
These are ouputs of np.linalg.eigh.
For instance to compute the energy difference \(E_n-E_m\) in the mu_kernel we augment the dimensions, resulting in an
array of shape (nk, n, n)
enm = energy[:, None, :] - energy[:, :, None]
Note that the order here is important as the kernel is defined as K_{mn}, m is intrpreted as row index and
n as column index. This matters when creating newaxes.
As we are interested in the inter-band contribution the diagonal terms are always zero using the invert
function which implements a safe division.
For the opera_kernel we also need to take psi as an argument which finally goas to the relevent operators.
Note that all Py4mulas operators have the signature (k_args, energy, psi, eta).
Prior to using a kernel, we always recomand to test the vectorized expression with a strightforward unvectorized implementation.