A band structure from plane-wave DFT

What the button does

It starts a Web Worker so the page stays responsive, and the worker runs one of two things. With Compiled selected it fetches /dft/si.wasm — the input file, the solver and the checks, compiled to one WebAssembly module — and calls its main. The module needs no JavaScript runtime: it imports only five WASI system calls, which carry its printf output. With Interpreted it loads viv.wasm and hands it the source. Either way, for silicon the program then does this, in order:

1. Builds the crystal. Two atoms on an fcc lattice, their Goedecker–Teter–Hutter pseudopotentials, and 16 k-points chosen so that the mesh has the crystal’s full symmetry.
2. Lists the plane waves. Every reciprocal-lattice vector with kinetic energy under 10 hartree: up to 420 per k-point. They are the basis; a wavefunction is 420 complex numbers.
3. Guesses a density. A uniform one, and from it the Hartree and exchange-correlation potentials on a 24³ grid.
4. Finds the lowest 8 bands at every k-point. The Hamiltonian is never formed. It is applied: kinetic energy is diagonal in plane waves, the potential is diagonal in real space, and an FFT goes between them. LOBPCG, a block eigensolver, needs only this product of H with a wavefunction.
5. Fills the bands, builds the new density, and mixes it. Anderson mixing, over the last three densities. Then back to step 3 — 11 times for silicon, until the density stops changing at the 10⁻⁹ level.
6. Walks L → Γ → X. 61 points, in that now-fixed potential. The bands are printed as plot lines, and the figure above is drawn from them.

Step 4 is the core of the method, and it is short. This is the function, as it stands in the repo today:

def hpsi(slot, x) {
    tk = slab(TK, slot)
    idx = slab(IDX, slot)
    mask = slab(MASK, slot)
    out = czeros(len(x), NPWMAX)
    for b to len(x) {
        xb = slab(x, b)
        # local potential: to the grid, multiply in real space, back
        fill(BOX, 0.0)
        put(BOX, idx, xb)
        ifft3(BOX)
        w = BOX * VR
        fft3(w)
        setslab(out, b, take(w, idx) * mask + tk * xb)
    }
    # nonlocal: sum_IJ |beta_I> D_IJ <beta_J|x>, for every row at once
    beta = slab(BETA, slot)
    pr = x @ transpose(conj(beta))
    return out + (pr @ DMAT) @ beta
}

put scatters the coefficients into the FFT box, ifft3 takes them to real space, BOX * VR multiplies a complex array by a real one, and take gathers the result back. The last two lines are the nonlocal pseudopotential, for every band at once, as two matrix products.

How it is checked

A DFT code can be wrong and still draw a plausible band structure, so there are three layers of checks. The test list under the plot runs the last two.

The formulas. The pseudopotential’s closed-form Fourier transforms are tested against brute-force radial quadrature of their real-space definitions, and the exchange-correlation functional against libxc. A formula copied from a paper can be wrong and still run; this way the same mistake has to be made twice, independently.

The implementation. The repo has a second, independent implementation in numpy, which diagonalizes H densely with LAPACK — a different algorithm, not only a different language. The Viv code has to land on its total energy to 10⁻⁷ hartree and its band energies to 10⁻⁵ eV. The expected numbers are stored in a check file generated by running the numpy code, so an input file and its expected values always match.

The physics. Published LDA results, with tolerances loose enough to cover the spread between papers. Silicon’s gap comes out at 0.47 eV against a measured 1.17: that is LDA’s well-known failure, and reproducing it is a sharper test than getting a right answer would be.

Across the six materials that is 60 tests, run with the program compiled to C and run natively. The five lighter ones also run through both browser engines when this site is built.

Copper

Copper is the hard one. Its 3d electrons are in the valence, the pseudopotential is sharp, and it takes 1,363 plane waves at each of 128 k-points to describe them — about 26 times silicon’s problem. Compiled to C and run natively it takes 180 seconds on one core. Compiled to WebAssembly it takes about 4 minutes and 370 MB, and passes the same 13 tests — so it is offered above, with that warning. What the page shows before you run it is the native run’s recorded output.

Band energies in eV from the Fermi level, against a published LDA calculation and experiment:

The published column is Marini, Onida and Del Sole, Phys. Rev. Lett. 88, 016403 (2002), Table I; the measured one is the same table’s, from Courths and Hüfner’s review. The d bands agree with the published calculation to 0.05 eV, and sit 0.5 to 0.7 eV above the measured ones, which is what LDA does to copper.

One number is printed and not asserted. The gap above L₂′ comes out at 4.87 eV here against 5.40 in that paper, whose pseudopotential keeps the 3s and 3p shells in the valence where this one does not. Every occupied level agrees to 0.2 eV or better; why that empty one differs by half an electron-volt has not been explained, so the check file claims no tolerance for it.

Bugs fixed during development

Each of these produced a reasonable-looking band structure. The checks above are how they were found.

k-mesh symmetry

The top of silicon’s valence band is three states with one energy. It came out as two and one, 56 meV apart. The standard shifted Monkhorst–Pack mesh, written in the primitive fcc reciprocal vectors, is not closed under the cubic point group — some rotations send the shift (½ ½ ½) to (½ 0 0) — and with no explicit symmetrization the density is only as symmetric as the mesh. The fix is the four-shift mesh, which is that shift’s whole orbit: the split is zero, and 16 k-points land within 6 millihartree of 128. Every semiconductor here now tests the triplet.

LOBPCG search direction

LOBPCG reached a residual of 10⁻⁶ and then grew by a factor of five per step. Its search direction was being formed as the difference of two nearly equal vectors, which leaves no relative accuracy; once normalized, the direction and its image under H no longer matched, and the small eigenproblem at the center of the method was solving for a matrix that was not H. Building the direction directly from the expansion coefficients — Knyazev’s original form — has no subtraction to cancel. The compiled run had been passing by luck; the interpreter’s contract on eig_herm, which refuses a matrix that is not Hermitian, is what caught it.

Smearing width in copper

Metals need fractional occupations near the Fermi level. At a width of 0.02 hartree the Fermi–Dirac tail reached the top of copper’s d bands, 1.5 eV down, and emptying them pushed the Fermi level up by a quarter of an electron-volt — every band energy in the table above was off by that much. At 0.01 the shift is 0.04 eV.

Language features this needed

Viv already had arrays of floats and a symmetric eigensolver. This program needed complex arrays, fft and fft3, a Hermitian eigensolver, whole-array complex arithmetic, and slab, take and put to move rows and scattered elements without a loop. Each went into both the interpreter and the C backend at once, with a test that the two agree to the last bit, the FFT included for every length.

It also turned up two things about the compiler. Generated C never freed an array, which was fine for silicon but ran out of memory on copper; functions now release what they allocated when they return. And clang was fusing a * b + c into one instruction that rounds once where the interpreter rounds twice, a last-bit difference that an iterative solver amplifies into a different iteration count. Generated C now turns that off.

Run it yourself

git clone https://github.com/AdamZettel/viv && cd viv
cargo build --release -p viv

dft/run.sh si            # silicon, dense engine: H built and diagonalized
dft/run.sh si it         # silicon, iterative engine: the one this page runs
dft/run.sh cu it         # copper, about three minutes
dft/check.sh             # five materials, every test
dft/build-wasm.sh        # Viv -> C -> wasm, what this page runs (needs emscripten)

dft/run.sh compiles the input file, the solver and the check file to one C program and runs it; dft/build-wasm.sh hands that same C to emscripten. The same files go through the interpreter with viv run. The slowest of the five materials that run in a browser (all but copper), gallium arsenide, takes 11 seconds compiled in a browser and 34.4 interpreted. An input file is ordinary Viv that sets the globals the solver reads, and it is generated — pseudopotential parameters are read from CP2K’s table by a script, never typed. The table is CP2K’s GTH_POTENTIALS as PySCF ships it (Goedecker, Teter and Hutter 1996; Hartwigsen, Goedecker and Hutter 1998); only the numbers are used, no CP2K or PySCF code. Every file can be read here, alongside what viv doc extracts from its annotations.

What it does not do

• Total energies you can compare between structures. The ion–ion (Ewald) term is left out. It is a constant for fixed atoms and does not move a band.
• Gradient corrections. LDA only; PBE needs the density’s gradient.
• Spin, or spin–orbit coupling. Germanium’s and gallium arsenide’s valence-band tops are split by it in reality and not here.
• Symmetry reduction. Every k-point of the full mesh is computed, halved only by time reversal. That is why the density is symmetric without being symmetrized, and also why copper needs 128 of them.
• Forces, or anything but fcc lattices.