docs/tutorials/tutorial_siesta_2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "os.chdir('siesta_2')\n",
    "import numpy as np\n",
    "from sisl import *\n",
    "import sisl.viz\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Siesta ---  graphene\n",
    "\n",
    "This tutorial will describe a complete walk-through of a large fraction of the `sisl` functionalities that may be related to the [Siesta code](https://gitlab.com/siesta-project/siesta).\n",
    "\n",
    "Contrary to the $\\mathrm H_2\\mathrm O$ system this tutorial will emphasize the usefulness of performing bandstructures etc. directly in Python using sisl.\n",
    "\n",
    "## Creating the geometry\n",
    "\n",
    "Our system of interest will be the smallest graphene cell. Instead of defining the atomic positions, the carbon atoms and supercell for graphene, we use a default implementation of graphene in `sisl`. There are a small selection of the typical geometries, including graphene."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "graphene = geom.graphene(1.44)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us plot the geometry."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "graphene.plot(axes=\"xy\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we need to create the input fdf file for Siesta:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "open('RUN.fdf', 'w').write(\"\"\"%include STRUCT.fdf\n",
    "SystemLabel siesta_2\n",
    "PAO.BasisSize SZP\n",
    "MeshCutoff 250. Ry\n",
    "CDF.Save true\n",
    "CDF.Compress 9\n",
    "SaveHS true\n",
    "SaveRho true\n",
    "%block kgrid.MonkhorstPack\n",
    "  61  1 1 0.\n",
    "   1 61 1 0.\n",
    "   0  0 1 0.\n",
    "%endblock\n",
    "\"\"\")\n",
    "graphene.write('STRUCT.fdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creating the electronic structure\n",
    "\n",
    "Before proceeding, run Siesta to calculate the ground state electronic structure.\n",
    "\n",
    "After having completed the Siesta run we may read the Hamiltonian to manipulate and extract different information.\n",
    "After reading the Hamiltonian it is obvious that a great deal of new data has been associated with the Hamiltonian."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fdf = get_sile('RUN.fdf')\n",
    "H = fdf.read_hamiltonian()\n",
    "print(H)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculating DOS and PDOS\n",
    "\n",
    "When we are dealing with periodic structures (such as graphene) it is imperative to calculate the density of states in a simple and efficient manner. Below we will calculate the DOS for a variety of Monkhorst-Pack grids to check the convergence of the DOS (it shouldn't take more than a minute):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "E = np.linspace(-6, 4, 500)\n",
    "for nk in [21, 41, 61, 81]:\n",
    "    bz = MonkhorstPack(H, [nk, nk, 1])\n",
    "    plt.plot(E, bz.apply.average.eigenvalue(wrap=lambda ev: ev.DOS(E)), label='nk={}'.format(nk));\n",
    "plt.xlim(E[0], E[-1]); plt.ylim(0, None)\n",
    "plt.xlabel(r'$E - E_F$ [eV]')\n",
    "plt.ylabel(r'DOS [1/eV]')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the k-points are converged for approximately 61 k-points using the default smearing method. The default smearing method is the Gaussian smearing technique with $\\sigma=0.1\\mathrm{eV}$. Note that intrinsically the `MonkhorstPack` grid assumes time-reversal symmetry. I.e. $\\mathbf k \\equiv -\\mathbf k$.\n",
    "\n",
    "Now we may use the Monkhorst-Pack grid for 81 points to find the projected DOS for some of the orbitals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bz = MonkhorstPack(H, [81, 81, 1])\n",
    "idx_s = list()\n",
    "idx_pxy = list()\n",
    "idx_pz = list()\n",
    "for i, orb in enumerate(H.geometry.atoms[0]):\n",
    "    if orb.l == 0:\n",
    "        idx_s.append(i)\n",
    "    elif orb.l == 1 and (orb.m in [-1, 1]):\n",
    "        idx_pxy.append(i)\n",
    "    elif orb.l == 1 and orb.m == 0:\n",
    "        idx_pz.append(i)\n",
    "print('Orbital index of s: {}'.format(idx_s))\n",
    "print('Orbital index of p_x and p_y: {}'.format(idx_pxy))\n",
    "print('Orbital index of p_z: {}'.format(idx_pz))\n",
    "# Get all orbitals\n",
    "all_s = np.add.outer(H.geometry.a2o([0, 1]), idx_s).ravel()\n",
    "all_pxy = np.add.outer(H.geometry.a2o([0, 1]), idx_pxy).ravel()\n",
    "all_pz = np.add.outer(H.geometry.a2o([0, 1]), idx_pz).ravel()\n",
    "\n",
    "def wrap(es):\n",
    "    PDOS = es.PDOS(E)[0]\n",
    "    pdos_s = PDOS[all_s, :].sum(0)\n",
    "    pdos_pxy = PDOS[all_pxy, :].sum(0)\n",
    "    pdos_pz = PDOS[all_pz, :].sum(0)\n",
    "    return np.stack((pdos_s, pdos_pxy, pdos_pz))\n",
    "pDOS = bz.apply.average.eigenstate(wrap=wrap)\n",
    "plt.plot(E, pDOS[0, :], label='$s$');\n",
    "plt.plot(E, pDOS[1, :], label='$p_x+p_y$');\n",
    "plt.plot(E, pDOS[2, :], label=r'$p_z$');\n",
    "plt.xlim(E[0], E[-1]); plt.ylim(0, None)\n",
    "plt.xlabel(r'$E - E_F$ [eV]')\n",
    "plt.ylabel(r'DOS [1/eV]')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As seen the $p_z$ orbitals are responsible for the DOS in a broad range of energies around the Fermi-level. This is one reason for the tight-binding models success with respect to graphene.\n",
    "\n",
    "Another way to gain information is via the so-called *fat-bands* which basically is the PDOS scaling (no broadening) on each band for the quantities we are interested in. To plot the fat-bands we need the band-structure and a projection of each state onto the requested orbitals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "weight_s = list()\n",
    "weight_pxy = list()\n",
    "weight_pz = list()\n",
    "def wrap_fatbands(eigenstate):\n",
    "    # The eigenstate object has several features.\n",
    "    # For now we will simply calculate the <psi_i| S(k) | psi_i> weight for\n",
    "    # the orbitals we are interested in.\n",
    "    norm2 = eigenstate.norm2(sum=False)\n",
    "    weight_s.append(norm2[:, all_s].sum(-1))\n",
    "    weight_pxy.append(norm2[:, all_pxy].sum(-1))\n",
    "    weight_pz.append(norm2[:, all_pz].sum(-1))\n",
    "    return eigenstate.eig\n",
    "# Define the band-structure\n",
    "bz = BandStructure(H, [[0] * 3, [2./3, 1./3, 0], [0.5, 0.5, 0], [1] * 3], 400, \n",
    "                   name=[r'$\\Gamma$', r'$K$', r'$M$', r'$\\Gamma$'])\n",
    "\n",
    "# Calculate all eigenvalues\n",
    "eig = bz.apply.array.eigenstate(wrap=wrap_fatbands).T\n",
    "weight_s = np.array(weight_s).T\n",
    "weight_pxy = np.array(weight_pxy).T\n",
    "weight_pz = np.array(weight_pz).T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally we plot the fat-bands. Our plots split each of the orbitals contributions into their separate fat band. Thus one can see details on multiple orbital contributions for each band."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "linear_k, k_tick, k_label = bz.lineark(True)\n",
    "\n",
    "Emin, Emax = -21, 10\n",
    "# This is to determine the width of the fat-bands\n",
    "# The width of the fat-bands is dependent on the energy range and also on the variety\n",
    "# of contributions.\n",
    "dE = (Emax - Emin) / 20.\n",
    "plt.ylabel(r'$E-E_F$ [eV]')\n",
    "plt.xlim(linear_k[0], linear_k[-1]);\n",
    "plt.xticks(k_tick, k_label);\n",
    "plt.ylim(Emin, Emax);\n",
    "\n",
    "# Now plot the bands\n",
    "for i, e in enumerate(eig):\n",
    "    s = np.abs(weight_s[i, :] * dE)\n",
    "    pxy = np.abs(weight_pxy[i, :] * dE)\n",
    "    pz = np.abs(weight_pz[i, :] * dE)\n",
    "    plt.plot(linear_k, e, color='k'); # black-line (band-structure)\n",
    "    # Full fat-band\n",
    "    plt.fill_between(linear_k, e - dE, e + dE, color='k', alpha=0.1);\n",
    "    # pz\n",
    "    plt.fill_between(linear_k, e + (s + pxy), e + (s + pxy + pz), color='g', alpha=0.5);\n",
    "    plt.fill_between(linear_k, e - (s + pxy + pz), e - (s + pxy), color='g', alpha=0.5);\n",
    "    # pxy\n",
    "    plt.fill_between(linear_k, e + (s), e + (s + pxy), color='b', alpha=0.5);\n",
    "    plt.fill_between(linear_k, e - (s), e - (s + pxy), color='b', alpha=0.5);\n",
    "    # s\n",
    "    plt.fill_between(linear_k, e - (s), e + (s), color='r', alpha=0.5);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- all, gray\n",
    "- $p_z$, green\n",
    "- $s$, red\n",
    "- $p_x+p_y$, blue"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hamiltonian eigenstates\n",
    "\n",
    "At this point we have plotted the $k$-averaged DOS, PDOS. We have also plotted the fat-bands (and thus the band-structure). \n",
    "\n",
    "In addition to these things we can plot the real-space eigenstates. We first plot the $\\Gamma$-point for the first 2x2 unit-cell. This $k$-point has complete unit-cell periodicity and thus the plotted wavefunction should be fully periodic along all directions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "es = H.eigenstate()\n",
    "idx_valence = (es.eig > 0).nonzero()[0][0] - 1\n",
    "# Only select the valence band state\n",
    "es = es.sub(idx_valence)\n",
    "\n",
    "# Generate a grid encompassing a 2x2 graphene unit-cell\n",
    "g = Grid(0.2, sc=H.geometry.sc.tile(2, 0).tile(2, 1))\n",
    "# Calculate the real-space wavefunctions\n",
    "es.wavefunction(g)\n",
    "\n",
    "# Extract the wavefunction a few Ang above the graphene plane\n",
    "# To do this we need to find the index of the corresponding z-plane.\n",
    "# The Grid.index method is useful in this regard.\n",
    "xyz = H.geometry.xyz[0, :].copy()\n",
    "xyz[2] += 1.\n",
    "z_idx = g.index(xyz, axis=2)\n",
    "x, y = np.mgrid[:g.shape[0], :g.shape[1]]\n",
    "x, y = x * g.dcell[0, 0] + y * g.dcell[1, 0], x * g.dcell[0, 1] + y * g.dcell[1, 1]\n",
    "plt.contourf(x, y, g.grid[:, :, z_idx]);\n",
    "xyz = H.geometry.tile(2, 0).tile(2, 1).xyz\n",
    "plt.scatter(xyz[:, 0], xyz[:, 1], 20, c='k');\n",
    "plt.xlabel(r'$x$ [Ang]');\n",
    "plt.ylabel(r'$y$ [Ang]');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will now try and plot the real-space wavefunction for a finite $k$. By choosing the $[1/2, 0, 0]$ point we know it must have a periodicity of 2 along the first lattice vector (this lattice vector is pointing right-up), and full periodicity along the second lattice vector. Since we have a finite $k$ the grid data-type *must* be complex because the eigenstates have complex components. And thus we will plot both the real and imaginary part."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_, axs = plt.subplots(1, 2, figsize=(16, 5))\n",
    "es = H.eigenstate([1./2, 0, 0])\n",
    "idx_valence = (es.eig > 0).nonzero()[0][0] - 1\n",
    "es = es.sub(idx_valence)\n",
    "g = Grid(0.2, dtype=np.complex128, sc=H.geometry.sc.tile(4, 0).tile(4, 1))\n",
    "es.wavefunction(g)\n",
    "x, y = np.mgrid[:g.shape[0], :g.shape[1]]\n",
    "x, y = x * g.dcell[0, 0] + y * g.dcell[1, 0], x * g.dcell[0, 1] + y * g.dcell[1, 1]\n",
    "axs[0].contourf(x, y, g.grid[:, :, z_idx].real);\n",
    "axs[1].contourf(x, y, g.grid[:, :, z_idx].imag);\n",
    "xyz = H.geometry.tile(4, 0).tile(4, 1).xyz\n",
    "for ax in axs:\n",
    "    ax.scatter(xyz[:, 0], xyz[:, 1], 20, c='k');\n",
    "    ax.set_xlabel(r'$x$ [Ang]');\n",
    "    ax.set_ylabel(r'$y$ [Ang]');"
   ]
  }
 ],
 "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
}