{
 "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": [
    "# From a geometry to a Hamiltonian\n",
    "\n",
    "This tutorial sets up tight-binding Hamiltonians from a `Geometry`.\n",
    "It explains the supercell, shows how to set matrix elements by hand,\n",
    "and then how to set them from the distances between atoms with `Hamiltonian.construct`.\n",
    "\n",
    "It assumes that you know how to create a geometry, see the [geometry quickstart](../quickstart/intro_tutorials/01_geometry.ipynb).\n",
    "\n",
    "## Supercells\n",
    "\n",
    "A Hamiltonian can only couple orbitals within the *supercell* of the geometry.\n",
    "The number of supercells, `nsc`, is a property of the `Lattice`.\n",
    "It is the number of cells along each lattice vector: the unit cell, plus the neighboring cells on both sides.\n",
    "Thus `nsc` is *always* an odd number. A molecule has `nsc=[1, 1, 1]`, i.e. only the unit cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "geometry = si.Geometry([[0, 0, 0]], lattice=si.Lattice(5, nsc=[3, 3, 3]))\n",
    "print(geometry.nsc)\n",
    "geometry.set_nsc([3, 1, 5])\n",
    "print(geometry.nsc)\n",
    "print(\"Number of supercells:\", geometry.n_s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The final geometry connects to the closest neighboring cell on each side along the first lattice vector (`1 + 2 == 3`),\n",
    "and to the 2 closest neighboring cells on each side along the third lattice vector (`1 + 2 + 2 == 5`).\n",
    "Methods that look for neighbors, such as `Geometry.close`, only find atoms in these supercells.\n",
    "\n",
    "The figure below shows a square 2D lattice with one atom in the unit cell, and `nsc=[5, 3, 1]`.\n",
    "A coupling within the supercell can be defined (full green arrow).\n",
    "A coupling that reaches farther than the supercell cannot be defined, and is thus *always* zero (dashed red arrow).\n",
    "\n",
    "![Square lattice with a 5x3 supercell](figures/04_supercell.png)\n",
    "\n",
    "Even though the geometry is 2D, sisl *requires* the third lattice vector.\n",
    "The non-periodic direction must have `nsc` equal to 1.\n",
    "Below we create the square lattice with nearest neighbor couplings, i.e. `nsc=[3, 3, 1]`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "square = si.Geometry([[0.5, 0.5, 0]], lattice=si.Lattice([1, 1, 10], nsc=[3, 3, 1]))\n",
    "square"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The geometries from `sisl.geom` have the `nsc` for nearest neighbor couplings.\n",
    "For graphene, the green and blue arrows show the couplings of the two atoms in the unit cell:\n",
    "\n",
    "![Nearest neighbor couplings in graphene](figures/04_graphene_couplings.png)\n",
    "\n",
    "`Geometry.find_nsc` calculates the smallest `nsc` for an interaction range.\n",
    "Even the 3rd nearest neighbors of graphene (2.84 Å) are in the closest neighboring cells:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "graphene = si.geom.graphene()\n",
    "print(graphene.nsc)\n",
    "print(graphene.find_nsc(R=2.9))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use `Geometry.set_nsc` to apply the result.\n",
    "\n",
    "## Setting matrix elements by hand\n",
    "\n",
    "A `Hamiltonian` takes the orbitals and the lattice from a `Geometry`.\n",
    "In sisl all Hamiltonians are sparse matrices, and you assign matrix elements as in a regular matrix.\n",
    "A third index selects the supercell of the second orbital, as an offset in units of the lattice vectors.\n",
    "\n",
    "Below we set the on-site energy of the square lattice to $-4$ eV and the 4 nearest neighbor couplings to $1$ eV:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "H = si.Hamiltonian(square)\n",
    "H[0, 0] = -4\n",
    "H[0, 0, (1, 0)] = 1\n",
    "H[0, 0, (-1, 0)] = 1\n",
    "H[0, 0, (0, 1)] = 1\n",
    "H[0, 0, (0, -1)] = 1\n",
    "H"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`nnz=5` is the number of non-zero matrix elements.\n",
    "sisl does not enforce that the Hamiltonian is Hermitian. *That is the responsibility of the user*.\n",
    "\n",
    "`Hamiltonian.eigh` calculates the eigenvalues at a $\\mathbf k$ point, in units of the reciprocal lattice vectors.\n",
    "The default is the $\\Gamma$ point:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(H.eigh())\n",
    "print(H.eigh(k=[0.5, 0.5, 0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The nearest neighbor model of graphene has zero on-site energy and a hopping parameter of $-2.7$ eV.\n",
    "Each atom couples to 3 neighbors (see the figure above):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "H = si.Hamiltonian(graphene)\n",
    "H[0, 1] = -2.7\n",
    "H[0, 1, (-1, 0)] = -2.7\n",
    "H[0, 1, (0, -1)] = -2.7\n",
    "H[1, 0] = -2.7\n",
    "H[1, 0, (1, 0)] = -2.7\n",
    "H[1, 0, (0, 1)] = -2.7\n",
    "print(\"Gamma:\", H.eigh())\n",
    "print(\"K:\", H.eigh(k=[2 / 3, 1 / 3, 0]).round(8))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The two bands are at $\\pm 3 \\times 2.7$ eV at $\\Gamma$, and touch at the $K$ point (the Dirac point).\n",
    "\n",
    "### Band structure\n",
    "\n",
    "`BandStructure` creates a path through the Brillouin zone.\n",
    "The integer is the total number of $k$ points on the path.\n",
    "`apply.array` returns the eigenvalues at all $k$ points as one array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "band = si.BandStructure(\n",
    "    H,\n",
    "    [[0, 0, 0], [2 / 3, 1 / 3, 0], [0.5, 0.5, 0], [0, 0, 0]],\n",
    "    300,\n",
    "    names=[r\"$\\Gamma$\", \"K\", \"M\", r\"$\\Gamma$\"],\n",
    ")\n",
    "eigs = band.apply.array.eigh()\n",
    "\n",
    "plt.plot(band.lineark(), eigs)\n",
    "plt.xticks(*band.lineartick())\n",
    "plt.ylabel(\"Eigenvalue [eV]\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setting matrix elements from distances\n",
    "\n",
    "Setting each matrix element by hand does not scale to large geometries.\n",
    "In place, set the elements from the distances between atoms.\n",
    "For this, sisl must know the range of the orbitals.\n",
    "Here the atom is a hydrogen atom with one orbital with a range of 1 Å:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hydrogen = si.Atom(1, R=1.0)\n",
    "square = si.Geometry(\n",
    "    [[0.5, 0.5, 0]], hydrogen, lattice=si.Lattice([1, 1, 10], nsc=[3, 3, 1])\n",
    ")\n",
    "print(\"Maximum orbital range:\", square.maxR())\n",
    "square.close(0, R=[0.1, 1.1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`Geometry.close` returns the indices of the atoms within a distance.\n",
    "With a list of radii, the first array is the atoms in the sphere with the first radius\n",
    "(here only the atom itself), and each following array is the atoms in the shell between\n",
    "the radius and the previous radius (here the 4 nearest neighbors).\n",
    "Indices larger than `na` are atoms in the neighboring cells.\n",
    "\n",
    "With `close` we set the parameters for all atoms in a loop.\n",
    "All atoms at the same distance get the same parameter, thus the Hamiltonian is Hermitian:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "H = si.Hamiltonian(square)\n",
    "for ia in square:\n",
    "    idx = square.close(ia, R=[0.1, 1.1])\n",
    "    H[ia, idx[0]] = -4.0\n",
    "    H[ia, idx[1]] = 1.0\n",
    "H"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using `construct`\n",
    "\n",
    "The loop above is slow for large geometries, because each call to `close` searches the full geometry.\n",
    "`Hamiltonian.construct` does the same, but it searches only small blocks of the geometry at a time.\n",
    "It is the recommended method.\n",
    "For atoms with one orbital, pass the radii and the parameter for each shell:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "H = si.Hamiltonian(square)\n",
    "H.construct([[0.1, 1.1], [-4.0, 1.0]])\n",
    "print(H.eigh(k=[0.5, 0.5, 0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For atoms with more than one orbital, pass a function.\n",
    "It gets the Hamiltonian, the atom, the block of atoms to search, and their coordinates:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def func(H, ia, atoms, atoms_xyz=None):\n",
    "    idx = H.geometry.close(ia, R=[0.1, 1.1], atoms=atoms, atoms_xyz=atoms_xyz)\n",
    "    H[ia, idx[0]] = -4.0\n",
    "    H[ia, idx[1]] = 1.0\n",
    "\n",
    "\n",
    "H = si.Hamiltonian(square)\n",
    "H.construct(func)\n",
    "print(H.eigh(k=[0.5, 0.5, 0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For very large systems, construct the Hamiltonian of a small cell and then use\n",
    "`Hamiltonian.tile` or `Hamiltonian.repeat`. This is much faster than `construct` on the large geometry:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "H.tile(100, 0).tile(100, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Next steps\n",
    "\n",
    "The [electronic structure walk-through](tutorial_es_1.ipynb) uses these Hamiltonians to calculate\n",
    "eigenstates, density of states, band structures and real-space wavefunctions."
   ]
  }
 ],
 "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
}