Silicon's band structure, from the crystal up

A notebook. Each cell's output appears under it; press Run all, or run a cell with Shift+Enter. The code is a small version of the plane-wave DFT program in dft/ (the dense engine, dft/pwdft.viv), cut into the order the physics goes in. Numbers are in hartree atomic units: lengths in bohr, energies in hartree (1 Ha = 27.2114 eV) until the last cells.

The whole calculation solves one equation, the Kohn-Sham equation above: each electron moves in the potential of the ions, plus the classical repulsion of the electron density $\rho$ (Hartree), plus an exchange-correlation term for what the classical repulsion misses. The potential depends on $\rho$, and $\rho$ is built from the solutions, so it is solved by iteration until the density stops changing.

The crystal

The cell above sets up silicon's crystal. Silicon has the diamond structure: an fcc lattice with lattice constant $a=10.26$ bohr (5.43 Å), and two atoms in each primitive cell, at $\pm \frac{a}{8}(1,1,1)$. The rows of A are the three primitive vectors; the cell they span holds a quarter of the cube, $a^3/4=270.01$ bohr³.

The rows of B are the reciprocal vectors, defined by ${\mathbf{b}}_i\cdot {\mathbf{a}}_j=2\pi {\delta }_{ij}$. Every function with the period of the crystal is a sum of plane waves $e^{i\mathbf{G}\cdot \mathbf{r}}$ with $\mathbf{G}=m_1{\mathbf{b}}_1+m_2{\mathbf{b}}_2+m_3{\mathbf{b}}_3$ for integers $m_i$. That is the reason to work with plane waves: the potential, the density and the periodic part of every orbital are Fourier series on this lattice.

The next cell describes one silicon ion to the valence electrons. A silicon atom has 14 electrons; the 10 core electrons stay put when atoms bond, so they are replaced, together with the nucleus, by a pseudopotential acting on the 4 valence electrons. The parameters are the Goedecker-Teter-Hutter ones for silicon, copied from CP2K's table by dft/pseudo/make_pseudo.py. Its local part is a function of $|\mathbf{G}|$ only. Expect it to be negative (attractive) and to vanish at large $|\mathbf{G}|$, since the ion looks smooth to the valence electrons.

The pseudopotential, local part

The first term, $-4\pi Z{e}^{-x^2/2}/(\Omega g^2)$, is the Coulomb attraction of an ion of charge $Z=4$ whose charge is spread over a Gaussian of radius $r_{\text{loc}}=0.44$ bohr; the second is a short-range correction. At $g=0$ the Coulomb term diverges. So do the electrons' own repulsion and the repulsion between ions, and in a neutral crystal the three cancel, so the function returns only the finite remainder there.

The local part is the same for every electron. The nonlocal part is not: it acts differently on the $s$ ($l=0$) and $p$ ($l=1$) parts of a wave around each ion, which is how it imitates the core electrons that were removed. It is built from projectors, Gaussians times polynomials in $|\mathbf{k}+\mathbf{G}|$, two for $l=0$ and one for $l=1$, coupled by the small matrices HS and HP. The next cell defines them and evaluates each at three values of $|\mathbf{k}+\mathbf{G}|$.

The nonlocal part, and the k-points

The $p$ projector starts at zero and grows with $|\mathbf{k}+\mathbf{G}|$: a plane wave with small $|\mathbf{k}+\mathbf{G}|$ is nearly constant across the ion, and a constant has no $p$ ($l=1$) component.

An electron in a crystal is a Bloch wave, $e^{i\mathbf{k}\cdot \mathbf{r}}$ times a periodic function, and each wave vector $\mathbf{k}$ in the Brillouin zone gives its own set of energy bands. The density averages over $\mathbf{k}$, and that average is replaced by a sum over a few points. The next cell uses the k-mesh module of the full program, dft/kmesh.viv, at its smallest setting. It is a four-shift mesh: its points are closed under the crystal's rotations, so the density comes out with the crystal's symmetry without any explicit symmetrization. Expect four points, the four $L$ points $(\frac{1}{2},\frac{1}{2},\frac{1}{2})$, $(\frac{1}{2},0,0)$, $(0,\frac{1}{2},0)$ and $(0,0,\frac{1}{2})$ in units of the reciprocal vectors, each with weight $\frac{1}{4}$.

The plane-wave basis

The basis at one $\mathbf{k}$ is every plane wave whose kinetic energy $\frac{1}{2}|\mathbf{k}+\mathbf{G}{|}^2$ is at most the cutoff. The full program uses 10 Ha for silicon; this notebook uses 5 Ha so it runs in a few seconds. The cutoff is the basis's one tunable parameter (the k-point mesh is the other convergence setting): raising it adds shorter waves and lowers the energy.

plane_waves loops over a box of integer triples $m$ large enough to hold every $\mathbf{G}$ that can pass, and keeps the ones that do. Called with slot $-1$ it only counts, and records the largest $|m_i|$ in MAXM. The count is done at the four k-points and at the 21 points of the band path used later, from $L$ through $\Gamma$ to $X$, so the arrays and the FFT grid are big enough for all of them. Expect 137 plane waves at $\Gamma$ and a few more at the other points.

The cell above only defines things. The next one counts, sizes the arrays (one slot per k-point, and one spare slot for the band path), and fills each k-point's slot.

The FFT grid and the ions' potential

The Hamiltonian couples two plane waves $\mathbf{G}$ and ${\mathbf{G}}^{\prime }$ through the potential's Fourier component at $\mathbf{G}-{\mathbf{G}}^{\prime }$, so the grid must hold differences up to twice the largest $|m_i|$ without wrapping around. Any side of at least $4\max |m_i|+2$ points does; fft_nice rounds up to a size with no prime factor above 5, which the FFT handles fastest.

On that grid the cell tabulates $1/|\mathbf{G}{|}^2$ (for the Hartree term later) and the ions' potential. Each atom contributes its local pseudopotential times a phase $e^{-i\mathbf{G}\cdot \mathbf{\tau }}$, the structure factor, which says where in the cell the atom sits. Silicon's two atoms sit at $\pm \mathbf{\tau }$, so the phases add to $2\cos (\mathbf{G}\cdot \mathbf{\tau })$, a real number: the ions' potential has no imaginary part, which the cell checks. Expect a grid of 18 points on a side.

Exchange and correlation

Two electrons repel, so each keeps the others at a distance, and electrons of the same spin are kept apart further by the Pauli principle. The Hartree term treats the density as a smooth charge cloud and misses both. The local density approximation (LDA) adds, at each point, the exchange-correlation energy per electron ${\epsilon }_{xc}(\rho )$ of a uniform electron gas with the density found there. The cell uses the Teter-93 Padé fit, a ratio of two polynomials in the Wigner-Seitz radius $r_s=(3/4\pi \rho {)}^{1/3}$, the radius of a sphere that holds one electron. It is the functional the pseudopotential was fitted with. The potential is the derivative, $v_{xc}=d(\rho {\epsilon }_{xc})/d\rho ={\epsilon }_{xc}-\frac{r_s}{3}d{\epsilon }_{xc}/d{r}_s$.

Expect both to be negative and to grow in size with the density: silicon's valence density is about 0.03 electrons per bohr³.

The Hamiltonian at one k-point

In the plane-wave basis at $\mathbf{k}$, the Hamiltonian is a matrix with one row per plane wave:

The kinetic energy is diagonal. The local potential $V$ couples two plane waves through its Fourier component at their difference, read off the grid. The nonlocal part is a sum over projectors, $V_{\text{nl}}={\sum }_{IJ}|{\beta }_I⟩D_{IJ}⟨{\beta }_J|$, where each ${\beta }_I$ is one projector of one atom times an angular factor (a real spherical harmonic of the direction of $\mathbf{k}+\mathbf{G}$) times the atom's structure factor. Five per atom: two $s$ projectors and one $p$ projector with three orientations. D holds the coupling matrices.

Every k-point has its own number of plane waves, so the matrices are padded to NPWMAX rows. The padding rows are decoupled and given diagonal values near 1000 Ha, far above any band, so they never mix in.

To see the Hamiltonian before any self-consistency, the cell diagonalizes it at $\Gamma$ in the starting potential: the ions plus the exchange-correlation potential of a uniform density (whose Hartree potential is zero). Expect one low $s$-like state and then a triple degeneracy: at $\Gamma$ the three $p$ directions are equivalent.

The next cell builds $H$ itself. pair_index gives, for every pair of plane waves, where their difference $\mathbf{G}-{\mathbf{G}}^{\prime }$ sits in the FFT grid, using whole-matrix arithmetic: an outer difference of the integer triples, wrapped onto the grid (a negative component $m$ is stored at $m+N$). wrap does the same for one number. take then reads the potential at all those places at once.

Now the starting potential and one diagonalization at $\Gamma$, in the spare slot. potential assembles the total local potential on the grid from a density: the ions, the Hartree potential and exchange-correlation. The Hartree potential comes from Poisson's equation, ${\nabla }^2{V}_H=-4\pi \rho$, which in reciprocal space is algebra: $V_H(\mathbf{G})=4\pi \rho (\mathbf{G})/|\mathbf{G}{|}^2$. An FFT takes the density to reciprocal space, one multiplication solves Poisson's equation, and the $\mathbf{G}=0$ term is left out, as the ions' $g=0$ term was.

The self-consistent loop

Silicon has 8 valence electrons per cell, two in each of the lowest 4 bands at every k-point. One pass of the loop:

1. build the potential from the density RHO;
2. at each k-point, build $H$, diagonalize it, and add the lowest 4 bands' density $2|\psi (\mathbf{r}){|}^2$ (an inverse FFT takes each band to the grid) with the k-point's weight;
3. compute the total energy and how far the new density is from the old one;
4. if they have both stopped changing (below $10^{-8}$), stop; otherwise mix: move the density 60% of the way toward the new one.

The energy is the sum of the occupied band energies, corrected for double counting: each band energy contains the electron's interaction with the whole density, which counts every pair of electrons twice, so the Hartree energy is subtracted once and $\int \rho v_{xc}$ is replaced by $\int \rho {\epsilon }_{xc}$.

The full program mixes with the last three iterations (Anderson mixing, dft/pw_common.viv), which takes fewer iterations; plain mixing is one line and reaches the same density. Plain mixing converges geometrically: once the first few iterations are past, expect the density change to shrink by about the same factor at every iteration.

The loop itself.

In the loop's output the density change fell by a factor of 0.42 per iteration, from 0.46 to $5\times 10^{-9}$ in 21 iterations.

The bands along a path

With the potential converged, the bands at any $\mathbf{k}$ come from one diagonalization each, with no further iteration. The path runs from $L$, the center of a hexagonal face of the Brillouin zone, to its center $\Gamma$, and on to $X$, the center of a square face. Energies are shifted so that the top of the valence bands is zero, and converted to eV.

What to look for: the valence-band top is at $\Gamma$ and threefold degenerate. The lowest conduction state is not at $\Gamma$ but most of the way to $X$, so silicon's gap is indirect: an electron at the bottom of the conduction band has a different crystal momentum from a hole at the top of the valence band. The LDA gap is also too small, about 0.5 eV where the measured gap is 1.17 eV. That is a known failure of the LDA, not of the code.

The next cell finds the valence-band top and the conduction-band bottom along the path and prints the six lowest bands at $L$, $\Gamma$ and $X$ in eV.

How far from converged

The full program runs silicon at 10 Ha with 16 k-points (dft/si.viv), and its check file pins the band energies at $L$, $\Gamma$ and $X$. The next cell compares against those: the largest difference among the valence bands, and the difference for the lowest conduction state at $L$. Two settings differ, the cutoff and the k-points, so each difference mixes the two effects. The energies on both sides are measured from the valence-band top, so the comparison is of spacings between bands, and the differences need not have one sign. Expect a few tenths of an eV either way: the lowest valence band higher here than at 10 Ha, the second pair of valence states at $X$ lower.

The picture

The six lowest bands along the path. The four at or below zero are the valence bands, which hold silicon's 8 valence electrons; the two above are the lowest conduction bands. The 10 Ha run has the same degeneracies at $L$, $\Gamma$ and $X$ and the same order of states; the cell above measured how far the energies differ.