This site is version 0.0.1. Most of the code on this site is drafted by an agent and needs review. I apologize for the AI-speak, we are working on rewriting it!
viv-lang.org

Notebooks · posts/silicon-bands-from-the-crystal-up

Silicon's band structure, from the crystal up

Adam Zettel · September 25, 2026 · written with an agent

A density-functional calculation of silicon, one cell at a time: the crystal, a plane-wave basis, the pseudopotential, the Hamiltonian, the self-consistent loop and the bands from L through Gamma to X. It runs at a 5 Ha cutoff on four k-points so it fits in a browser tab; its tests check each stage, from the basis sizes to the band energies, against an independent numpy implementation at the same settings.

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 (Hartree), plus an exchange-correlation term for what the classical repulsion misses. The potential depends on , and is built from the solutions, so it is solved by iteration until the density stops changing.

Look inside this program

  • Atlas

    The program’s map: every function and which calls which.

  • Explore

    Run it and step through what each part did, with the code, values and plots beside it.

Shift+Enter in a cell runs it and moves to the next. Running a cell runs every cell above it too, as one program.

[1]lines 25–68

The crystal

The cell above sets up silicon's crystal. Silicon has the diamond structure: an fcc lattice with lattice constant bohr (5.43 Å), and two atoms in each primitive cell, at . The rows of A are the three primitive vectors; the cell they span holds a quarter of the cube, bohr³.

The rows of B are the reciprocal vectors, defined by . Every function with the period of the crystal is a sum of plane waves with for integers . 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 only. Expect it to be negative (attractive) and to vanish at large , since the ion looks smooth to the valence electrons.

[2]lines 97–124

The pseudopotential, local part

The first term, , is the Coulomb attraction of an ion of charge whose charge is spread over a Gaussian of radius bohr; the second is a short-range correction. At 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 () and () 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 , two for and one for , coupled by the small matrices HS and HP. The next cell defines them and evaluates each at three values of .

[3]lines 145–180

The nonlocal part, and the k-points

The projector starts at zero and grows with : a plane wave with small is nearly constant across the ion, and a constant has no () component.

An electron in a crystal is a Bloch wave, times a periodic function, and each wave vector in the Brillouin zone gives its own set of energy bands. The density averages over , 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 points , , and in units of the reciprocal vectors, each with weight .

[4]lines 202–212

The plane-wave basis

The basis at one is every plane wave whose kinetic energy 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 large enough to hold every that can pass, and keeps the ones that do. Called with slot it only counts, and records the largest in MAXM. The count is done at the four k-points and at the 21 points of the band path used later, from through to , so the arrays and the FFT grid are big enough for all of them. Expect 137 plane waves at and a few more at the other points.

[5]lines 232–280

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.

[6]lines 288–310

The FFT grid and the ions' potential

The Hamiltonian couples two plane waves and through the potential's Fourier component at , so the grid must hold differences up to twice the largest without wrapping around. Any side of at least 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 (for the Hartree term later) and the ions' potential. Each atom contributes its local pseudopotential times a phase , the structure factor, which says where in the cell the atom sits. Silicon's two atoms sit at , so the phases add to , a real number: the ions' potential has no imaginary part, which the cell checks. Expect a grid of 18 points on a side.

[7]lines 330–382

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 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 , the radius of a sphere that holds one electron. It is the functional the pseudopotential was fitted with. The potential is the derivative, .

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

[8]lines 402–432

The Hamiltonian at one k-point

In the plane-wave basis at , the Hamiltonian is a matrix with one row per plane wave:

The kinetic energy is diagonal. The local potential couples two plane waves through its Fourier component at their difference, read off the grid. The nonlocal part is a sum over projectors, , where each is one projector of one atom times an angular factor (a real spherical harmonic of the direction of ) times the atom's structure factor. Five per atom: two projectors and one 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 in the starting potential: the ions plus the exchange-correlation potential of a uniform density (whose Hartree potential is zero). Expect one low -like state and then a triple degeneracy: at the three directions are equivalent.

[9]lines 462–505

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

[10]lines 516–558

Now the starting potential and one diagonalization at , 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, , which in reciprocal space is algebra: . An FFT takes the density to reciprocal space, one multiplication solves Poisson's equation, and the term is left out, as the ions' term was.

[11]lines 571–611

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 , diagonalize it, and add the lowest 4 bands' density (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 ), 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 is replaced by .

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.

[12]lines 639–658

The loop itself.

[13]lines 664–723

In the loop's output the density change fell by a factor of 0.42 per iteration, from 0.46 to in 21 iterations.

The bands along a path

With the potential converged, the bands at any come from one diagonalization each, with no further iteration. The path runs from , the center of a hexagonal face of the Brillouin zone, to its center , and on to , 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 and threefold degenerate. The lowest conduction state is not at but most of the way to , 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.

[14]lines 745–761

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

[15]lines 768–814

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 , and . The next cell compares against those: the largest difference among the valence bands, and the difference for the lowest conduction state at . 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 lower.

[16]lines 830–853

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 , and and the same order of states; the cell above measured how far the energies differ.

[17]lines 865–884