{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import sisl as si\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Anomalous Hall conductivity (AHC) for graphene\n",
    "\n",
    "This tutorial will describe a complete walk-through of how to calculate the anomalous Hall conductivity for graphene.\n",
    "\n",
    "## Creating the geometry to investigate\n",
    "\n",
    "Our system of interest will be the pristine graphene system, from a DFT (SIESTA) calculation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "H = si.get_sile(\"siesta_2/RUN.fdf\").read_hamiltonian()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The anomalous Hall conductivity (AHC) requires a rather dense $k$ grid. It is related to the Berry curvature defined as:\n",
    "$$\n",
    "    \\boldsymbol\\Omega_{i,\\alpha\\beta} = 2i\\hbar^2\\sum_{j\\neq i}\n",
    "                \\frac{\\hat v^{\\alpha}_{ij} \\hat v^\\beta_{ji}}\n",
    "                     {[\\epsilon_j - \\epsilon_i]^2 + i\\eta^2}\n",
    "$$\n",
    "where $\\hat v$ is the velocity operator. One can determine that the units of this quantity is $\\mathrm{Ang}^2$.\n",
    "The AHC can then be calculated via:\n",
    "$$\n",
    "    \\sigma_{\\alpha\\beta} = \\frac{-e^2}{\\hbar}\\int\\,\\mathrm d\\mathbf k\\sum_i f_i\\Omega_{i,\\alpha\\beta}(\\mathbf k).\n",
    "$$\n",
    "This method is implemented in `sisl.physics.electron.ahc`.\n",
    "The units of AHC is $\\mathrm S / \\mathrm{Ang}^{2 - D}$ which for 2D systems is just $\\mathrm S$.  \n",
    "Its API looks like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "help(si.physics.electron.ahc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "-----\n",
    "\n",
    "We will be interested in calculating the `ahc` for a set of different energies (equivalent to different chemical potentials).\n",
    "So we need to define an energy-range, and a distribution function, here the simple step-function is used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "E = np.linspace(-5, 2, 251)\n",
    "# When calculating for a variety of energy-points, we have to have an available axis for the eigenvalue distribution\n",
    "# calculation.\n",
    "dist = si.get_distribution(\"step\", x0=E.reshape(-1, 1))\n",
    "# Generally you want a *very* dense k-point grid\n",
    "bz = si.MonkhorstPack(H, [33, 33, 1], trs=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since we are only interested in the $xy$ plane (there is no periodicity along $z$, hence superfluous calculation), we will try and reduce the computation by specifying which axes we want to calculate the AHC along.\n",
    "Additionally, we can speed up the calculation for matrices that are small, by explicitly specifying them to be calculated in the `numpy.ndarray` format, as opposed to the default `scipy.sparse.csr_matrix` format (slower, but much less memory consuming).\n",
    "When dealing with these conductivities it can also be instructive to view the $k$-resolved conductivities at certain energies.\n",
    "\n",
    "\n",
    "Lastly, `eta=True` specifies we want to show the progressbar.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ahc = si.physics.electron.ahc(\n",
    "    bz,\n",
    "    # yield a k-resolved AHC array\n",
    "    k_average=False,\n",
    "    distribution=dist,\n",
    "    # Speed up by format='array'\n",
    "    eigenstate_kwargs={\"format\": \"array\", \"dtype\": np.complex128, \"eta\": True},\n",
    "    # Speed up by only calculating the xx, xy, yx, yy contributions\n",
    "    derivative_kwargs={\"axes\": \"xy\"},\n",
    ").real  # we don't need the imaginary part."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have a `(len(bz), 2, 2, len(E))` AHC array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(E, ahc.sum(0)[0, 1], label=r\"$\\sigma_{xy}$\")\n",
    "plt.plot(E, ahc.sum(0)[1, 0], label=r\"$\\sigma_{yx}$\")\n",
    "plt.xlabel(\"Energy [eV]\")\n",
    "plt.ylabel(r\"$\\sigma$ [Ang^2]\")\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see there are some contributions near $E_F$, and at the lower bands.  \n",
    "The numbers though are not converged, and likely requires more scrutiny (outside the scope of this tutorial).\n",
    "\n",
    "We can now plot the $k$-resolved AHC values at $E_F$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "E0 = np.argmin(np.fabs(E))\n",
    "kx = np.unique(bz.k[:, 0])\n",
    "ky = np.unique(bz.k[:, 1])\n",
    "plt.contourf(\n",
    "    # unique values along x, y\n",
    "    kx,\n",
    "    ky,\n",
    "    ahc[:, 0, 1, E0].reshape(len(kx), len(ky)),\n",
    ")\n",
    "plt.colorbar();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This concludes a simple tutorial on how to calculate the AHC for a given  system, and also how to calculate the $k$-resolved AHC."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}