{ "cells": [ { "cell_type": "markdown", "id": "latest-correspondence", "metadata": {}, "source": [ "# Periodic problems and plane-wave discretisations\n", "\n", "Before digging into big real-world problems, i.e. solving density-functional theory in 3D, we will first take a step back and briefly look at some specifics of periodic problems and plane-wave discretisations in one dimension." ] }, { "cell_type": "markdown", "id": "afc6f9b1", "metadata": {}, "source": [ "## Periodicity and lattices\n", "A periodic problem is characterized by being invariant to certain translations.\n", "For example the $\\sin$ function is periodic with periodicity $2π$, i.e.\n", "$$\n", " \\sin(x) = \\sin(x + 2πm) \\quad ∀ m ∈ \\mathbb{Z},\n", "$$\n", "This is nothing else than saying that any translation by an integer multiple of $2π$\n", "keeps the $\\sin$ function invariant. In a more formal way one can use the\n", "translation operator $T_{-2πm}$ to write this as:\n", "$$\n", " T_{-2πm} \\, \\sin(x) = \\sin(x + 2πm) = \\sin(x).\n", "$$\n", "\n", "Whenever such periodicity exists one can exploit it to save computational work.\n", "Consider a problem in which we want to find a function $f : \\mathbb{R} → \\mathbb{R}$,\n", "but *a priori* the solution is known to be periodic with periodicity $a$. As a consequence\n", "of said periodicity it is sufficient to determine the values of $f$ for all $x$ from the\n", "interval $[-a/2, a/2)$ to uniquely define $f$ on the full real axis. Naturally exploiting\n", "periodicity in a computational procedure thus greatly reduces the required amount of work.\n", "\n", "Let us introduce some jargon: The periodicity of our problem implies that we may define\n", "a **lattice**\n", "```\n", " -3a/2 -a/2 +a/2 +3a/2\n", " ... |---------|---------|---------| ...\n", " a a a\n", "```\n", "with lattice constant $a$. Each cell of the lattice is an identical periodic image of\n", "any of its neighbors. For finding $f$ it is thus sufficient to consider only the\n", "problem inside a **unit cell** $[-a/2, a/2)$. In passing we note that the definition\n", "of the unit cell is itself only unique up to translations. A choice $[0, a)$,\n", "for example, would have done just as well." ] }, { "cell_type": "markdown", "id": "ccf8bda9", "metadata": {}, "source": [ "## Periodic operators and the Bloch transform\n", "Not only functions, but also operators can feature periodicity.\n", "Consider for example the **free-electron Hamiltonian**\n", "$$\n", " H = -\\frac12 Δ.\n", "$$\n", "In free-electron model (which gives rise to this Hamiltonian) electron motion is only\n", "by their own kinetic energy. As this model features no potential which could make one point\n", "in space more preferred than another, we would expect this model to be periodic.\n", "If an operator is periodic with respect to a lattice such as the one defined above,\n", "than it commutes with all lattice translations. For the free-electron case $H$\n", "one can easily show exactly that, i.e.\n", "$$\n", " T_{ma} H = H T_{ma} \\quad ∀ m ∈ \\mathbb{Z}.\n", "$$\n", "We note in passing that the free-electron model is actually very special in the sense that\n", "the choice of $a$ is completely arbitrary here. In other words $H$ is periodic\n", "with respect to any translation. In general, however, periodicity is only\n", "attained with respect to a finite number of translations $a$ and we will take this\n", "viewpoint here.\n", "\n", "**Bloch's theorem** now tells us that for periodic operators,\n", "the solutions to the eigenproblem\n", "$$\n", " H ψ_{kn} = ε_{kn} ψ_{kn}\n", "$$\n", "satisfy a factorization\n", "$$\n", " ψ_{kn}(x) = e^{i k⋅x} u_{kn}(x)\n", "$$\n", "into a plane wave $e^{i k⋅x}$ and a lattice-periodic function\n", "$$\n", " T_{ma} u_{kn}(x) = u_{kn}(x - ma) = u_{kn}(x) \\quad ∀ m ∈ \\mathbb{Z}.\n", "$$\n", "In this $n$ is a labeling integer index and $k$ is a real number,\n", "whose details will be clarified in the next section.\n", "The index $n$ is sometimes also called the **band index** and\n", "functions $ψ_{kn}$ satisfying this factorization are also known as\n", "**Bloch functions** or **Bloch states**.\n", "\n", "Consider the application of $2H = -Δ = - \\frac{d^2}{d x^2}$\n", "to such a Bloch wave. First we notice for any function $f$\n", "$$\n", " -i∇ \\left( e^{i k⋅x} f \\right)\n", " = -i\\frac{d}{dx} \\left( e^{i k⋅x} f \\right)\n", " = k e^{i k⋅x} f + e^{i k⋅x} (-i∇) f = e^{i k⋅x} (-i∇ + k) f.\n", "$$\n", "Using this result twice one shows that applying $-Δ$ yields\n", "$$\n", "\\begin{aligned}\n", " -\\Delta \\left(e^{i k⋅x} u_{kn}(x)\\right)\n", " &= -i∇ ⋅ \\left[-i∇ \\left(u_{kn}(x) e^{i k⋅x} \\right) \\right] \\\\\n", " &= -i∇ ⋅ \\left[e^{i k⋅x} (-i∇ + k) u_{kn}(x) \\right] \\\\\n", " &= e^{i k⋅x} (-i∇ + k)^2 u_{kn}(x) \\\\\n", " &= e^{i k⋅x} 2H_k u_{kn}(x),\n", "\\end{aligned}\n", "$$\n", "where we defined\n", "$$\n", " H_k = \\frac12 (-i∇ + k)^2.\n", "$$\n", "The action of this operator on a function $u_{kn}$ is given by\n", "$$\n", " H_k u_{kn} = e^{-i k⋅x} H e^{i k⋅x} u_{kn},\n", "$$\n", "which in particular implies that\n", "$$\n", " H_k u_{kn} = ε_{kn} u_{kn} \\quad ⇔ \\quad H (e^{i k⋅x} u_{kn}) = ε_{kn} (e^{i k⋅x} u_{kn}).\n", "$$\n", "To seek the eigenpairs of $H$ we may thus equivalently\n", "find the eigenpairs of *all* $H_k$.\n", "The point of this is that the eigenfunctions $u_{kn}$ of $H_k$\n", "are periodic (unlike the eigenfunctions $ψ_{kn}$ of $H$).\n", "In contrast to $ψ_{kn}$ the functions $u_{kn}$ can thus be fully\n", "represented considering the eigenproblem only on the unit cell.\n", "\n", "A detailed mathematical analysis shows that the transformation from $H$\n", "to the set of all $H_k$ for a suitable set of values for $k$ (details below)\n", "is actually a unitary transformation, the so-called **Bloch transform**.\n", "This transform brings the Hamiltonian into the symmetry-adapted basis for\n", "translational symmetry, which are exactly the Bloch functions.\n", "Similar to the case of choosing a symmetry-adapted basis for other kinds of symmetries\n", "(like the point group symmetry in molecules), the Bloch transform also makes\n", "the Hamiltonian $H$ block-diagonal:\n", "$$\n", " T_B H T_B^{-1} ⟶ \\left( \\begin{array}{cccc} H_1&&&0 \\\\ &H_2\\\\&&H_3\\\\0&&&\\ddots \\end{array} \\right)\n", "$$\n", "with each block $H_k$ taking care of one value $k$.\n", "This block-diagonal structure under the basis of Bloch functions lets us\n", "completely describe the spectrum of $H$ by looking only at the spectrum\n", "of all $H_k$ blocks.\n", "\n", "**Aside:** Block-diagonal is a bit an abuse of terms here, since the Hamiltonian\n", " is not a matrix but an operator and the number of blocks is essentially infinite.\n", " The mathematically precise term is that the Bloch transform reveals the fibers\n", " of the Hamiltonian." ] }, { "cell_type": "markdown", "id": "4ab5317d", "metadata": {}, "source": [ "## The Brillouin zone\n", "\n", "We now consider the parameter $k$ of the Hamiltonian blocks in detail.\n", "\n", "- As discussed $k$ is a real number. It turns out, however, that some of\n", " these $k∈\\mathbb{R}$ give rise to operators related by unitary transformations\n", " (again due to translational symmetry).\n", "- Since such operators have the same eigenspectrum, only one version needs to be considered.\n", "- The smallest subset from which $k$ is chosen is the **Brillouin zone** (BZ).\n", "\n", "- The BZ is the unit cell of the **reciprocal lattice**, which may be constructed from\n", " the **real-space lattice** by a Fourier transform.\n", "- In our simple 1D case the reciprocal lattice is just\n", " ```\n", " ... |--------|--------|--------| ...\n", " 2π/a 2π/a 2π/a\n", " ```\n", " i.e. like the real-space lattice, but just with a different lattice constant\n", " $b = 2π / a$.\n", "- The BZ in our example is thus $B = [-π/a, π/a)$. The members of $B$\n", " are typically called $k$-points." ] }, { "cell_type": "markdown", "id": "27da28eb", "metadata": {}, "source": [ "## Discretization and plane-wave basis sets\n", "\n", "With what we discussed so far the strategy to find all eigenpairs of a periodic\n", "Hamiltonian $H$ thus reduces to finding the eigenpairs of all $H_k$ with $k ∈ B$.\n", "This requires *two* discretisations:\n", "\n", " - $B$ is infinite (and not countable). To discretize we first only pick a finite number\n", " of $k$-points. Usually this **$k$-point sampling** is done by picking $k$-points\n", " along a regular grid inside the BZ, the **$k$-grid**.\n", " - Each $H_k$ is still an infinite-dimensional operator.\n", " Following a standard Ritz-Galerkin ansatz we project the operator into a finite basis\n", " and diagonalize the resulting matrix.\n", "\n", "For the second step multiple types of bases are used in practice (finite differences,\n", "finite elements, Gaussians, ...). We will conider here so-called plane-wave\n", "discretizations.\n", "\n", "For our 1D example normalized plane waves are defined as the functions\n", "$$\n", "e_{G}(x) = \\frac{e^{i G x}}{\\sqrt{a}} \\qquad G \\in b\\mathbb{Z}\n", "$$\n", "and typically one forms basis sets from these by specifying a\n", "**kinetic energy cutoff** $E_\\text{cut}$:\n", "$$\n", "\\left\\{ e_{G} \\, \\big| \\, (G + k)^2 \\leq 2E_\\text{cut} \\right\\}\n", "$$" ] }, { "cell_type": "markdown", "id": "swedish-sherman", "metadata": {}, "source": [ "## Let's crunch some numbers!\n", "\n", "One typical approach to get physical insight into a Hamiltonian $H$ is to plot a so-called **band structure**, that is the eigenvalues of $H_k$ versus $k$. The rough steps to plot a band structure are:\n", "\n", "1. Select a set of $k$-points ($k$-point sampling)\n", "1. Select plane-wave cutoff, thus the basis for discretising $H_k$\n", "1. Build and diagonalise $H_k$ for each $k$.\n", "1. Plot eigenvalues versus $k$\n", "\n", "Notice that the free-electrion $H_k = \\frac12 (-i∇ + k)^2$ is already diagonal within the plane-wave basis\n", "\\begin{align}\n", "\\langle e_{G} | H_k e_{G'} \\rangle\n", "&= \\frac12 \\left\\langle e_{G} \\middle| (-i \\nabla + k)^2 \\, e_{G'} \\right\\rangle \\\\\n", "&= \\frac12 (G + k)^2 \\left\\langle e_{G} \\middle|e_{G'}\\right\\rangle \\\\\n", "&= \\frac12 (G + k)^2 \\delta_{GG'}.\n", "\\end{align}\n", "\n", "Therefore obtaining such a band structure from scratch can be done in only a few lines of code:" ] }, { "cell_type": "code", "execution_count": null, "id": "incorporated-earthquake", "metadata": {}, "outputs": [], "source": [ "using LinearAlgebra\n", "using Plots\n", "\n", "# Parameters\n", "a = 100 # Lattice constant\n", "Ecut = 300 # Cutoff in Hartree\n", "kgrid = 50 # Number of points on equally-spaced grid for k\n", "\n", "# Derived quantities\n", "b = 2π / a # Reciprocal lattice\n", "\n", "# Step 1: k-Points\n", "kpoints = b * (collect(1:kgrid) .- ceil(Int, kgrid / 2)) ./ kgrid\n", "\n", "# Step 2: Basis for H_k\n", "# Represented as one array of all valid G*b per kpoint\n", "Gmax = ceil(Int, sqrt(2Ecut) + b) # Rough upper bound for G\n", "Gs = [[Gidx*b for Gidx in -Gmax:Gmax if abs2(Gidx*b + k) ≤ 2Ecut]\n", " for k in kpoints]\n", "\n", "# Step 3: Build the discretised Hk. In this case it is diagonal,\n", "# i.e. its diagonal values (== eigenvalues) are all we need.\n", "# We directly determine them and sort them ascendingly\n", "ev_Hk = [sort([abs2(G + k)/2 for G in Gs[ik]])\n", " for (ik, k) in enumerate(kpoints)]\n", "\n", "# Step 4: Plot the bands\n", "n_bands = 6\n", "bands = [[ev_Hk[ik][iband] for ik in 1:length(kpoints)]\n", " for iband in 1:n_bands]\n", "\n", "p = plot()\n", "for iband in 1:n_bands\n", " plot!(p, kpoints, bands[iband], color=:blue, label=\"\")\n", "end\n", "p" ] }, { "cell_type": "markdown", "id": "80cac9c5", "metadata": {}, "source": [ "## Introducing the density-functional toolkit (DFTK)" ] }, { "cell_type": "markdown", "id": "e4d3ff6a", "metadata": {}, "source": [ "For such a simple problem coding it up by hand is arguably still managable. But considering that we not only want to work on a free-electron problem, but consider density-functional theory on 3D materials, the story is totally different. Still, being able to control and access computations and intermediate quantities to explore the physics and design novel algorithms is a desirable criterion in mathematically-oriented research.\n", "\n", "For this reason in 2019 we started the **density-functional toolkit** (DFTK, https://dftk.org). A few key characteristics:" ] }, { "cell_type": "markdown", "id": "6edbc1bb", "metadata": {}, "source": [ "\n", "\n", "- Lowers barriers for interdisciplinary research:\n", " * Restriction to relevant model problems\n", " * Scale-up to application regime (1000 electrons)\n", " * Support for parallelisation, GPU (planned), ...\n", "- 7000 lines of julia code, open-source\n", " * Low entrance barrier\n", " * Scripting *and* algorithm development\n", " * Only one language required!" ] }, { "cell_type": "markdown", "id": "36c1761d", "metadata": {}, "source": [ "Returning to our free-electron example from before, what we coded up in about 30 lines of code, can be done in just about 5 lines using DFTK:" ] }, { "cell_type": "markdown", "id": "68d316ec", "metadata": {}, "source": [ "**Step 1:** Build the 1D lattice. DFTK is mostly tailored for 3D problems.\n", "Therefore quantities related to the problem space are have a fixed\n", "dimension 3. The convention is that for 1D / 2D problems the\n", "trailing entries are always zero and ignored in the computation.\n", "For the lattice we therefore construct a 3x3 matrix with only one entry." ] }, { "cell_type": "code", "execution_count": null, "id": "db17cc85", "metadata": {}, "outputs": [], "source": [ "using DFTK\n", "\n", "lattice = zeros(3, 3)\n", "lattice[1, 1] = 20.;" ] }, { "cell_type": "markdown", "id": "b4b22be8", "metadata": {}, "source": [ "**Step 2:** Select a model. In this case we choose again the free-electron model,\n", "which is the same as saying that there is only a Kinetic term\n", "(and no potential) in the model." ] }, { "cell_type": "code", "execution_count": null, "id": "bec57887", "metadata": {}, "outputs": [], "source": [ "model = Model(lattice; terms=[Kinetic()])" ] }, { "cell_type": "markdown", "id": "f3bab17d", "metadata": {}, "source": [ "**Step 3:** Define a plane-wave basis using this model and a cutoff $E_\\text{cut}$\n", "of 300 Hartree. The $k$-point grid is given as a regular grid in the BZ\n", "(a so-called **Monkhorst-Pack** grid). Here we select only one $k$-point (1x1x1),\n", "see the note below for some details on this choice." ] }, { "cell_type": "code", "execution_count": null, "id": "d11fa193", "metadata": {}, "outputs": [], "source": [ "basis = PlaneWaveBasis(model; Ecut=300, kgrid=(1, 1, 1))" ] }, { "cell_type": "markdown", "id": "bbf512e6", "metadata": {}, "source": [ "**Step 4:** Plot the bands! Select a density of $k$-points for the $k$-grid to use\n", "for the bandstructure calculation, discretize the problem and diagonalize it.\n", "Afterwards plot the bands." ] }, { "cell_type": "code", "execution_count": null, "id": "b572ece1", "metadata": {}, "outputs": [], "source": [ "using Unitful\n", "using UnitfulAtomic\n", "using Plots\n", "\n", "plot_bandstructure(basis; n_bands=6, kline_density=100)" ] }, { "cell_type": "markdown", "id": "74fd020b", "metadata": {}, "source": [ "**Note:** Selection of k-point grids in `PlaneWaveBasis` construction \n", " You might wonder why we only selected a single $k$-point (clearly a very crude\n", " and inaccurate approximation). In this example the `kgrid` parameter specified\n", " in the construction of the `PlaneWaveBasis`\n", " is not actually used for plotting the bands. It is only used when solving more\n", " involved models like density-functional theory (DFT) where the Hamiltonian is\n", " non-linear. In these cases before plotting the bands the self-consistent field\n", " equations (SCF) need to be solved first. This is typically done on\n", " a different $k$-point grid than the grid used for the bands later on.\n", " In our case we don't need this extra step and therefore the `kgrid` value passed\n", " to `PlaneWaveBasis` is actually arbitrary." ] }, { "cell_type": "markdown", "id": "e0e1532f", "metadata": {}, "source": [ "## Adding potentials\n", "So far so good. But free electrons are actually a little boring,\n", "so let's add a potential interacting with the electrons.\n", "\n", "- The modified problem we will look at consists of diagonalizing\n", " $$\n", " H_k = \\frac12 (-i \\nabla + k)^2 + V\n", " $$\n", " for all $k \\in B$ with a periodic potential $V$ interacting with the electrons.\n", "\n", "- A number of \"standard\" potentials are readily implemented in DFTK and\n", " can be assembled using the `terms` kwarg of the model.\n", " This allows to seamlessly construct\n", "\n", " * density-functial theory (DFT) models for treating electronic structures\n", " (see the Tutorial).\n", " * Gross-Pitaevskii models for bosonic systems\n", " (see Gross-Pitaevskii equation in one dimension)\n", " * even some more unusual cases like anyonic models.\n", "\n", "In this tutorial we will go a little more low-level and directly provide\n", "an analytic potential describing the interaction with the electrons to DFTK.\n", "\n", "First we define a custom Gaussian potential as a new \"element\" inside DFTK:" ] }, { "cell_type": "code", "execution_count": null, "id": "6d151663", "metadata": {}, "outputs": [], "source": [ "struct ElementGaussian <: DFTK.Element\n", " α # Prefactor\n", " L # Extend\n", "end\n", "\n", "# Some default values\n", "ElementGaussian() = ElementGaussian(0.3, 10.0)\n", "\n", "# Real-space representation of a Gaussian\n", "function DFTK.local_potential_real(el::ElementGaussian, r::Real)\n", " -el.α / (√(2π) * el.L) * exp(- (r / el.L)^2 / 2)\n", "end\n", "\n", "# Fourier-space representation of the Gaussian\n", "function DFTK.local_potential_fourier(el::ElementGaussian, q::Real)\n", " # = ∫ -α exp(-(r/L)^2 exp(-ir⋅q) dr\n", " -el.α * exp(- (q * el.L)^2 / 2)\n", "end" ] }, { "cell_type": "markdown", "id": "5ab11b64", "metadata": {}, "source": [ "A single potential looks like:" ] }, { "cell_type": "code", "execution_count": null, "id": "b4c2f9fa", "metadata": {}, "outputs": [], "source": [ "using Plots\n", "using LinearAlgebra\n", "nucleus = ElementGaussian()\n", "plot(r -> DFTK.local_potential_real(nucleus, norm(r)), xlims=(-50, 50))" ] }, { "cell_type": "markdown", "id": "fa60dc41", "metadata": {}, "source": [ "With this element at hand we can easily construct a setting\n", "where two potentials of this form are located at positions\n", "``20`` and ``80`` inside the lattice ``[0, 100]``:" ] }, { "cell_type": "code", "execution_count": null, "id": "3259f7f0", "metadata": {}, "outputs": [], "source": [ "using LinearAlgebra\n", "\n", "# Define the 1D lattice [0, 100]\n", "lattice = diagm([100., 0, 0])\n", "\n", "# Place them at 20 and 80 in *fractional coordinates*,\n", "# that is 0.2 and 0.8, since the lattice is 100 wide.\n", "nucleus = ElementGaussian()\n", "atoms = [nucleus, nucleus]\n", "positions = [[0.2, 0, 0], [0.8, 0, 0]]\n", "\n", "# Assemble the model, discretize and build the Hamiltonian\n", "model = Model(lattice, atoms, positions; terms=[Kinetic(), AtomicLocal()])\n", "basis = PlaneWaveBasis(model; Ecut=300, kgrid=(1, 1, 1));\n", "ham = Hamiltonian(basis)\n", "\n", "# Extract the total potential term of the Hamiltonian and plot it\n", "potential = DFTK.total_local_potential(ham)[:, 1, 1]\n", "rvecs = collect(r_vectors_cart(basis))[:, 1, 1] # slice along the x axis\n", "x = [r[1] for r in rvecs] # only keep the x coordinate\n", "plot(x, potential, label=\"\", xlabel=\"x\", ylabel=\"V(x)\")" ] }, { "cell_type": "markdown", "id": "88862f56", "metadata": {}, "source": [ "This potential is the sum of two \"atomic\" potentials (the two \"Gaussian\" elements).\n", "Due to the periodic setting we are considering interactions naturally also occur\n", "across the unit cell boundary (i.e. wrapping from `100` over to `0`).\n", "The required periodization of the atomic potential is automatically taken care,\n", "such that the potential is smooth across the cell boundary at `100`/`0`.\n", "\n", "With this setup, let's look at the bands:" ] }, { "cell_type": "code", "execution_count": null, "id": "78219d38", "metadata": {}, "outputs": [], "source": [ "using Unitful\n", "using UnitfulAtomic\n", "\n", "plot_bandstructure(basis; n_bands=6, kline_density=500)" ] }, { "cell_type": "markdown", "id": "448afb5a", "metadata": {}, "source": [ "The bands are noticeably different.\n", " - The bands no longer overlap, meaning that the spectrum of $H$ is no longer continuous\n", " but has gaps.\n", "\n", " - The two lowest bands are almost flat. This is because they represent\n", " two tightly bound and localized electrons inside the two Gaussians.\n", "\n", " - The higher the bands are in energy, the more free-electron-like they are.\n", " In other words the higher the kinetic energy of the electrons, the less they feel\n", " the effect of the two Gaussian potentials. As it turns out the curvature of the bands,\n", " (the degree to which they are free-electron-like) is highly related to the delocalization\n", " of electrons in these bands: The more curved the more delocalized. In some sense\n", " \"free electrons\" correspond to perfect delocalization.\n", " \n", "**Live exercise 2.1:** Try playing with the parameters of the Gaussian potentials by setting\n", "```julia\n", "nucleus = ElementGaussian(α, L)\n", "```\n", "with different $\\alpha$ and $L$ in the above procedure." ] }, { "cell_type": "markdown", "id": "economic-wagon", "metadata": {}, "source": [ "#### Takeaways\n", "\n", "- For periodic problems the Bloch transform allows to find the eigenpairs of an operator $k$-Point by $k$-Point\n", "- The $k$ points are taken in a discrete mesh from the Brillouin zone.\n", "- A common way to visualise the eigenvalues is as a plot versus $k$, the band structure." ] } ], "metadata": { "kernelspec": { "display_name": "Julia 1.7.3", "language": "julia", "name": "julia-1.7" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.7.3" } }, "nbformat": 4, "nbformat_minor": 5 }