6_Floating_point_error.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "polyphonic-disabled",
   "metadata": {},
   "source": [
    "# Tracking the floating-point error\n",
    "\n",
    "In the previous notebooks we were concerned with the convergence of the SCF procedure itself.\n",
    "Now we change gears a little and take a look at analysing the error within an obtained SCF solution.\n",
    "\n",
    "There are multiple sources of error in DFT (see [this paper](https://doi.org/10.1039/D0FD00048E) for details):\n",
    "1. An obvious one is the DFT model itself\n",
    "1. The error due to the discretisation (`Ecut` and `kgrid`)\n",
    "1. The error due to aborting the SCF before getting the residual to zero (`tol`)\n",
    "1. The error due to using a finite-precision floating-point arithmetic, the **floating-point error**\n",
    "\n",
    "In this notebook we will focus on the floating-point error and to try and understand the floating-point error in a converged SCF solution. Or in other words:\n",
    "$$ \\text{How many digits of my solution are trustworthy?} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "competent-census",
   "metadata": {},
   "source": [
    "## Solving DFT with different floating-point types\n",
    "\n",
    "One approximate way to get an idea of the floating-point error\n",
    "in a procedure is to solve the problem using higher precision\n",
    "floating-point types and compare the matching digits.\n",
    "\n",
    "To make things a little more interesting we will use 32bit floating-point arithmetic:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "radio-association",
   "metadata": {},
   "outputs": [],
   "source": [
    "using DFTK\n",
    "\n",
    "# Silicon lattice\n",
    "a = 10.26\n",
    "lattice = a / 2 .* [[0 1 1.]; [1 0 1.]; [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",
    "# DFTK will use the floating-point type used to represent the lattice\n",
    "# to deduce the floating-point type for the computation\n",
    "model = model_DFT(Array{Float32}(lattice), atoms, [:lda_x, :lda_c_vwn])\n",
    "basis = PlaneWaveBasis(model; Ecut=7, kgrid=[4, 4, 4], fft_size=(16, 16, 16))\n",
    "scfres = self_consistent_field(basis, tol=1e-10)\n",
    "\n",
    "results = Dict{DataType, Any}(\n",
    "    Float32 => scfres.energies.total\n",
    ")\n",
    "\n",
    "@show results[Float32];"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "direct-idaho",
   "metadata": {},
   "source": [
    "**Exercise:** To get a rough idea how many of these energy digits can be trusted,\n",
    "re-run the computation using `Float64` and `Double64` (a floating point type offering even higher accuracy that 64bits). For the `Double64` case also converge the SCF to higher precision (e.g. `tol=1e-16`).\n",
    "How many digits of the energy can be trusted at `Float32` and `Float64` level?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "traditional-round",
   "metadata": {},
   "outputs": [],
   "source": [
    "using GenericLinearAlgebra  # Needed to enable generic fallbacks in DFTK\n",
    "using DoubleFloats  # Defines Double64\n",
    "\n",
    "# Your solution here\n",
    "\n",
    "results[Float64]  = zero(Float64)\n",
    "results[Double64] = zero(Double64)\n",
    "\n",
    "@show results"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "iraqi-measurement",
   "metadata": {},
   "source": [
    "## Using IntervalArithmetic to get a solid answer\n",
    "\n",
    "Next we want to use the [`IntervalArithmetic`](https://github.com/JuliaIntervals/IntervalArithmetic.jl) package to rigorously track the error of the `Float32` computation.\n",
    "\n",
    "The idea of interval arithmetic is to represent a number by two floats in form of an interval $[a, b]$. This interval is constructed in a way that the exact number is always inside the interval. Moreover this representation is tracked through the full flow of the program. I.e. all operations (addition, multiplication, more involved calls) are performed on both numbers $a$ and $b$ simultanously, but making sure that operations on $a$ always round down and those on $b$ always round up.\n",
    "The result is thus again obtained as an interval,\n",
    "which is moreover guaranteed to bracket the exact answer. As a result: The tighter the interval, the smaller the floating-point error.\n",
    "\n",
    "To apply this to `DFTK` we first do the `Float32` computation as normal:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "confident-score",
   "metadata": {},
   "outputs": [],
   "source": [
    "using DFTK\n",
    "\n",
    "# Silicon lattice\n",
    "a = 10.26\n",
    "lattice = a / 2 .* [[0 1 1.]; [1 0 1.]; [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",
    "model = model_DFT(Array{Float32}(lattice), atoms, [:lda_x, :lda_c_vwn], symmetries=false)\n",
    "basis = PlaneWaveBasis(model; Ecut=7, kgrid=[4, 4, 4], fft_size=(16, 16, 16))\n",
    "scfres = self_consistent_field(basis, tol=1e-10);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "generic-alliance",
   "metadata": {},
   "source": [
    "Next we check how far the resulting `scfres` is from being a fixed point.\n",
    "Since iterative algorithms (e.g. diagonalisation) in interval arithmetic are tricky,\n",
    "we will ensure $\\rho = F(V(\\rho))$ by checking\n",
    "$$ H[\\rho] \\psi_i = \\varepsilon_i \\psi_i \\qquad \\text{for $\\rho = \\sum_i |\\psi_i|^2$},$$\n",
    "i.e. that the orbitals used for building a $\\rho$ and from this the Hamiltonian $H[\\rho]$\n",
    "still diagonalise $H[\\rho]$.\n",
    "If this is the case this implies that a subsequent diagonalisation of the $H[\\rho]$\n",
    "can only get me the same $\\rho$ back.\n",
    "\n",
    "The floating-point error in this computation of this check we will fully track using `IntervalArithmetic`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "interpreted-corner",
   "metadata": {},
   "outputs": [],
   "source": [
    "using GenericLinearAlgebra\n",
    "using IntervalArithmetic\n",
    "\n",
    "# TODO Missing in IntervalArithmetic for Float32:\n",
    "#      Interval ^ Interval not defined\n",
    "#      Fixing PR: https://github.com/JuliaIntervals/IntervalArithmetic.jl/pull/482\n",
    "import Base: ^\n",
    "^(a::Interval{Float32}, x::Interval) = IntervalArithmetic.atomic(Interval{Float32}, IntervalArithmetic.bigequiv(a)^x)\n",
    "# end \n",
    "\n",
    "# Get interval equivalents of key quantities\n",
    "intModel = model_DFT(Array{Interval{Float32}}(lattice), atoms, [:lda_x, :lda_c_vwn], symmetries=false)\n",
    "intBasis = PlaneWaveBasis(intModel; Ecut=7, kgrid=[4, 4, 4], fft_size=(16, 16, 16))\n",
    "intOccupation = [Interval.(occk) for occk in scfres.occupation]\n",
    "intEigvals = [Interval.(λk) for λk in scfres.eigenvalues]\n",
    "intEigvecs = [Interval.(ψk) for ψk in scfres.ψ]\n",
    "\n",
    "# Form density in interval arithmetic (i.e. evaluate F)\n",
    "intρ = DFTK.compute_density(intBasis, intEigvecs, intOccupation)\n",
    "@show maximum(radius, intρ)\n",
    "\n",
    "# Compute energy and Hamiltionian in interval arithmetic using this density (i.e. evaluate V)\n",
    "intEne, intHam = energy_hamiltonian(intBasis, intEigvecs, intOccupation;\n",
    "                                    ρ=intρ, eigenvalues=intEigvals,\n",
    "                                    εF=Interval(scfres.εF))\n",
    "@show radius(intEne.total)\n",
    "\n",
    "# Check the obtained eigenpairs are eigenpairs of this Hamiltonian\n",
    "# (i.e. check we are at a fixed point)\n",
    "residual_norm = similar(intEigvals)\n",
    "for ik in 1:length(intBasis.kpoints)\n",
    "    # Form Ritz values via Rayleigh quotient in interval arithmetic\n",
    "    Λks = intEigvecs[ik]' * (intHam.blocks[ik] * intEigvecs[ik])\n",
    "    residual_norm[ik] = norm.(eachcol(Λks - Diagonal(intEigvals[ik])))\n",
    "end\n",
    "\n",
    "n_converged = length(intEigvals[1]) - 3  # SCF contains 3 bands, which are not converged.\n",
    "@show maximum(maximum(mid, knorm[1:n_converged]) for knorm in residual_norm)\n",
    "@show maximum(maximum(radius, knorm[1:n_converged]) for knorm in residual_norm)\n",
    "nothing"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fuzzy-database",
   "metadata": {},
   "source": [
    "Clearly this calculation looses quite a bit of precision ...\n",
    "Unfortunately IntervalArithmetic only guarantees that the first digit of the energy\n",
    "can be trusted. Moreover it is not at all guaranteed that calculation has even converged!\n",
    "\n",
    "However, one should notice that interval arithmetic in general overestimates the floating-point error. For example in this case one can already see that the density computation already looses about 3-4 significant digits. If we improve upon this (e.g. compute the density in `Float64`), then we still have a reasonable number of trustworthy digits in the energy (about 3-4).\n",
    "\n",
    "#### Takeaway\n",
    "- Trustworthy DFT calculations in pure `Float32` are tricky.\n",
    "- IntervalArithmetic allows to identify routines where much precision is lost.\n",
    "- If in doubt DFTK allows you to check your DFT calculation in higher precision (e.g. `Double64`) ... but of course that comes at a cost."
   ]
  }
 ],
 "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
}