{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from sisl import *\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Berry phase calculation for graphene\n",
    "\n",
    "This tutorial will describe a complete walk-through of how to calculate the Berry phase for graphene.\n",
    "\n",
    "## Creating the geometry to investigate\n",
    "\n",
    "Our system of interest will be the pristine graphene system with the on-site terms shifted by $\\pm\\delta$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "graphene = geom.graphene()\n",
    "H = Hamiltonian(graphene)\n",
    "H.construct([(0.1, 1.44), (0, -2.7)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`H` now contains the pristine graphene tight-binding model. The anti-symmetric Hamiltonian is constructed like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "H_bp = H.copy()  # an exact copy\n",
    "H_bp[0, 0] = 0.1\n",
    "H_bp[1, 1] = -0.1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparing electronic structures\n",
    "\n",
    "Before proceeding to the Berry phase calculation lets compare the band structure and DOS of the two models. The anti-symmetric Hamiltonian opens a gap around the Dirac cone. A zoom on the Dirac cone shows this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "band = BandStructure(\n",
    "    H, [[0, 0.5, 0], [1 / 3, 2 / 3, 0], [0.5, 0.5, 0]], 400, [r\"$M$\", r\"$K$\", r\"$M'$\"]\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "band.set_parent(H)\n",
    "plot = band.plot(Erange=[-3, 3])\n",
    "band.set_parent(H_bp)\n",
    "plot.merge(band.plot(Erange=[-3, 3]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The gap opened is equal to the difference between the two on-site terms. In this case it equals $0.2\\mathrm{eV}$. Lets, for completeness sake calculate the DOS close to the Dirac point for the two systems. To resolve the gap the distribution function (in this case the Gaussian) needs to have a small smearing value to ensure the states are not too spread and the gap smeared out."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bz = MonkhorstPack(\n",
    "    H, [41, 41, 1], displacement=[1 / 3, 2 / 3, 0], size=[0.125, 0.125, 1]\n",
    ")\n",
    "bz_average = bz.apply.average  # specify the Brillouin zone to perform an average"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above `MonkhorstPack` grid initialization is creating a Monkhorst-Pack grid centered on the $K$ point with a reduced Brillouin zone size of $1/8$th of the entire Brillouin zone. Essentially this *only* calculates the DOS in a small $k$-region around the $K$-point. Since in this case we know the electronic structure of our system we can neglect all contributions from $k$-space away from the $K$-point because we are only interested in energies close to the Dirac-point.\n",
    "Here the sampled $k$-points are plotted. Note how they are concentrated around $[1/3, -1/3]$ which is the $K$-point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(bz.k[:, 0], bz.k[:, 1], 2)\n",
    "plt.xlabel(r\"$k_x$ [$b_x$]\")\n",
    "plt.ylabel(r\"$k_y$ [$b_y$]\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before proceeding to the Berry phase calculation we calculate the DOS in an energy region around the Dirac-point to confirm the band-gap."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "E = np.linspace(-0.5, 0.5, 1000)\n",
    "dist = get_distribution(\"gaussian\", 0.03)\n",
    "bz.set_parent(H)\n",
    "plt.plot(\n",
    "    E,\n",
    "    bz_average.eigenvalue(wrap=lambda ev: ev.DOS(E, distribution=dist)),\n",
    "    label=\"Graphene\",\n",
    ")\n",
    "bz.set_parent(H_bp)\n",
    "plt.plot(\n",
    "    E,\n",
    "    bz_average.eigenvalue(wrap=lambda ev: ev.DOS(E, distribution=dist)),\n",
    "    label=\"Graphene anti\",\n",
    ")\n",
    "plt.legend()\n",
    "plt.ylim([0, None])\n",
    "plt.xlabel(\"$E - E_F$ [eV]\")\n",
    "plt.ylabel(\"DOS [1/eV]\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Berry phase calculation\n",
    "\n",
    "To calculate the Berry phase we are going to perform a discretized integration of the Bloch states on a closed loop. We are going to calculate it around the $K$-point with a given radius. After having created the "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Number of discretizations\n",
    "N = 50\n",
    "# Circle radius in 1/Ang\n",
    "kR = 0.01\n",
    "# Normal vector (in units of reciprocal lattice vectors)\n",
    "normal = [0, 0, 1]\n",
    "# Origo (in units of reciprocal lattice vectors)\n",
    "origin = [1 / 3, 2 / 3, 0]\n",
    "circle = BrillouinZone.param_circle(H, N, kR, normal, origin)\n",
    "plt.plot(circle.k[:, 0], circle.k[:, 1])\n",
    "plt.xlabel(r\"$k_x$ [$b_x$]\")\n",
    "plt.ylabel(r\"$k_y$ [$b_y$]\")\n",
    "plt.gca().set_aspect(\"equal\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above plot shows a skewed circle because the $k$-points in the Brillouin zone object is stored in units of reciprocal lattice vectors. I.e. the circle is perfect in the reciprocal space. Note that the below Berry phase calculation ensures the loop is completed by also taking into account the first and last point.\n",
    "\n",
    "To confirm that the circle is *perfect* in reciprocal space, we convert the $k$-points and plot again. Note also that the radius of the circle is $0.01\\mathrm{Ang}^{-1}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k = circle.tocartesian(circle.k)\n",
    "plt.plot(k[:, 0], k[:, 1])\n",
    "plt.xlabel(r\"$k_x$ [1/Ang]\")\n",
    "plt.ylabel(r\"$k_y$ [1/Ang]\")\n",
    "plt.gca().set_aspect(\"equal\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we are ready to calculate the Berry phase. We calculate it for both graphene and the anti-symmetric graphene using only the first, second and both bands:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "circle.set_parent(H)\n",
    "print(\"Pristine graphene (0): {:.5f} rad\".format(electron.berry_phase(circle, sub=0)))\n",
    "print(\"Pristine graphene (1): {:.5f} rad\".format(electron.berry_phase(circle, sub=1)))\n",
    "print(\"Pristine graphene (:): {:.5f} rad\".format(electron.berry_phase(circle)))\n",
    "circle.set_parent(H_bp)\n",
    "print(\n",
    "    \"Anti-symmetric graphene (0): {:.5f} rad\".format(\n",
    "        electron.berry_phase(circle, sub=0)\n",
    "    )\n",
    ")\n",
    "print(\n",
    "    \"Anti-symmetric graphene (1): {:.5f} rad\".format(\n",
    "        electron.berry_phase(circle, sub=1)\n",
    "    )\n",
    ")\n",
    "print(\"Anti-symmetric graphene (:): {:.5f} rad\".format(electron.berry_phase(circle)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now plot the Berry phase as a function of integration radius with a somewhat constant discretization. In addition we calculate the Berry phase in the skewed circle in reciprocal space and perfectly circular in the units of the reciprocal lattice vectors. This enables a comparison of the integration path."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "kRs = np.linspace(0.001, 0.2, 70)\n",
    "dk = 0.0001\n",
    "bp = np.empty([4, len(kRs)])\n",
    "for i, kR in enumerate(kRs):\n",
    "    circle = BrillouinZone.param_circle(H_bp, dk, kR, normal, origin)\n",
    "    bp[0, i] = electron.berry_phase(circle, sub=0)\n",
    "    circle_other = BrillouinZone.param_circle(\n",
    "        utils.mathematics.fnorm(H_bp.rcell), dk, kR, normal, origin\n",
    "    )\n",
    "    circle.k[:, :] = circle_other.k[:, :]\n",
    "    bp[1, i] = electron.berry_phase(circle, sub=0)\n",
    "plt.plot(kRs, bp[0, :] / np.pi, label=r\"1/Ang\")\n",
    "plt.plot(kRs, bp[1, :] / np.pi, label=r\"$b_i$\")\n",
    "plt.legend()\n",
    "plt.xlabel(r\"Integration radius [1/Ang]\")\n",
    "plt.ylabel(r\"Berry phase [$\\phi/\\pi$]\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.9.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}