Crude import of some material from DFTK school (#941)

docs/.gitignore

@@ -5,6 +5,12 @@ src/guide/tutorial.md
 src/guide/tutorial.ipynb
 src/guide/periodic_problems.ipynb
 src/guide/periodic_problems.md
+src/guide/atomic_chains.ipynb
+src/guide/atomic_chains.md
+src/guide/self_consistent_field.ipynb
+src/guide/self_consistent_field.md
+src/guide/discretisation.ipynb
+src/guide/discretisation.md
 src/tricks/scf_checkpoints.ipynb
 src/tricks/scf_checkpoints.md
 /Artifacts.toml

docs/make.jl

@@ -17,11 +17,17 @@ PAGES = [
     "Home" => "index.md",
     "features.md",
     "Getting started" => [
-        # Installing DFTK, tutorial, theoretical background
         "guide/installation.md",
-        "Tutorial" => "guide/tutorial.jl",
-        "guide/periodic_problems.jl",
+        "guide/tutorial.jl",
+    ],
+    "Background" => [
+        # Theoretical background
         "guide/introductory_resources.md",
+        "guide/periodic_problems.jl",
+        "guide/discretisation.jl",
+        "guide/atomic_chains.jl",
+        "guide/density_functional_theory.md",
+        "guide/self_consistent_field.jl",
         "school2022.md",
     ],
     "Basic DFT calculations" => [
@@ -40,7 +46,6 @@ PAGES = [
     "Response and properties" => [
         "examples/polarizability.jl",
         "examples/forwarddiff.jl",
-        "examples/dielectric.jl",
     ],
     "Ecosystem integration" => [
         # This concerns the discussion of interfaces, IO and integration
@@ -60,6 +65,7 @@ PAGES = [
         "examples/custom_solvers.jl",
         "examples/scf_callbacks.jl",
         "examples/compare_solvers.jl",
+        "examples/analysing_scf_convergence.jl",
     ],
     "Nonstandard models" => [
         "examples/gross_pitaevskii.jl",

docs/src/guide/atomic_chains.jl

@@ -0,0 +1,163 @@
+# # Modelling atomic chains
+#
+# In [Periodic problems and plane-wave discretisations](@ref periodic-problems) we already
+# summarised the net effect of Bloch's theorem.
+# In this notebook, we will explore some basic facts about periodic systems,
+# starting from the very simplest model, a tight-binding monoatomic chain.
+# The solutions to the hands-on exercises are given at the bottom of the page.
+
+# ## Monoatomic chain
+#
+# In this model, each site of an infinite 1D chain is a degree of freedom, and
+# the Hilbert space is $\ell^2(\mathbb Z)$, the space of square-summable
+# biinfinite sequences $(\psi_n)_{n \in \mathbb Z}$.
+#
+# Each site interacts by a "hopping term" with its neighbors, and the
+# Hamiltonian is
+# ```math
+# H = \left(\begin{array}{ccccc}
+#   \dots&\dots&\dots&\dots&\dots \\
+#   \dots& 0 & 1 & 0 & \dots\\
+#   \dots&1 & 0 &1&\dots \\
+#   \dots&0 & 1 & 0& \dots  \\
+#   \dots&\dots&\dots&\dots&…
+# \end{array}\right)
+# ```
+#
+# !!! tip "Exercise 1"
+#     Find the eigenstates and eigenvalues of this Hamiltonian by
+#     solving the second-order recurrence relation.
+#
+# !!! tip "Exercise 2"
+#     Do the same when the system is truncated to a finite number of $N$
+#     sites with periodic boundary conditions.
+#
+# We are now going to code this:
+
+function build_monoatomic_hamiltonian(N::Integer, t)
+    H = zeros(N, N)
+    for n = 1:N-1
+        H[n, n+1] = H[n+1, n] = t
+    end
+    H[1, N] = H[N, 1] = t  # Periodic boundary conditions
+    H
+end
+
+# !!! tip "Exercise 3"
+#     Compute the eigenvalues and eigenvectors of this Hamiltonian.
+#     Plot them, and check whether they agree with theory.
+
+# ## Diatomic chain
+# Now we are going to consider a diatomic chain `A B A B ...`, where the coupling
+# `A<->B` ($t_1$) is different from the coupling `B<->A` ($t_2$). We will use a new
+# index $\alpha$ to denote the `A` and `B` sites, so that wavefunctions are now
+# sequences $(\psi_{\alpha n})_{\alpha \in \{1, 2\}, n \in \mathbb Z}$.
+#
+# !!! tip "Exercise 4"
+#     Show that eigenstates of this system can be looked for in the form
+#     ```math
+#        \psi_{\alpha n} = u_{\alpha} e^{ikn}
+#     ```
+#
+# !!! tip "Exercise 5"
+#     Show that, if $\psi$ is of the form above
+#     ```math
+#        (H \psi)_{\alpha n} = (H_k u)_\alpha e^{ikn},
+#     ```
+#     where
+#     ```
+#     H_k = \left(\begin{array}{cc}
+#     0                & t_1 + t_2 e^{-ik}\\
+#     t_1 + t_2 e^{ik} & 0
+#     \end{array}\right)
+#     ```
+#
+# Let's now check all this numerically:
+
+function build_diatomic_hamiltonian(N::Integer, t1, t2)
+    ## Build diatomic Hamiltonian with the two couplings
+    ## ... <-t2->   A <-t1-> B <-t2->   A <-t1-> B <-t2->   ...
+    ## We introduce unit cells as such:
+    ## ... <-t2-> | A <-t1-> B <-t2-> | A <-t1-> B <-t2-> | ...
+    ## Thus within a cell the A<->B coupling is t1 and across cell boundaries t2
+
+    H = zeros(2, N, 2, N)
+    A, B = 1, 2
+    for n = 1:N
+        H[A, n, B, n] = H[B, n, A, n] = t1  # Coupling within cell
+    end
+    for n = 1:N-1
+        H[B, n, A, n+1] = H[A, n+1, B, n] = t2  # Coupling across cells
+    end
+    H[A, 1, B, N] = H[B, N, A, 1] = t2  # Periodic BCs (A in cell1 with B in cell N)
+    reshape(H, 2N, 2N)
+end
+
+function build_diatomic_Hk(k::Integer, t1, t2)
+    ## Returns Hk such that H (u e^ikn) = (Hk u) e^ikn
+    ##
+    ## intra-cell AB hopping of t1, plus inter-cell hopping t2 between
+    ## site B (no phase shift) and site A (phase shift e^ik)
+    [0                 t1 + t2*exp(-im*k);
+     t1 + t2*exp(im*k) 0                 ]
+end
+
+using Plots
+function plot_wavefunction(ψ)
+    p = plot(real(ψ[1:2:end]), label="Re A")
+    plot!(p, real(ψ[2:2:end]), label="Re B")
+end
+
+# !!! tip "Exercise 6"
+#     Check the above assertions. Use a $k$ of the form
+#     $2 π \frac{l}{N}$ in order to have a $\psi$ that has the periodicity
+#     of the supercell ($N$).
+
+# !!! tip "Exercise 7"
+#     Plot the band structure, i.e. the eigenvalues of $H_k$ as a function of $k$
+#     Use the function `build_diatomic_Hk` to build the Hamiltonians.
+#     Compare with the eigenvalues of the ("supercell") Hamiltonian from
+#     `build_diatomic_hamiltonian`. In the case $t_1 = t_2$, how do the bands follow
+#     from the previous study of the monoatomic chain?
+
+# !!! tip "Exercise 8"
+#     Repeat the above analysis in the case of a finite-difference
+#     discretization of a continuous Hamiltonian $H = - \frac 1 2 \Delta + V(x)$
+#     where $V$ is periodic
+#     *Hint:* It is advisable to work through [Comparing discretization techniques](@ref)
+#     before tackling this question.
+
+# ## Solutions
+#
+# ### Exercise 1
+# !!! note "TODO"
+#     This solution has not yet been written. Any help with a PR is appreciated.
+#
+# ### Exercise 2
+# !!! note "TODO"
+#     This solution has not yet been written. Any help with a PR is appreciated.
+#
+# ### Exercise 3
+# !!! note "TODO"
+#     This solution has not yet been written. Any help with a PR is appreciated.
+#
+# ### Exercise 4
+# !!! note "TODO"
+#     This solution has not yet been written. Any help with a PR is appreciated.
+#
+# ### Exercise 5
+# !!! note "TODO"
+#     This solution has not yet been written. Any help with a PR is appreciated.
+#
+# ### Exercise 6
+# !!! note "TODO"
+#     This solution has not yet been written. Any help with a PR is appreciated.
+#
+# ### Exercise 7
+# !!! note "TODO"
+#     This solution has not yet been written. Any help with a PR is appreciated.
+#
+# ### Exercise 8
+# !!! note "TODO"
+#     This solution has not yet been written. Any help with a PR is appreciated.
+

docs/src/guide/density_functional_theory.md

@@ -0,0 +1,64 @@
+# Introduction to density-functional theory
+
+!!! note "More details"
+    This chapter only gives a very rough overview for now. Further details
+    can be found in the [summary of DFT theory](https://michael-herbst.com/teaching/2022-mit-workshop-dftk/2022-mit-workshop-dftk/DFT_Theory.pdf)
+    or in the [Introductory Resources](@ref introductory-resources).
+
+Density-functional theory is a simplification of the full electronic
+Schrödinger equation leading to an effective single-particle model.
+Mathematically this can be cast as an energy minimisation problem
+in the electronic density $\rho$. In the Kohn-Sham variant
+of density-functional theory the corresponding first-order stationarity
+conditions can be written as the non-linear eigenproblem
+```math
+\begin{aligned}
+&\left( -\frac12 \Delta + V\left(\rho\right) \right) \psi_i = \varepsilon_i \psi_i, \\
+V(\rho) = &\, V_\text{nuc} + v_C \rho + V_\text{XC}(\rho), \\
+\rho = &\sum_{i=1}^N f(\varepsilon_i)  \abs{\psi_i}^2, \\
+\end{aligned}
+```
+where $\{\psi_1,\ldots, \psi_N\}$ are $N$ orthonormal orbitals,
+the one-particle eigenfunctions and $f$ is the occupation function
+(e.g. [`DFTK.Smearing.FermiDirac`](@ref), [`DFTK.Smearing.Gaussian`](@ref)
+or [`DFTK.Smearing.MarzariVanderbilt`](@ref))
+chosen such that the integral of the density is equal to the number of electrons in the system.
+Further the potential terms that make up $V(\rho)$ are
+- the nuclear attraction potential $V_\text{nuc}$ (interaction of electrons and nuclei)
+- the exchange-correlation potential $V_\text{xc}$,
+  depending on $\rho$ and potentially its derivatives.
+- The Hartree potential $v_C \rho$, which is obtained by solving the Poisson equation
+  ```math
+  -\Delta \left(v_C \rho\right) = 4\pi \rho
+  ```
+The non-linearity is such due to the fact that the DFT Hamiltonian
+```math
+H = -\frac12 \Delta + V\left(\rho\right)
+```
+depends on the density $\rho$, which is built from its own eigenfunctions.
+Often $H$ is also called the Kohn-Sham matrix or Fock matrix.
+
+Introducing the *potential-to-density map*
+```math
+D(V) = \sum_{i=1}^N f(\varepsilon_i)  \abs{\psi_i}^2
+\qquad \text{where} \qquad \left(-\frac12 \Delta + V\right) \psi_i = \varepsilon_i \psi_i
+```
+allows to write the DFT problem in the short-hand form
+```math
+\rho = D(V(\rho)),
+```
+which is a fixed-point problem in $\rho$, also known as the
+**self-consistent field problem** (SCF problem).
+Notice that computing $D(V)$ requires the diagonalisation of the operator
+$-\frac12 \Delta + V$, which is usually the most
+expensive step when solving the SCF problem.
+
+To solve the above SCF problem $ \rho = D(V(\rho)) $
+one usually follows an iterative procedure.
+That is to say that starting from an initial guess $\rho_0$ one then
+computes $D(V(\rho_n))$ for a sequence of iterates $\rho_n$ until input and output
+are close enough. That is until the residual
+```math
+R(\rho_n) = D(V(\rho_n)) - \rho_n
+```
+is small. Details on such SCF algorithms will be discussed in [Self-consistent field methods](@ref).