7_Exercises.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "configured-charity",
   "metadata": {},
   "source": [
    "# Exercises \n",
    "\n",
    "These exercises provide a short look at the \"numerical\" and \"computational\" challenges of density-functional theory (DFT). For the work here we will use the [density-functional toolkit](https://dftk.org) (DFTK), a small Julia code for plane-wave DFT simulations. \n",
    "\n",
    "Some modification of Julia code is required to do the exercises, but apart from that you can do the analysis of your numerical results in any form you want. But naturally this would be a good opportunity to learn some Julia ;).\n",
    " \n",
    "Some help in case you get stuck and feel you would like to learn more:\n",
    "- https://michael-herbst.com/learn-julia: Two-day Julia introductory course. Highly recommended are the notebooks on\n",
    "    * [Variables and control structures](https://nbviewer.jupyter.org/github/mfherbst/2022-rwth-julia-workshop/blob/master/01_Variables_Control_Packages.ipynb)\n",
    "    * [Functions](https://nbviewer.jupyter.org/github/mfherbst/2022-rwth-julia-workshop/blob/master/02_Functions_Types_Dispatch.ipynb)\n",
    "    * [Arrays](https://nbviewer.jupyter.org/github/mfherbst/2022-rwth-julia-workshop/blob/master/03_Arrays_Parametric_Types.ipynb)\n",
    "- https://julialang.org/learning/: Various learning resources for Julia\n",
    "- https://docs.julialang.org/: Julia documentation\n",
    "\n",
    "\n",
    "## Installation of Julia and required packages\n",
    "\n",
    "For installation instructions on the packages required to run these exercises and some tips on using Julia on the RWTH Jupyter lab setup, see [1_Installation.ipynb](1_Installation).\n",
    "\n",
    "By default Julia only using a single thread. But you can also create a *Jupyter kernel* for multithreaded Julia. It is highly recommended you do this for the exercises, see [1_Installation.ipynb](1_Installation) for details."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "equivalent-signal",
   "metadata": {},
   "source": [
    "## Setting up DFTK \n",
    "\n",
    "Do this every time you start the notebook to initialise the threading inside DFTK:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "integrated-pregnancy",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Setup DFTK\n",
    "using DFTK\n",
    "setup_threading(n_blas=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "round-lindsay",
   "metadata": {},
   "source": [
    "## Exercise 1: Why do we need pseudopotentials? (2 points)\n",
    "\n",
    "It was discussed in the lecture that plane-wave DFT codes usually employ so-called pseudopotentials (PSP) to model the interaction of electrons and nuclei. Replacing the true Coulombic interaction by a PSP is an approximation. Multiple types of PSPs exist and from a mathematical perspective little is understood about the errors introduced by PSPs. However, from a physical point of view using PSPs is reasonable as it only effects the electron density close to the nuclei, which for most phaenomena only plays a minor role. On the upside PSPs turn out to make plane-wave calculations much more feasible.\n",
    "\n",
    "This aspect we will investigate numerically in this exercise. Our goal will be to determine the energy of a single neon atom using the so-called LDA approximation up to an absolute error of `0.1`. We will do this by an explicit convergence study, i.e. by increasing the basis size and improving the numerics until we are sure we have found the energy to this tolerance.\n",
    "\n",
    "Our computational setup is simple: We will put a single neon atom into a cubic box of length $a$. As DFTK uses periodic boundary conditions, this atom interacts spuriously with the neighbouring cells, introducing an error to our calculation. As $a \\to \\infty$ this error disappears, thus in principle convergence with increasing $a$ should be one parameter to study. For simplicitity (as large values of $a$ lead to very costly calculations) we will ignore this aspect here.\n",
    "\n",
    "The convergence parameter we will not ignore, however, is the $E_\\text{cut}$ parameter discussed in the lecture. This parameter determines the basis set size, thus the accuracy of the calculation. Again calculations get better for $E_\\text{cut} \\to \\infty$.\n",
    "\n",
    "In the language of DFTK the calculation we want to perform is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "monthly-colony",
   "metadata": {},
   "outputs": [],
   "source": [
    "using LinearAlgebra\n",
    "\n",
    "# Numerical parameters\n",
    "#\n",
    "Ecut = 20  # Plane-wave basis cutoff\n",
    "\n",
    "# Use a HGH pseudopotential:\n",
    "Ne = ElementPsp(:Ne, psp=load_psp(\"hgh/lda/ne-q8\"))\n",
    "\n",
    "# Use the true Coulomb interaction:\n",
    "# Ne = ElementCoulomb(:Ne)\n",
    "#\n",
    "# End numerical parameters \n",
    "\n",
    "# Setup Cubic box and place neon atom in the middle\n",
    "a = 10    # Box size \n",
    "lattice   = a * Matrix(I, 3, 3)\n",
    "atoms     = [Ne]\n",
    "positions = [[0.5, 0.5, 0.5]]\n",
    "\n",
    "# Use LDA model, discretise and run SCF algorithm\n",
    "model  = model_LDA(lattice, atoms, positions)\n",
    "basis  = PlaneWaveBasis(model; Ecut, kgrid=[1, 1, 1])\n",
    "scfres = self_consistent_field(basis)\n",
    "Etot   = scfres.energies.total\n",
    "\n",
    "# Print total energy\n",
    "println(\"\\nTotal DFT energy is $Etot\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "industrial-shuttle",
   "metadata": {},
   "source": [
    "**Warning:** DFT calculations can be both time and memory consuming. When doing these convergence studies therefore start with small values for `Ecut` (like the ones shown here) and increase slowly, especially if you are running these on a laptop with 8GB of main memory or less. You are not expected to re-compute the provided reference results."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ignored-medium",
   "metadata": {},
   "source": [
    "(a) First stick with the pseudpotential (PSP) version of the neon atom (`ElementPsp`). Converge the energy to an absolute error of `0.1`. The easiest way to do this is to run multiple calculations for different values of `Ecut` and to plot the error against a reference on a semilog scale. A good reference is at `Ecut = 600`, where the DFT energy is `-34.85376`.<br />\n",
    "*Hint:* If you want to do the ploting within Julia good, take a look at the [Plots.jl](http://docs.juliaplots.org/latest/generated/gr/) documentation for some nice examples. Another good plotting package could be [PyPlot.jl](https://github.com/JuliaPy/PyPlot.jl), which has exactly the same interface as the `matplotlib` Python package.\n",
    "\n",
    "(b) Repeat the exercise for the all-electron case (`ElementCoulomb`). Here the reference result at `Ecut = 600` is `-123.5983`. Explore the convergence all the way up to `Ecut=100`. What differences in the obtained total energy as well as the convergence with `Ecut` do you observe? Suggest an explanaition keeping in mind that part of the idea of the PSP is to avoid the explicit treatment of the core electrons."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "incoming-niagara",
   "metadata": {},
   "source": [
    "## Exercise 2: SCF algorithms to solve DFT (3 points)\n",
    "\n",
    "As we discussed in [5_Solving_the_SCF_problem.ipynb](5_Solving_the_SCF_problem.ipynb) the self-consistent field procedure required to solve the DFT problem can be written as a fixed-point problem\n",
    "$$ F(\\rho) = \\rho $$\n",
    "where $F$ represents the basic SCF step. That is the construction of the Kohn-Sham Hamiltonian $H(\\rho)$ given the density $\\rho$, followed its diagonalisation to obtain its eigenpairs $(\\varepsilon_{k i}, \\psi_{ki})$\n",
    "and from these a new density\n",
    "$$ \\rho(r) = \\sum_{k\\in\\text{BZ}} \\sum_i f_{\\varepsilon_F}(\\varepsilon_{ki}) \\, \\psi_{ki}(r) \\, \\psi_{ki}^\\ast(r)$$\n",
    "with the Fermi level $\\varepsilon_F$ chosen to conserve the number of electrons.\n",
    "\n",
    "In this exercise we will take the function $F$ for \"granted\" (i.e. delivered by DFTK) and we will investigate multiple algorthms to find the fixed-point density satisfying $F(\\rho) = \\rho $. Our setting is defined by the following function, which builds and discretises an LDA model for an aluminium supercell, see [5_Solving_the_SCF_problem.ipynb](5_Solving_the_SCF_problem.ipynb) for details."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "272b1196",
   "metadata": {},
   "outputs": [],
   "source": [
    "using DFTK\n",
    "using LinearAlgebra\n",
    "\n",
    "function aluminium_setup_ex2(repeat=1; Ecut=13.0, kgrid=[2, 2, 2])\n",
    "    a = 7.65339\n",
    "    lattice = a * Matrix(I, 3, 3)\n",
    "    Al = ElementPsp(:Al, psp=load_psp(\"hgh/lda/al-q3\"))\n",
    "    atoms     = [Al, Al, Al, Al]\n",
    "    positions = [[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",
    "    # Make supercell\n",
    "    supercell = ase_atoms(lattice, atoms, positions) * (repeat, 1, 1)\n",
    "    lattice   = load_lattice(supercell)\n",
    "    positions = load_positions(supercell)\n",
    "    atoms = fill(Al, length(positions))\n",
    "\n",
    "    # Construct an LDA model and discretise\n",
    "    model = model_LDA(lattice, atoms, positions; temperature=1e-3, symmetries=false)\n",
    "    PlaneWaveBasis(model; Ecut, kgrid)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "differential-spyware",
   "metadata": {},
   "source": [
    "**(a) Fixed-point iterations.** The easiest way to solve $F(\\rho) = \\rho$ are fixed-point iterations, i.e.\n",
    "$$ \\rho_{n+1} = F(\\rho_n), $$\n",
    "starting from a hopefully good initial guess $\\rho_0$. DFTK automatically provides a reasonable\n",
    "guess density, such that we only need to take care of the iterations themselves.\n",
    "In the language of DFTK this algorithm is written as:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "advance-brunswick",
   "metadata": {},
   "outputs": [],
   "source": [
    "function fixed_point_iteration(F, ρ₀, maxiter; tol)\n",
    "    # F:        The SCF step function\n",
    "    # ρ₀:       The initial guess density\n",
    "    # maxiter:  The maximal number of iterations to be performed\n",
    "    # tol:      The selected convergence tolerance\n",
    "    \n",
    "    ρ  = ρ₀\n",
    "    Fρ = F(ρ)\n",
    "    for n = 1:maxiter \n",
    "        # If change less than tolerance, break iterations:\n",
    "        if norm(Fρ - ρ) < tol\n",
    "            break\n",
    "        end\n",
    "        ρ  = Fρ\n",
    "        Fρ = F(ρ)\n",
    "    end\n",
    "    \n",
    "    # Return some stuff DFTK needs ...\n",
    "    (fixpoint=ρ, converged=norm(Fρ-ρ) < tol)\n",
    "end;\n",
    "\n",
    "# use this algorithm inside DFTK's SCF for solving the silicon problem\n",
    "# (the other parameters are needed to overwrite some DFTK defaults we don't want to use yet).\n",
    "self_consistent_field(aluminium_setup_ex2(1); solver=fixed_point_iteration, damping=1.0,\n",
    "                                              maxiter=40, mixing=SimpleMixing());"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "blocked-acting",
   "metadata": {},
   "source": [
    "As can be observed this algorithm is not very good and in fact even fails to converge albeit we are only looking at a very simple system.\n",
    "\n",
    "This is a known limitation of this algorithm, which is why it is not used in practice. Instead one introduces a so-called damping parameter $\\alpha$, which is given a value between $0$ and $1$. One now iterates as follows:\n",
    "$$ \\rho_{n+1} = \\rho_{n} + \\alpha (F(\\rho_n) - \\rho_n) $$\n",
    "In other words the update $F(\\rho_n) - \\rho_n$ proposed in the $n$-th SCF step is not fully taken, but scaled-down by the damping $\\alpha$.\n",
    "\n",
    "Modify `fixed_point_iteration` such that it supports this *damped* fixed-point iteration. In other words implement damping *inside* your algorithm and not by changing the `damping` parameter of the `self_consistent_field` function driving the SCF.  \n",
    "Using your algorithm try different values for $\\alpha$ between $0$ and $1$ and estimate roughly the $\\alpha$ which gives fastest convergence. For which $\\alpha$ do you observe no convergence at all?\n",
    "\n",
    "**Remark:** The observations you make here are general. One can proove that every SCF converges (locally) if a small enough $\\alpha > 0$ is chosen."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4276a981",
   "metadata": {},
   "source": [
    "**(b) Anderson acceleration.** The `fixed_point_iteration` function above (with the damping extension) actually already covers the main gist of the DFT algorithms used in practice. One additional idea to make things converge faster is Anderson acceleration, where not only $\\rho_n$ and $F(\\rho_n)$, but also older iterates are used to propose the next density.\n",
    "\n",
    "For Anderson one exploits that the update $R(\\rho) = F(\\rho) - \\rho$ is also the residual of the fixed-point problem $F(\\rho) = \\rho$, i.e. how far away we are from the fixed-point density. A good next density $\\rho_{n+1}$ therefore should be found by minimising an approximation for $R(\\rho_{n+1})$. Assuming the SCF was linear in the density (which it is not), a good idea is to find a linear combination of residuals\n",
    "$$\\min_{\\beta_i} \\left\\| \\sum_i \\beta_i R(\\rho_i) \\right\\|^2$$\n",
    "which has the smallest possible norm and to use these coefficients $\\beta_i$ to extrapolate the next\n",
    "density\n",
    "$$ \\rho_{n+1} =  \\sum_i \\beta_i (\\rho_i + \\alpha R(\\rho_i)) $$\n",
    "where you notice the \"standard\" damped fixed-point iteration in the summed terms.\n",
    "\n",
    "In terms of an algorithm this is"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0f3b2c79",
   "metadata": {},
   "outputs": [],
   "source": [
    "function anderson_iteration(F, ρ₀, maxiter; tol)\n",
    "    # F:        The SCF step function\n",
    "    # ρ₀:       The initial guess density\n",
    "    # maxiter:  The maximal number of iterations to be performed\n",
    "    # tol:      The selected convergence tolerance\n",
    "    \n",
    "    converged = false\n",
    "    ρ  = ρ₀\n",
    "    ρs = []\n",
    "    Rs = []\n",
    "    for n = 1:maxiter\n",
    "        Fρ = F(ρ)\n",
    "        Rρ = Fρ - ρ\n",
    "        converged = norm(Rρ) < tol\n",
    "        converged && break\n",
    "        \n",
    "        ρnext = vec(ρ) .+ vec(Rρ)\n",
    "        if !isempty(Rs)\n",
    "            M = hcat(Rs...) .- vec(Rρ)\n",
    "            βs = -(M \\ vec(Rρ))\n",
    "            \n",
    "            for (iβ, β) in enumerate(βs)\n",
    "                ρnext .+= β .* (ρs[iβ] .- vec(ρ) .+ Rs[iβ] .- vec(Rρ))\n",
    "            end\n",
    "        end\n",
    "                    \n",
    "        push!(ρs, vec(ρ))\n",
    "        push!(Rs, vec(Rρ))\n",
    "        ρ = reshape(ρnext, size(ρ₀)...)\n",
    "    end\n",
    "    \n",
    "    # Return some stuff DFTK needs ...\n",
    "    (fixpoint=ρ, converged=converged)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c405d1b3",
   "metadata": {},
   "source": [
    "To work with this algorithm we will use DFTK's intrinsic mechanism to choose a damping. The syntax for this is\n",
    "\n",
    "```julia\n",
    "repeat = 1\n",
    "self_consistent_field(aluminium_setup_ex2(repeat);\n",
    "                      solver=anderson_iteration,\n",
    "                      damping=0.8, maxiter=40,\n",
    "                      mixing=SimpleMixing());\n",
    "```\n",
    "to choose a damping of $\\alpha = 0.8$ and run for at most `maxiter` iterations.\n",
    "\n",
    "(i) Based on this Anderson implementation verify (by making a few experiments) that the algorithm converges for `repeat=1` for any $0 < \\alpha \\leq 2$. You may now use the `damping` parameter for changing the value $\\alpha$ used by the SCF. State the number of iterations and runtimes you observe.\n",
    "\n",
    "(ii) Pick $\\alpha = 0.8$ and make the problem harder by increasing `repeat` (e.g. `2`, `4`, `6`, `8`). Can you make Anderson fail to converge? What do you notice in terms of the number of iterations and runtimes?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "prescription-morrison",
   "metadata": {},
   "source": [
    "**(c) Mixing methods.** Anderson allows us to push the boundary for SCF methods, but for larger or more challenging systems it is not fully sufficient. The next ingredient for a stable SCF procedure is based on the insight that the convergence properties of an SCF provably depend on the dielectric properties of materials, which is simulated. Amongst others this is to say that insulators (like glass), semiconductors (like silicon) or metals (like aluminium) have rather differing SCF behaviours. As a result the ideal SCF procedure should be slightly different for each material.\n",
    "\n",
    "The standard approach to include material specificity into the SCF is to employ *preconditioned* damped fixed-point iterations.\n",
    "To explain the idea, let us consider again a framework without Anderson acceleration.\n",
    "Preconditioned iterations are then\n",
    "$$ \\rho_{n+1} = \\rho_{n} + \\alpha P^{-1} (F(\\rho_n) - \\rho_n), $$\n",
    "where $P^{-1}$ is a preconditioner that hopefully improve convergence.\n",
    "To re-introduce Anderson around this iteration\n",
    "just replace the previous definition of $R$ by\n",
    "$R(\\rho) = P^{-1} (F(\\rho_n) - \\rho_n)$.\n",
    "\n",
    "Finding a good preconditioner $P$ is not always easy and for some cases satisfactory options are not yet known. For our aluminium case, however, we are lucky. The `KerkerMixing` implemented in DFTK provides a good $P$ for aluminium.\n",
    "\n",
    "You might wonder about the term *mixing*. In an interdisciplinary community like DFT, different scientists use different vocabulary and \"mixing\" is the physicists' term used for preconditioning.\n",
    "\n",
    "To use `KerkerMixing` with DFTK run the SCF as follows\n",
    "```julia\n",
    "self_consistent_field(basis; damping=0.8, mixing=KerkerMixing());\n",
    "```\n",
    "\n",
    "Try this setup for different values of `repeat` and check the number of iterations needed. Other mixings DFTK has to offer are `DielectricMixing` (best for semiconductors), `SimpleMixing` (which is $P = I$, i.e. no preconditioner at all, best for insulators) or `LdosMixing` (self-adapting, suitable for both metals *or* insulators *or* inhomogeneous mixtures). Note that `LdosMixing` is the default in DFTK (i.e. used if the `mixing` parameter is *not* supplied to `self_consistent_field`. Try these mixings (`SimpleMixing`, `DielectricMixing`, `LdosMixing` and `KerkerMixing`) and summarise your findings.\n",
    "\n",
    "You should notice that choosing a preconditioner matching the material under study aids a fast SCF convergence, but that sometimes being off does not seem to do much harm for our case. For larger values of `repeat` (beyond what you can probably effort on your laptop) this is no longer true and one needs to be very careful in selecting the right preconditioner. See for example the investigation in [this paper](https://michael-herbst.com/publications/2020.09.03_ldos_preconditioning.pdf). "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1bded2a9",
   "metadata": {},
   "source": [
    "## Exercise 3: Analysing SCF convergence (2 points)\n",
    "\n",
    "The goal of this exercise is to explain the differing convergence behaviour of SCF algorithms depending on the choice of the preconditioner $P^{-1}$ and the underlying material.\n",
    "\n",
    "For this we will find the largest and smallest eigenvalue of $(P^{-1} \\varepsilon^\\dagger)$ and $\\varepsilon^\\dagger$, where $\\varepsilon^\\dagger$ is the dielectric operator (see [6_Analysing_SCF_convergence.ipynb](6_Analysing_SCF_convergence.ipynb)). The ratio of largest to smallest eigenvalue is the condition number\n",
    "  $$ \\kappa = \\frac{\\lambda_\\text{max}}{\\lambda_\\text{min}},$$\n",
    "which can be related to the rate of convergence. The smaller the condition number, the faster the convergence.\n",
    "\n",
    "**(a) Aluminium metal.** We start by taking a look at a slightly cruder (thus computationally cheaper) version of our aluminium setup from above: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b383f4b8",
   "metadata": {},
   "outputs": [],
   "source": [
    "using DFTK\n",
    "using LinearAlgebra\n",
    "\n",
    "function aluminium_setup_ex3(repeat=1; Ecut=7.0, kgrid=[1, 1, 1])\n",
    "    a = 7.65339\n",
    "    lattice = diagm(fill(a, 3))\n",
    "    Al = ElementPsp(:Al, psp=load_psp(\"hgh/lda/al-q3\"))\n",
    "    atoms = [Al, Al, Al, Al]\n",
    "    positions = [[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",
    "    # Make supercell\n",
    "    supercell = ase_atoms(lattice, atoms, positions) * (repeat, 1, 1)\n",
    "    lattice   = load_lattice(supercell)\n",
    "    positions = load_positions(supercell)\n",
    "    atoms     = fill(Al, length(positions))\n",
    "\n",
    "    # Construct the model\n",
    "    model = model_LDA(lattice, atoms, positions, temperature=1e-3, symmetries=false)\n",
    "    PlaneWaveBasis(model; Ecut, kgrid)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c16ce4d0",
   "metadata": {},
   "source": [
    "Additionally we will need the function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6082fea1",
   "metadata": {},
   "outputs": [],
   "source": [
    "function construct_Pinv_epsilon(scfres)\n",
    "    basis = scfres.basis\n",
    "    \n",
    "    Pinv_Kerker(δρ) = DFTK.mix_density(KerkerMixing(), basis, δρ)\n",
    "    \n",
    "    function epsilon(δρ)  # Apply ε† = 1 - χ0 Khxc\n",
    "        δV   = apply_kernel(basis, δρ; ρ=scfres.ρ)\n",
    "        χ0δV = apply_χ0(scfres.ham, scfres.ψ, scfres.εF, scfres.eigenvalues, δV)\n",
    "        δρ - χ0δV   \n",
    "    end    \n",
    "    \n",
    "    epsilon, Pinv_Kerker\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "957f4589",
   "metadata": {},
   "source": [
    "To construct functions representing $\\varepsilon^\\dagger$ and the Kerker preconditioner $P^{-1}$.\n",
    "\n",
    "(i) Find the largest eigenvalue of $\\varepsilon^\\dagger$ for this aluminium case using `KrylovKit`. For this use the following code snippet:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "561c9477",
   "metadata": {},
   "outputs": [],
   "source": [
    "using KrylovKit\n",
    "\n",
    "scfres = self_consistent_field(aluminium_setup_ex3(3); tol=1e-12);\n",
    "epsilon, Pinv_Kerker = construct_Pinv_epsilon(scfres)\n",
    "\n",
    "λ_large, X_large, info = eigsolve(epsilon, randn(size(scfres.ρ)), 4, :LM;\n",
    "                                  tol=1e-4, eager=true, verbosity=2)\n",
    "@assert info.converged ≥ 4\n",
    "λ_max = maximum(real.(λ_large))\n",
    "\n",
    "println(\"Largest eigenvalue: $(λ_max)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b8b0ce41",
   "metadata": {},
   "source": [
    "You can assume the smallest eigenvalue is $λ_{min} = 0.952$. What is the condition number in this case?\n",
    "\n",
    "(ii) Find the largest eigenvalue for the SCF of the aluminium supercell (`repeat=3`) in case the Kerker preconditioner is used.  \n",
    "*Hint:* You can construct the operator $P^{-1} \\varepsilon^\\dagger$ by simply chaining the functions (`Pinv_Kerker ∘ epsilon`). Assuming that the smallest eigenvalue is about $0.8$, what is the condition number now? \n",
    "\n",
    "**(b) Helium chain (insulator).** We will use the following setup:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "57dd7ef2",
   "metadata": {},
   "outputs": [],
   "source": [
    "using DFTK\n",
    "using LinearAlgebra\n",
    "\n",
    "function helium_setup(repeat=30; Ecut=7.0, kgrid=[1, 1, 1])\n",
    "    a = 5.0\n",
    "    lattice = a * Matrix(I, 3, 3)\n",
    "    He = ElementPsp(:He, psp=load_psp(\"hgh/lda/he-q2\"))\n",
    "    atoms = [He, ]\n",
    "    positions = [[0.0, 0.0, 0.0], ]\n",
    "\n",
    "    # Make supercell \n",
    "    supercell = ase_atoms(lattice, atoms, positions) * (repeat, 1, 1)\n",
    "    lattice   = load_lattice(supercell)\n",
    "    positions = load_positions(supercell)\n",
    "    atoms     = fill(He, length(positions))\n",
    "\n",
    "    # Construct the model\n",
    "    model = model_LDA(lattice, atoms, positions, temperature=1e-3, symmetries=false)\n",
    "    PlaneWaveBasis(model; Ecut, kgrid)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc053c65",
   "metadata": {},
   "source": [
    "Repeat the analysis from (a) for a Helium chain with `repeat=30`. To find the smallest and largest eigenvalues of $\\varepsilon^\\dagger$ and $P^{-1} \\varepsilon^\\dagger$ use"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cebba7f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "using KrylovKit\n",
    "\n",
    "scfres = self_consistent_field(helium_setup(30); tol=1e-12);\n",
    "epsilon, Pinv_Kerker = construct_Pinv_epsilon(scfres)\n",
    "\n",
    "operator = epsilon\n",
    "\n",
    "λ_large, X_large, info = eigsolve(operator, randn(size(scfres.ρ)), 2, :LM;\n",
    "                                  tol=1e-3, eager=true, verbosity=2)\n",
    "@assert info.converged ≥ 2\n",
    "λ_max = maximum(real.(λ_large))\n",
    "    \n",
    "λ_small, X_small, info = eigsolve(operator, randn(size(scfres.ρ)), 2, EigSorter(abs, rev=false);\n",
    "                                  tol=1e-3, eager=true, verbosity=2)\n",
    "@assert info.converged ≥ 2\n",
    "λ_min = minimum(real.(λ_small))\n",
    "\n",
    "println(\"Smallest eigenvalue: $(λ_min)\")\n",
    "println(\"Largest eigenvalue:  $(λ_max)\")\n",
    "println(\"Condition number:    $(λ_max / λ_min)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a09712d7",
   "metadata": {},
   "source": [
    "Then run the two SCFs with and without Kerker preconditioning, that is"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e35334b7",
   "metadata": {},
   "outputs": [],
   "source": [
    "scfres = self_consistent_field(helium_setup(30); tol=1e-12, mixing=SimpleMixing());"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae8252d2",
   "metadata": {},
   "source": [
    "as well as"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9c80a0d2",
   "metadata": {},
   "outputs": [],
   "source": [
    "scfres = self_consistent_field(helium_setup(30); tol=1e-12, mixing=KerkerMixing());"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3a8afa19",
   "metadata": {},
   "source": [
    "and explain the observations with respect to convergence, taking your findings on the eigenvalues of $\\varepsilon^\\dagger$ and $P^{-1} \\varepsilon^\\dagger$ into account."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "electrical-campaign",
   "metadata": {},
   "source": [
    "## Exercise 4: Why iron is magnetic and silicon is not (3 points)\n",
    "\n",
    "Only few compounds exhibit a natural permanent magnetism. One of these is iron, while most (like silicon) are not magnetic. This exercise tries to provide a little insight how one could understand, using DFT simulations, why this is the case.\n",
    "\n",
    "The key assumption in this exercise will be that DFT is a realistic model and that the SCF therefore finds a good approximation to the true ground state of a compound. If this ground state turns out to be magnetic, the compound should therefore feature a permanant magnetism.\n",
    "\n",
    "We use this computational setup to simulate silicon:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "rough-fluid",
   "metadata": {},
   "outputs": [],
   "source": [
    "using DFTK\n",
    "\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, Si]\n",
    "positions = [-ones(3)/8, ones(3)/8]\n",
    "\n",
    "# Guess some inital magnetic moments\n",
    "# (Need to be != 0 otherwise SCF makes the assumption that the ground state is not magnetic\n",
    "#  to gain some performance ...)\n",
    "magnetic_moments = [2, 2]\n",
    "\n",
    "model  = model_LDA(lattice, atoms, positions; temperature=0.01, magnetic_moments)\n",
    "basis  = PlaneWaveBasis(model; Ecut=13, kgrid=[2, 2, 2]);\n",
    "ρ0     = guess_density(basis, magnetic_moments)\n",
    "scfres = self_consistent_field(basis; ρ=ρ0);\n",
    "scfres.energies.total"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "hundred-contact",
   "metadata": {},
   "source": [
    "Even though we have forced some magnetism into the initial density guess, this magnetisation (indicated by `Magnet`) disappears as the SCF converges. Therefore silicon appears to have a non-magnetic ground state."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "daily-white",
   "metadata": {},
   "source": [
    "(i) Repeat the calculation for iron. In this case the system setup is\n",
    "```julia\n",
    "a = 5.42352  # iron lattice constant in bohr\n",
    "lattice = a / 2 * [[-1  1  1];\n",
    "                   [ 1 -1  1];\n",
    "                   [ 1  1 -1]]\n",
    "Fe = ElementPsp(:Fe, psp=load_psp(\"hgh/lda/Fe-q8.hgh\"))\n",
    "atoms = [Fe]\n",
    "positions = [zeros(3)]\n",
    "magnetic_moments = [3]\n",
    "```\n",
    "otherwise use the same settings as in the silicon calculation. Based on this calculation, what would you conclude?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "pleased-translation",
   "metadata": {},
   "source": [
    "(ii) As it turns out too small basis sets and Brilouin-zone sampling (too small `kgrid`s) are not able to support magnetic ground states. Repeat both the silicon as well as the iron calculations for different values of `Ecut` and `kgrid` (i.e. `[1,1,1]`,`[3,3,3]`, `[4,4,4]` etc.), where in both cases larger means better. Play with these parameters to determine the first two digits of the ground state energy of silicon and iron. Based on these numerical parameters what do you conclude now?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "heated-property",
   "metadata": {},
   "source": [
    "(iii) To show that a non-magnetic iron structure is actually higher in energy re-run the iron calculation with the `Ecut` and `kgrid` parameters determined in (ii), but this time set `magnetic_moments = [0]`. This enforces the SCF to converge to a non-magnetic ground state. This is why the magnetisation `Magnet` is no longer printed ... it is `0` by assumption. What is the energy difference between this non-magnetic iron ground state and the ground state you determined in (ii)?"
   ]
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "Julia (4 threads) 1.7.3",
   "language": "julia",
   "name": "julia-(4-threads)-1.7"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}