.. _advanced: ==================================== Build Your Own Formula ==================================== ``Py4mulas`` implements a number of built-in responses, which are the common formulas we use for computing transport responses (See. :ref:`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 :math:`\mu`, :math:`T` and :math:`\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 .. 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. We denote :math:`L_{mn}=\langle m|\beta|n \rangle`, :math:`R_{nm} = \langle n|\alpha|m \rangle` We then write it as a product of a :math:`\mu` dependent kernel (``mu_kernel``) :math:`K` and another kernel :math:`O` which involves ``Py4mulas`` operators. This is what we call an ``opera_kernel`` .. math:: \sigma_{\alpha,\beta}=\sum_k \sum_{mn}K_{mn} O_{mn} with :math:`O = L R^T` and .. math:: K_{mn} = \frac{f_m - f_n}{(E_n - E_m) (E_n - E_m + i \eta)}. Using ``np.einsum``, this formula can be simply written as : .. code-block:: python 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 (:math:`\mu`, :math:`T`, :math:`\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``. .. code-block:: python 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`` .. code-block:: python 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. .. code-block:: python 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 :class:`~py4mulas.mu_kernels.KuboKernel`. So it can be also accessed as .. code-block:: python 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. .. code-block:: python 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 : .. code-block:: python 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. :download:`advancedKernels.py <../../../examples/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 :class:`~py4mulas.formulas.KuboFormula`. However, in this case, one should implement :meth:`~py4mulas.formulas.KuboFormula._mu_kernel` and :meth:`~py4mulas.formulas.KuboFormula._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 .. code-block:: python def kernel1(energy, mu, temperature, eta): pass def kernel2(energy, mu, temperature, eta): pass mu_kernel = [kernel1, kernel2] .. code-block:: python 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. :ref:`summary` and :ref:`transport`), precomputation and memeory use optimization (see. :ref:`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 :class:`~py4mulas.mu_kernels.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 :math:`\eta`, the way we construct the formula matters for performance. In fact, if the ``opera_kernel`` is :math:`\eta` dependent the use of the precomputation strategy would not help, as ``Py4mulas`` would recompute the kernel (which is supposed to be frozen) whenever :math:`\eta` is varied. See. :ref:`options`. In this case, it is more interesting to define the kernels such that the kernel to be stored (the ``opera_kernel``) is :math:`\eta` independent. This means that ``mu_kernel`` will then depend on :math:`\eta`. One could in fact, make :math:`\eta` always belong to the ``mu_kernel``. However, in most cases treating the :math:`\eta` dependence in ``opera_kernel`` comes with an important reduction of the kernel matrix dimensions, in a manner that drastically improves performance when :math:`\mu` and/or :math:`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 .. math:: \sigma_{\alpha,\beta} = \sum_k\sum_{m} f_m V_m, where .. math:: V_m = \sum_n (\Gamma_{mn} - \Gamma_{nm}) and .. math:: \Gamma_{mn} = \frac{O_{mn}}{(E_n - E_m) (E_n - E_m + i \eta)}. Now the kernel to store is dimension reduced with the cost of transferring the :math:`\eta` dependence into ``opera_kernel``. Therefore, this kernel cannot be stored if :math:`\eta` is being varied. Nonetheless, varying :math:`\mu` and/or :math:`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 :math:`\eta` dependent ``mu_kernel``, which makes the computation efficient for :math:`\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 :class:`~py4mulas.mu_kernels.FermiDirac` which is naturally :math:`\eta` independent. Here is how we would proceed using ``Py4mulas`` built-in Kubo response with the reduced kernel. .. code-block:: python 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 :math:`\eta` is fixed in the computation. Note taht setting ``mu_kernel=None`` is equivalent to setting .. code-block:: python from py4mulas.mu_kernels import FermiDirac mu_kernel = FermiDirac() Now, if we aim at varying :math:`\eta` it is more efficient to use an :math:`\eta` dependent ``mu_kernel`` .. code-block:: python 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 :math:`E_n-E_m` in the ``mu_kernel`` we augment the dimensions, resulting in an array of shape ``(nk, n, n)`` .. code-block:: python 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.