3_Density_functional_theory.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "divided-institution",
   "metadata": {},
   "source": [
    "# Density-functional theory\n",
    "\n",
    "In this notebook we will dive into density-functional theory (DFT).\n",
    "As I will describe simulating electrons in periodic systems with DFT\n",
    "is a common setting where the periodic one-particle Hamiltonians $H$\n",
    "described in the previous notebook arise."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "annual-alcohol",
   "metadata": {},
   "source": [
    "## Motivation and quantum-mechanical setting\n",
    "\n",
    "- The goal of electronic structure is to describe the behaviour of electrons in matter.\n",
    "- This allows to understand the chemical and physical properties of materials.\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "- For accurately predicting properties from first principles a quantum-mechanical description of electrons is required.\n",
    "- The **states** describing the motion and behaviour of electrons are complex-valued\n",
    "  functions in an infinite-dimensional function space, which is a Hilbert space.\n",
    "- The **operators** are mappings defined on this Hilbert space.\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "- The motion of electrons is highly correlated. Pictorially speaking the motion of one electron depends on the motion of all others.\n",
    "- This makes an accurate modelling of the electronic structure rather involved.\n",
    "- In this workshop we focus on **density-functional theory** (DFT).\n",
    "  While this approximate model has a number of limitations,\n",
    "  it has become one of the most widely used practical methods due to its balance \n",
    "  between cost and accuracy. \n",
    "- Moreover problems of the form of DFT are usually the first step in modelling\n",
    "  the electronic structure, even when a more accurate methods will eventually be\n",
    "  targeted. So understanding how DFT works represents a good step in either case.\n",
    "  \n",
    "**More details**\n",
    "- My [2020 Juliacon talk](https://www.youtube.com/watch?v=-RomkxjlIcQ) (motivation for electronic-structure theory and DFT)\n",
    "- Chapter 2 of my [PhD thesis](https://michael-herbst.com/publications/2018.05_phd_corrected.pdf) (Mathematical formulation of quantum mechanics)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "civil-science",
   "metadata": {},
   "source": [
    "## Formulation of the DFT problem\n",
    "\n",
    "- In DFT an electronic state is described by a set of orthonormal **orbitals**\n",
    "  $\\{\\psi_i\\}$ ($i = 1, 2, \\ldots N$).\n",
    "- Each of these orbitals is a function from $\\mathbb{R}^3$ to $\\mathbb{C}$.\n",
    "- To a good approximation one can think of one orbital $\\psi_i$ as the description\n",
    "  of the behaviour of one electron pair.\n",
    "- If there are $2 N$ electrons in the system, we therefore need $N$ orbitals."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "noted-render",
   "metadata": {},
   "source": [
    "- Of largest physical interest is the ground state of a system, that is the state\n",
    "  which features the lowest energy.\n",
    "- At the level of DFT this translates to the minimisation problem\n",
    "  $$ \\displaystyle \\text{min}_{\\{\\psi_i\\}} \\mathcal{E}_\\text{DFT}(\\{\\psi_i\\}) $$\n",
    "  where we seek the lowest-energy orbitals given a DFT energy functional.\n",
    "\n",
    "\n",
    "- The details of the energy functional depend on the chosen DFT model.\n",
    "  The common structure is as follows:\n",
    "\n",
    "    $$ \\mathcal{E}_\\text{DFT}(\\{\\psi_i\\}) = \\sum_{i=1}^N 2 \\int \\psi_i^\\ast \\left(-\\frac12 \\Delta\\right) \\psi_i\n",
    "    + \\int V_\\text{ext} \\rho + \\int V_\\text{H}[\\rho] \\rho + \\int V_\\text{xc}[\\rho] \\rho + E_\\text{nuclear}$$\n",
    "    \n",
    "  with\n",
    "     * the electron **density** $\\rho = \\sum_i^N |\\psi_i|^2$\n",
    "       being directly dependent on *all* orbitals.\n",
    "     * $\\sum_i 2 \\int \\psi_i^\\ast \\left(-\\frac12 \\Delta\\right) \\psi_i$ describing\n",
    "       the **kinetic** energy of the electrons\n",
    "       (i.e. the term we discussed in the previous notebook)\n",
    "     * $\\int V_\\text{ext} \\rho$ being the **external** potential energy,\n",
    "       i.e. the electron-nuclear interaction.\n",
    "     * $\\int V_\\text{H}[\\rho] \\rho$ being the **Hartree** energy,\n",
    "       i.e. the classical Coulombic electron-electron repulsion\n",
    "     * $\\int V_\\text{xc}[\\rho] \\rho$ being the **exchange-correlation**\n",
    "       (or XC) energy, which describes the electron-electron interaction beyond\n",
    "       the classical model.\n",
    "     * $E_\\text{nuclear}$: A fixed energy offset due to the repulsion of the nuclei\n",
    "       from each other.\n",
    "\n",
    "**Note:** The Hartree potential $V_\\text{H}$ and the XC potential $V_\\text{xc}$ themselves depend on the density $\\rho$. Therefore $\\mathcal{E}_\\text{DFT}$ is non-linear in the orbitals!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "tested-improvement",
   "metadata": {},
   "source": [
    "For our numerical experiments we will use the [density-functional toolkit](https://dftk.org) (DFTK).\n",
    "In DFTK such a DFT model can be conveniently assembled by collecting the respective energy terms:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "written-thesaurus",
   "metadata": {},
   "outputs": [],
   "source": [
    "using DFTK\n",
    "\n",
    "# Build up the LDA model by assembling the energy terms\n",
    "terms_LDA = [Kinetic(),           # -1/2 Δ\n",
    "             AtomicLocal(),       # part of Vₑₓₜ\n",
    "             AtomicNonlocal(),    # part of Vₑₓₜ\n",
    "             PspCorrection(),     # part of Vₑₓₜ\n",
    "             Hartree(),           # Vₕ\n",
    "             Xc(:lda_xc_teter93), # Vxc\n",
    "             Ewald()]             # part of Eₙᵤᵪₗₑₐᵣ"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "urban-effects",
   "metadata": {},
   "source": [
    "DFTK has many more terms to offer (e.g. external fields, magnetism, etc.). See [the source code](https://dftk.org/tree/master/src/terms) for the full list."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "competitive-smell",
   "metadata": {},
   "source": [
    "## Solving DFT: The self-consistent field problem\n",
    "\n",
    "The stationarity conditions (Euler-Lagrange equations) corresponding to the DFT minimisation problem are\n",
    "\n",
    "$$\\left\\{\\begin{aligned}\n",
    "    &\\forall 1 \\leq i \\leq N: &\\left(-\\frac12 \\Delta + V[\\rho] \\right) \\psi_i &= \\varepsilon_i \\psi_i\n",
    "    \\qquad \\text{(with $\\varepsilon_1, \\varepsilon_2, \\ldots, \\varepsilon_N$ smallest)} \\\\\n",
    "    &\\forall 1 \\leq i, j \\leq N: &\\int \\psi_i^\\ast \\psi_j &= \\delta_{ij} \\\\\n",
    "    &&V[\\rho] &= V_\\text{ext} + V_\\text{H}[\\rho] + V_\\text{xc}[\\rho] \\\\\n",
    "    &&\\rho &= \\sum_{i=1}^N 2 |\\psi_i|^2\n",
    "\\end{aligned}\\right.$$\n",
    "\n",
    "in which $H[\\rho] = -\\frac12 \\Delta + V[\\rho]$ is the **DFT Hamiltonian**. These equations are called the **self-consistent field** (SCF) problem.\n",
    "\n",
    "- The first and second line are basically an eigenvalue problem\n",
    "  with the usual orthogonality between the resulting eigenfunctions.\n",
    "- Since $V$ depends on $\\rho$ (3rd line), which itself depends on $\\psi_i$ (last line),\n",
    "  solving the SCF equations implies solving a *non-linear* eigenvalue problem.\n",
    "- In this sense the DFT Hamiltonian is a generalisation of the diagonalisation problems\n",
    "  we discussed in the previous notebook."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "amino-grounds",
   "metadata": {},
   "source": [
    "To allow for writing the SCF problem more compact\n",
    "we define the **potential-to-density map**\n",
    "$$\n",
    "D(V) = \\sum_{i=1}^N |\\psi_i|^2 \\qquad \\text{$\\psi_i$ are the $N$ lowest eigenvectors of $-\\frac12 \\Delta + V$}.\n",
    "$$\n",
    "With it the SCF problem can be written as\n",
    "$$ \\rho = D(V(\\rho)). $$\n",
    "Clearly an SCF is nothing else but a **fixed-point problem**."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "changing-maine",
   "metadata": {},
   "source": [
    "**Note:** In the picture I sketch in this notebook I make a few assumptions for simplicity (e.g. no unpaired electrons, no spin, only integer occupation, existence of a gap, etc.). Adding in these extra complications does not change much of the bottom line, so I've decided to leave these aspects out for clarity. You can find more details (at various levels) in the [references given in the DFTK documentation](https://juliamolsim.github.io/DFTK.jl/stable/guide/density_functional_theory/)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "metallic-newspaper",
   "metadata": {},
   "source": [
    "## Modelling silicon\n",
    "\n",
    "Puh ... that was quite a lot of equations. Let's see some numbers!\n",
    "\n",
    "First we will use the LDA model defined above to model silicon,\n",
    "a nicely periodic lattice.\n",
    "\n",
    "| <img src=\"img/Silicon_crystal.jpg\" width=180 height=180 /> |\n",
    "| --------------- |\n",
    "| silicon crystal |"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "celtic-colon",
   "metadata": {},
   "outputs": [],
   "source": [
    "using DFTK\n",
    "\n",
    "# Look up the structure of silicon in a book\n",
    "# and define the lattice and the atomic positions\n",
    "a = 10.26\n",
    "lattice = a / 2 * [[0 1 1.];\n",
    "                   [1 0 1.];\n",
    "                   [1 1 0.]]\n",
    "Si = ElementPsp(:Si, psp=load_psp(\"hgh/lda/si-q4\"))\n",
    "atoms = [Si => [ones(3)/8, -ones(3)/8]]\n",
    "\n",
    "# Setup a model in DFTK from the terms we defined above\n",
    "model = Model(lattice; atoms, terms=terms_LDA)\n",
    "# Note: For common models DFTK has shortcuts. In this case\n",
    "#            model = model_LDA(lattice, atoms)\n",
    "#       does exactly the same\n",
    "\n",
    "# Discretise the problem in a plane-wave basis\n",
    "# We use a mesh of 2x2x2 k-Points in the Brillouin zone\n",
    "# and an energy cutoff of 13 Hartree\n",
    "basis  = PlaneWaveBasis(model; Ecut=15, kgrid=[2, 2, 2])\n",
    "\n",
    "# Run the SCF and catch the result:\n",
    "scfres = self_consistent_field(basis);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "integrated-alpha",
   "metadata": {},
   "source": [
    "Ok so this did something, but what?\n",
    "\n",
    "- First DFTK inspects the molecule to guess an initial density $\\rho_0$.\n",
    "- Then we solve $ \\rho = D(V(\\rho)) $ by computing\n",
    "  $D(V(\\rho_n))$ for a sequence of iterates $\\rho_n$ until input and output\n",
    "  are close enough, i.e. the residual\n",
    "  $$ R(\\rho_n) = D(V(\\rho_n)) - \\rho_n$$\n",
    "  is small. Then DFTK flags convergence.\n",
    "- We will discuss these algorithms with more details in the next notebook.\n",
    "\n",
    "In the end we obtained the ground-state energy of silicon as"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "mighty-istanbul",
   "metadata": {},
   "outputs": [],
   "source": [
    "scfres.energies"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "single-avenue",
   "metadata": {},
   "source": [
    "We can also look at the density:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "searching-condition",
   "metadata": {},
   "outputs": [],
   "source": [
    "using Plots\n",
    "heatmap(scfres.ρ[:, :, 5], c=:Blues)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "genuine-aging",
   "metadata": {},
   "source": [
    "... or access the eigenpairs of the final Hamiltonian diagonalisation\n",
    "in `scfres.ψ` and `scfres.eigenvalues`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "going-thompson",
   "metadata": {},
   "outputs": [],
   "source": [
    "scfres.eigenvalues"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fifty-simple",
   "metadata": {},
   "source": [
    "It might come to a surprise for your why this does not return a flat list,\n",
    "but in fact a list of two lists. This will be clarified in the following\n",
    "\n",
    "\n",
    "**Remark:** Since silicon (and in fact all problems we will consider) are periodic,\n",
    "the DFT Hamiltonian $H$ is also periodic.\n",
    "Therefore we are in the periodic setting discussed in the previous notebook\n",
    "and one is able to apply the Bloch transform to $H$.\n",
    "This allows to reformulate the SCF problem as \n",
    "\n",
    "$$\\left\\{\\begin{aligned}\n",
    "    &\\forall k\\in\\text{BZ}, \\ \\  1 \\leq i \\leq N: &\\left(\\frac12 (-i \\nabla + k)^2 + V[\\rho] \\right) \\psi_{ki} &= \\varepsilon_{ki} \\psi_{ki}\\\\\n",
    "    &\\forall k\\in\\text{BZ}, \\ \\ 1 \\leq i, j \\leq N: &\\int \\psi_{ki}^\\ast \\psi_{kj} &= \\delta_{ij} \\\\\n",
    "    &&V[\\rho] &= V_\\text{ext} + V_\\text{H}[\\rho] + V_\\text{xc}[\\rho] \\\\\n",
    "    &&\\rho &= \\sum_{k\\in\\text{BZ}} \\sum_{i=1}^N 2 |\\psi_{ki}|^2\n",
    "\\end{aligned}\\right.$$\n",
    "\n",
    "In other words we can diagonalisate $H$ $k$-Point by $k$-Point\n",
    "and only at the level of computing the density we need to\n",
    "sum over $k$-Points.\n",
    "\n",
    "In DFTK most quantities in the `scfres` are therefore stored\n",
    "as a list over $k$-Points. Therefore\n",
    "```julia\n",
    "    scfres.eigenvalues[1]\n",
    "```\n",
    "gives access to the $\\varepsilon_{ki}$ for the first $k$-Point."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "allied-tract",
   "metadata": {},
   "source": [
    "In our problem we used `kgrid=[2,2,2]`. One would therefore expect 8 $k$-Points\n",
    "to be present in the discretisation. As it turns out there are only"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "delayed-installation",
   "metadata": {},
   "outputs": [],
   "source": [
    "length(basis.kpoints)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "charming-portugal",
   "metadata": {},
   "source": [
    "The reason for this is that the symmetries of the problem allow to make further reductions in effort. If we disable symmetries, we get the expected 8 $k$-Points:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "manual-jordan",
   "metadata": {},
   "outputs": [],
   "source": [
    "model_nosym = Model(lattice; atoms, terms=terms_LDA, symmetries=false)\n",
    "length(PlaneWaveBasis(model_nosym; Ecut=15, kgrid=(2, 2, 2)).kpoints)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "perceived-undergraduate",
   "metadata": {},
   "source": [
    "To conclude this section, let us plot the band structure of silicon,\n",
    "i.e. the relationship between the eigenvalues of $H$\n",
    "and the position $k$ inside the Brillouin zone.\n",
    "(For more details see the previous notebook)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "single-poultry",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_bandstructure(scfres)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "scheduled-finance",
   "metadata": {},
   "source": [
    "## Modelling aluminum\n",
    "\n",
    "**Exercise:** Try running a simulation of aluminium yourself.\n",
    "\n",
    "For aluminum a possible structural setup is\n",
    "```julia\n",
    "a = 7.65339\n",
    "lattice = diagm(fill(a, 3))\n",
    "Al = ElementPsp(:Al, psp=load_psp(\"hgh/lda/al-q3\"))\n",
    "atoms = [Al => [[0.0, 0.0, 0.0], [0.0, 0.5, 0.5], [0.5, 0.0, 0.5], [0.5, 0.5, 0.0]]]\n",
    "```\n",
    "\n",
    "First use the same LDA model as before, defined using the DFTK shortcut\n",
    "```julia\n",
    "model = model_LDA(lattice, atoms)\n",
    "```\n",
    "and use `Ecut=15` and `kgrid=[2, 2, 2]` to discretise the problem.\n",
    "\n",
    "When running the SCF like this you should observe either slow or no convergence.\n",
    "This points to the fact that modelling aluminium is clearly a little different\n",
    "than silicon.\n",
    "In fact the issue is due to aluminium being a metal,\n",
    "while silicon is an insulator / semiconductor.\n",
    "More on this below.\n",
    "\n",
    "To make progress and achieve SCF convergence add the parameter\n",
    "```julia\n",
    "model = model_LDA(lattice, atoms, temperature=1e-3)\n",
    "```\n",
    "This is called **smearing** and is a numerical trick\n",
    "to improve the convergence for quantities which involve\n",
    "sums over $k$-Points in metals.\n",
    "As the density $\\rho$ involves such a sum, the SCF typically fails\n",
    "to converge in metals without smearing.\n",
    "\n",
    "With this setup it the SCF should succeed. Use the obtained scfres to compute the band structure of aluminium with `plot_bandstructure(scfres)`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bearing-donna",
   "metadata": {},
   "outputs": [],
   "source": [
    "using DFTK\n",
    "using LinearAlgebra\n",
    "\n",
    "a = 7.65339\n",
    "lattice = diagm(fill(a, 3))\n",
    "Al = ElementPsp(:Al, psp=load_psp(\"hgh/lda/al-q3\"))\n",
    "atoms = [Al => [[0.0, 0.0, 0.0], [0.0, 0.5, 0.5], [0.5, 0.0, 0.5], [0.5, 0.5, 0.0]]]\n",
    "\n",
    "# your solution here "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "joint-purple",
   "metadata": {},
   "source": [
    "## Contrasting aluminium and silicon"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "exceptional-license",
   "metadata": {},
   "source": [
    "With these two calculations we have computed two band structures,\n",
    "one for a metal (aluminium) and one for a semiconductor (silicon).\n",
    "\n",
    "- All details of the band structure do not matter for our purposes.\n",
    "- The important difference is that silicon has a so-called **band gap**,\n",
    "  i.e. a region near the Fermi level (which is at $0$ in the shown band plots)\n",
    "  where there are no eigenvalues, i.e. no blue bands.\n",
    "  In contrast aluminium does not have this gap.\n",
    "- This makes aluminium a metal and silicon an insulator respectively semiconductor.\n",
    "\n",
    "For the numerical treatment the absence of a band gap has strong implications:\n",
    "- At the beginning of the SCF the DFT Hamiltonian $H[\\rho_n]$\n",
    "  at the current density $\\rho_n$ is only a rough approximation to the\n",
    "  converged Hamiltonian $H[\\rho_\\ast]$.\n",
    "- This means that the order of the $\\psi_{ki}$ near the Fermi level,\n",
    "  i.e. near $i = N$, may be easily swapped as many $\\psi_{ki}$\n",
    "  of similar eigenvalues $\\varepsilon_i$ are close.\n",
    "- From one iteration to the next the $\\psi_{ki}$ which are summed\n",
    "  over in the density expression\n",
    "  $$ \\rho = \\sum_{k\\in\\text{BZ}} \\sum_{i=1}^N 2 |\\psi_{ki}|^2 $$\n",
    "  may thus change (as $N$ becomes $N+1$, thus not summed over, but $N+1$ becomes $N$).\n",
    "- Without going into details the smearing temperature reduces this\n",
    "  problem by not enforcing a hard cut between the functions summed over in the density,\n",
    "  but employing a smoother cutoff function.\n",
    "\n",
    "- As becomes apparent already in these first few experiments\n",
    "  the numerics of an SCF and the physics of the system are closely linked.\n",
    "\n",
    "\n",
    "#### Takeaways\n",
    "- Insulators / semiconductors have a band gap\n",
    "- Metals have no band gap\n",
    "- Without a band gap a smearing is needed to aid convergence.\n",
    "- SCF convergence properties for metals and insulators differ quite drastically (more on this in the next notebooks)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.6.1",
   "language": "julia",
   "name": "julia-1.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}