.. _transport:

====================================
Complete Example
====================================
    
Defining a symbolic hamiltonian
-------------------------------

We start by writing down the Rashba hamiltonian at interest. For this we will need sympy
.. code-block:: python

    from sympy import Matrix, cos, sin, symbols
    

We define our Rashba hamiltonian

.. math::
    H_R = 2 t - 2 t (\cos k_x + \cos k_y) + \Delta \sigma_z +  \alpha  (\sin k_y \sigma_x - \sin k_x  \sigma_y)

.. code-block:: python

    sigma_0 = Matrix([[1, 0], [0, 1]])
    sigma_x = Matrix([[0, 1], [1, 0]])
    sigma_y = Matrix([[0, -1j], [1j, 0]])
    sigma_z = Matrix([[1, 0], [0, -1]])
    def rashba_model():
        # Rashba model as a sympy expression
        k_x, k_y = symbols('k_x k_y')
        t, alpha, delta = symbols('t alpha delta')
        model = 2*t * sigma_0 - 2 * t * (cos(k_x) + cos(k_y)) * sigma_0 

        # add ferromagnetism
        model += delta * sigma_z

        # add spin-orbit coupling
        model += alpha * (sin(k_y) * sigma_x - sin(k_x) * sigma_y)
        return model

We transform this hamiltonian into a py4mulas Kmodel

.. code-block:: python

    from py4mulas.models import Kmodel
    model = rashba_model()
    kxs = np.linspace(-np.pi, np.pi, 100)
    kys = kxs
    params = dict(t=1, alpha=1.4, delta=0.2)

    # define a py4mulas model
    model = Kmodel(model, k_1=kxs, k_2=kys, params=params, real_space=True)

This defines a square momentum space hamiltonian with 100x100 kpoints.
The parameter real_space=True means that the hamiltonian is writtern in cartesian momenta.
In this case it does not matter as lattice momenta corresponds exactly to the cartesian coordinates.
If the system was constructed using lattice momenta, real_space should be set to False. Otherwise, 
the appropriate Bruillon zone should be considered. We recomand the use of the reduced momenta, in this case
one can always integrate over a torus :math:`[-\pi, \pi]^d`.

Computing the transport response
--------------------------------

Once the model is established we want to compute its transport response. 
We import the KuboBastinKernel which
enables to use the KuboBastin formula. 
We also import the main class Kubo which computes the Kubo response
with a given kernel.

.. code-block:: python

    from py4mulas.energy_kernels import KuboBastinKernel
    from py4mulas.responses import Kubo

We now define options that will be used in the computation. We want to use a precomputation of the response
in a way that changing the fermi energy, temperature or broadening does not rediagonalize the hamiltonian
We alse divide the momentum space into chunks, this helps optimize the computation. 

.. code-block:: python

    options = dict(precomp=True, chunk_size=10000)

Further, we compute the Kubo formula :math:`\sigma_{\alpha, \beta}`. 
The Kubo Bastin Kernals has 4 names which can be selected to specify the type of response we would like to compute.
This includes the streda decomposition kernels (StredaI) and (StredaII) or the Bonbien Manchon
kernels (BMI) and (BMII). If name is set to the default ``None`` the total KuboBastin kernel is called.

.. code-block:: python

    formula = Kubo(model, alpha='x', beta='y', kspace_options=options, kernel=KuboBastinKernel(name=None))

This is the transport object. Now it should be called for different transport parameters, here we 
compute the response while varying the chemical potential :math:`\mu` for fixed temperature and broadening

Plotting the output
-------------------

.. code-block:: python

    import numpy as np
    energies = np.linspace(-3, -1, 50)
    results = []
    for energy in energies:
        kwargs = dict(mu=energy, temperature=0.1, eta=0.1)
        results.append(formula(**kwargs))

now one can plot the results

.. code-block:: python

    import matplotlib.pyplot as plt
    plt.plot(energies, results)
    plt.show()

See examples for more options.

This sets a very simple example of transport calculation in Py4mulas.

Defining the model form Kwant
-----------------------------

Alternatively, the model can be provided as a kwant Builder which will be transformed into a symbolic model
using the kwant builder-to-model transformer. For details on how to do this see ``py4mulas/examples/haldane.py``

Computing in parallel
---------------------

If the kspace contains a large number of kpoints it is often suitable to perform the computation in parallel
leveraging available resources such as HPCs. For this purpose, py4mulas includes a submodule :mod:`~py4mulas.mpi.computers`

To use this facility we need to have mpi4py installed. We then import it and get the rank which will be useful
for materializing the output of the computation

.. code-block:: python

    import mpi4py
    comm = MPI.COMM_WORLD
    rank = comm.Get_rank()

The module :mod:`~py4mulas.mpi.computers` has two type of computers 

(i) :mod:`~py4mulas.mpi.computers.MuTEtaComputer`
which is designed to be efficeint when :math:`\mu`, :math:`T` or :math:`\eta` are varied in the simulation.

(ii) :mod:`~py4mulas.mpi.computers.HamParamComputer` which is designed to compute efficiently in parrallel
when one or more parameters of the Hamiltonian are changed.

Here, we illustrate the use of :mod:`~py4mulas.mpi.computers.MuTEtaComputer`.
To run the parallel computer one need to define an executor which treats the kspace related inputs

.. code-block:: python

    from py4mulas.mpi.computers import MuTEtaComputer
    from py4mulas.mpi.executors import KspaceExecutor
    executor = KspaceExecutor()

This executor can in principle handle different parallelization scenarios.  
Here we use its default functionality where it simply distributes the available kspace vectors equally to the workers

.. code-block:: python

    from py4mulas.mpi.executors import KspaceExecutor
    executor = KspaceExecutor()

Now to make the computation parallel we simply apply the computer to the formula with the executor

.. code-block:: python

    pFormula = MuTEtaComputer(formula, executor=executor)

The parallel formula object can take a list of chemical potentials, a list of temperatures or a list of 
broadening values. Or even a list for each one of them.
Here, we harness and plot the resuls.

 .. code-block:: python   

    result = pFormula(mu=energies, temperature=[0], eta=[0])
    if rank==0:
        plt.plot(energies, result)
        plt.show()

if all parameters are given as lists or 1d arrays the output array result will be ordered according to 
to the ordering of the product of these lists. 
Internally py4mulas uses 

.. code-block:: python  

    itertools.product(mu, temperature, eta)

This order is deterministic and should be respected. Otherwise the result will be read in a wrong manner.