KRONOS: Building a Plane-Wave DFT Engine From Scratch in C++20 and CUDA

Since March my nights-and-weekends repo hasn’t been a web app or a CLI tool. It’s been KRONOS — Kohn-Residual Optimized Numerics Over Silicon — an ab initio plane-wave Density Functional Theory engine written from scratch in C++20, with a CUDA backend behind a GPU abstraction layer. Feed it a crystal — lattice vectors, atom positions, a pseudopotential per element — and it computes the ground-state total energy, the electron density, Kohn-Sham eigenvalues, forces on every atom, and band structures. On the systems I’ve validated, it agrees with Quantum ESPRESSO — the reference production code in this field — to fractions of a milli-electron-volt per atom. There’s a project page with the short version; this post is the long one.

I should be upfront about the why, because the README doesn’t carry a mission statement and I’m not going to retrofit one. I write backend systems for a living; DFT is about as far from a Django view as software gets — quantum mechanics, spectral methods, and iterative eigensolvers stacked on top of each other. I wanted to know whether I could understand that machinery all the way down — not use it, understand it — and the only test of understanding I trust is building the thing and checking it against people who do this for a living. Everything else in the repo — the validation matrix, the 450+ tests, the physics notes in docs/ — is downstream of taking that test seriously.

What a DFT engine actually computes#

Strip the physics vocabulary and the problem statement is simple. Input: the geometry of a crystal. Output: the energy of its electrons in their lowest-energy configuration, and the force that configuration exerts on each nucleus. Those two numbers are the currency of computational materials science — compare energies between two candidate structures and you know which one nature prefers; follow the forces downhill and you find the geometry a material actually relaxes into.

The reason this needs an engine rather than a formula is that electrons interact with each other. The exact quantum state of N electrons is a function of 3N coordinates, and solving for it directly is exponentially hard. DFT is the reformulation that made the problem computable: the ground-state energy turns out to be fully determined by the electron density, an ordinary 3D scalar field, and the Kohn-Sham construction swaps the interacting electrons for fictitious independent ones moving in an effective potential built from that density. The catch — and the thing that shapes the entire architecture — is circularity. The potential depends on the density; the density depends on the orbitals you get by solving the eigenproblem in that potential. So a DFT engine is, at its core, a fixed-point solver: guess a density, build the potential, solve the eigenproblem, compute the new density, repeat until input and output agree. That’s the self-consistent field (SCF) loop, and everything else is machinery serving it.

One piece is genuinely approximate: the “exchange-correlation” functional, which packages everything the independent-electron fiction gets wrong into one term. KRONOS ships the classic local-density approximations (Perdew-Zunger LDA, spin-polarized LSDA) built in, and PBE/PBEsol gradient-corrected functionals via libxc. Choosing a functional is choosing your error bar — the engine’s job is to make sure the numerics contribute nothing on top.

Why plane waves fit crystals#

Every electronic-structure code has to pick a basis: the set of functions you expand orbitals in so the eigenproblem becomes linear algebra. For molecules people use Gaussians centered on atoms. For crystals, the right answer falls out of the physics. A crystal is periodic, and Bloch’s theorem says its eigenfunctions are periodic functions times a phase factor — and the natural expansion of a periodic function is a Fourier series. Plane waves are the Fourier basis. One vector of complex coefficients per orbital, indexed by reciprocal-lattice vectors G, truncated by a single knob: keep every G with |k+G|² ≤ ecutwfc.

That single knob is the engineering win. Basis-set quality is one monotonic parameter — raise the cutoff, converge the answer, done. No basis-set superposition errors, no judgment calls about which functions to center where. The second win is structural: the kinetic-energy operator is diagonal in reciprocal space (it’s just |k+G|² per coefficient), while the local potential is diagonal in real space (pointwise multiply). An FFT converts between the two in O(N log N). So applying the Hamiltonian — the hot operation the eigensolver calls thousands of times — never materializes a matrix. It’s: multiply by |k+G|² in G-space, FFT to real space, multiply by V_eff(r), FFT back, add the nonlocal pseudopotential term (a couple of GEMMs against projector vectors). The entire engine is architected around making that FFT sandwich fast and exact.

The price of plane waves is that they resolve everything, everywhere, at uniform resolution — and near a nucleus, core electrons oscillate so sharply that representing them would need absurd cutoffs. The standard fix, which KRONOS implements, is pseudopotentials: replace the nucleus-plus-core-electrons with a smooth effective potential, and simulate only the valence electrons that do the chemistry. KRONOS reads the norm-conserving UPF v2 files that the Quantum ESPRESSO ecosystem publishes, and it validates them at parse time — the integrated atomic density has to match the declared valence charge to within 1e-3, and a failed check aborts with the file path and a suggested replacement rather than letting a corrupt pseudopotential produce quietly wrong energies downstream.

One fixed-point loop wrapped around FFTs#

Here’s the actual pipeline, as it runs. The outer flow is linear — parse, build, iterate, post-process — and the inner loop is the fixed point:

flowchart TD
    A[Parse YAML input - strict schema] --> B[Load and validate UPF pseudopotentials]
    B --> C[Generate Monkhorst-Pack k-points - spglib symmetry reduction]
    C --> D[Build plane-wave basis and FFT grid]
    D --> E[Initialize density from atomic superposition]
    E --> SCF
    subgraph SCF [SCF loop]
        L1[FFT density to G-space] --> L2[Hartree potential - Poisson solve, diagonal in G]
        L2 --> L3[Exchange-correlation potential in real space]
        L3 --> L4[Assemble V_eff]
        L4 --> L5[Davidson eigensolver per k-point]
        L5 --> L6[Fermi level by bisection - occupations]
        L6 --> L7[New density from occupied orbitals]
        L7 --> L8{Energy AND density converged?}
        L8 -- no --> L9[Pulay/DIIS density mixing] --> L1
    end
    L8 -- yes --> P[Ewald ion-ion energy - Hellmann-Feynman forces]
    P --> O[JSON summary + HDF5 output]

A few of the boxes deserve a sentence each. The Hartree potential — the classical electron-electron repulsion — is a Poisson equation, which in G-space collapses to a pointwise division: one line of code instead of a PDE solver. The eigensolver is a blocked Davidson iteration: it never diagonalizes the full Hamiltonian (catastrophic at ~10⁴–10⁵ plane waves per k-point), it iteratively refines a subspace three times the number of bands you asked for, preconditioned by the kinetic-energy diagonal. When Davidson diverges — bad pseudopotentials, pathological inputs — KRONOS detects the blown residual and auto-switches that k-point to LOBPCG rather than dying. And the new density is never fed straight back into the next iteration: raw fixed-point iteration oscillates for most real systems, so a Pulay/DIIS mixer extrapolates from the history of previous densities, with Kerker damping for metals.

Convergence is a dual criterion — the energy change and the G-space density residual both have to fall below threshold — and when the loop can’t get there, KRONOS notices the classic failure mode (energy seesawing between iterations) and aborts with a diagnostic instead of burning CPU forever. After convergence: Ewald summation for the ion-ion energy, Hellmann-Feynman forces symmetrized by spglib, and optionally a BFGS loop on top that moves the atoms downhill and re-runs SCF until the structure relaxes. Results land as a JSON summary plus HDF5 binaries — written atomically via temp-file-and-rename, because a checkpoint corrupted by a killed job is worse than no checkpoint — with structured JSON logs (timings, GPU memory) on stderr.

Running it looks like this:

cmake -B build -S .
cmake --build build -j$(nproc)
cd build && ctest --output-on-failure   # 450+ GoogleTest cases

./build/src/kronos examples/si_bulk.yaml

The input is a contract, not a suggestion#

Inputs are YAML, and the schema is strict — unknown keys are rejected, hard:

system:
  lattice: [[5.43, 0, 0], [0, 5.43, 0], [0, 0, 5.43]]
  atoms:
    - {symbol: Si, position: [0.0, 0.0, 0.0]}
    - {symbol: Si, position: [0.25, 0.25, 0.25]}

calculation:
  type: scf          # scf | relax | bands | dos
  ecutwfc: 60.0      # wavefunction cutoff (Ry)
  ecutrho: 240.0     # density cutoff (Ry), default = 4 * ecutwfc
  kpoints: [8, 8, 8, 0, 0, 0]  # Monkhorst-Pack grid + shift
  xc: PBE            # LDA_PZ | LDA_PW | PBE | PBEsol
  smearing: marzari-vanderbilt

pseudopotentials:
  Si: pseudopotentials/Si_ONCV_PBE-1.2.upf

convergence:
  energy: 1e-8       # Ry
  density: 1e-9
  max_scf_steps: 100

The conventional move in scientific software is the opposite: Fortran-namelist-style inputs where a misspelled key is silently ignored and the code proceeds on defaults. I rejected that deliberately, because the worst failure mode in this domain is not a crash — it’s a plausible number. Typo ecutwfc into something the parser doesn’t recognize and a lenient code happily runs at its default cutoff, hands you an energy that looks fine to four digits, and poisons every comparison you build on it. KRONOS aborts on the unknown key instead. A web engineer would call this schema validation at the API boundary and consider it table stakes; in this field it still counts as an opinion.

Double precision is not a style choice#

The document I’d point any engineer at first is the numerical-precision notes, because it explains the constraint that shapes everything: a DFT total energy is a difference of large numbers. Kinetic, Hartree, exchange-correlation, and ion-ion terms are each tens of Rydberg; the chemistry you actually care about lives at the milli-eV-per-atom level, which means every term needs relative accuracy better than about 1e-7 just to leave the signal standing. The error budget stacks up fast: FFT round-trip noise, eigensolver residuals at 1e-8 Ry, SCF density tolerance at 1e-6 — summed, that’s roughly the 1e-6 Ry/atom you can afford. A 64³ FFT in fp32 alone blows that budget by about five orders of magnitude. So KRONOS hard-codes std::complex<double> for wavefunction coefficients through every physics path, and the unit tests assert the fp64 FFT round-trip error stays at machine-precision levels on 16³ through 64³ grids.

Determinism gets the same treatment as precision. Davidson starts from random vectors — seeded with a fixed 42, so the same input reproduces the same trajectory and the same final energy to all digits on the same hardware. On GPU builds, CUBLAS_WORKSPACE_CONFIG=:4096:8 pins cuBLAS reduction order, because ±1 ULP run-to-run wobble is harmless in production and fatal to bit-for-bit regression tests.

My favorite counterintuitive detail lives in the Fermi-level solver. The Fermi level is found by bisection, and the tolerance is 1e-10 Ry — deliberately not tighter. Crank it to 1e-13 and, on toy pseudopotentials with no smearing, the Fermi level starts shifting by floating-point noise between SCF iterations; occupations seesaw, the density never settles, and the loop oscillates forever. Tighter tolerance, worse convergence. It’s the kind of empirically-tuned constant that looks arbitrary in a diff and is load-bearing in practice — which is exactly why it’s documented instead of left as folklore.

One namespace is allowed to know about CUDA#

The GPU design decision I care most about is a boundary, not a kernel. The kronos::gpu namespace is the single point of coupling between physics code and vendor APIs. Physics code in hamiltonian/, basis/, and solver/ calls three thin interfaces — a GPUFFTGrid with forward/inverse, a gemm/zdotc pair, and gpu_malloc/gpu_free/gpu_memcpy_* — and never sees cuda_runtime.h. Behind that boundary sit four backends: CUDA (cuFFT/cuBLAS), HIP (rocFFT/rocBLAS), Metal, and a CPU stub build where every GPU call throws GPUNotAvailableError so a CPU-only binary can’t silently wander down a device path. Adding a backend means touching src/gpu/ and nothing else.

Device memory ownership is explicit and RAII all the way down. DeviceBuffer<T> owns exactly one GPU allocation, is move-only, and frees in its destructor — there is no code path where a raw cudaMalloc pointer escapes into physics code to leak or double-free:

template <typename T>
class DeviceBuffer {
public:
    explicit DeviceBuffer(size_t n) : size_(n) {
        if (n > 0) ptr_ = static_cast<T*>(gpu_malloc(n * sizeof(T)));
    }
    DeviceBuffer(DeviceBuffer&& other) noexcept
        : ptr_(other.ptr_), size_(other.size_) {
        other.ptr_ = nullptr;
        other.size_ = 0;
    }
    DeviceBuffer(const DeviceBuffer&) = delete;            // no implicit
    DeviceBuffer& operator=(const DeviceBuffer&) = delete; // device copies
    ~DeviceBuffer() { free(); }

    void upload(const std::vector<T>& host_data);
    std::vector<T> download() const;
    T* data() { return ptr_; }
    // ...
};

The custom CUDA kernels themselves are deliberately boring, and that’s the point. The flops live in cuFFT and cuBLAS; my kernels exist to keep data resident on the device between library calls, so wavefunctions don’t ping-pong across PCIe inside the Hamiltonian’s FFT sandwich. Scatter plane-wave coefficients onto the FFT grid, gather them back, pointwise-multiply by the potential, apply the diagonal kinetic term:

// Apply kinetic energy: hpsi[i] = |k+G|^2 * psi[i]  (diagonal in G-space)
__global__ void apply_kinetic_kernel(
    cuDoubleComplex* hpsi,
    const cuDoubleComplex* psi,
    const double* kinetic_energies,
    int npw)
{
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < npw) {
        double ke = kinetic_energies[i];
        hpsi[i] = make_cuDoubleComplex(
            ke * cuCreal(psi[i]),
            ke * cuCimag(psi[i]));
    }
}

Every element is independent, every access explicit — the explicit-over-clever version of what a fused kernel might do faster. If profiling ever says this glue matters, it can get clever then. Honesty requires the status label too: the CUDA and HIP paths are fp64 and functional, but experimental — they have not yet completed the external validation matrix against Quantum ESPRESSO that the CPU path has. That’s the very next roadmap line, not a footnote.

Shipping a GPU path I refuse to trust#

The Metal backend is my clearest “conventional move versus what I did” call in the project. Apple GPUs have no fp64 ALUs — Metal Shading Language will not compile double, on any toolchain, with no emulation path. Given everything the precision section just said, the conventional move is obvious: skip Apple hardware entirely, since an fp32 DFT engine cannot produce validation-grade energies, full stop.

I built the Metal backend anyway — as an explicitly fenced fp32 research tier — because the machine on my desk is Apple Silicon, and most GPU bugs in this codebase are not precision bugs. They’re orchestration bugs: a scatter map built for the wrong grid, a buffer uploaded before it’s sized, an FFT normalization applied twice. An fp32 device path that runs the whole SCF loop locally catches every one of those at laptop iteration speed, without waiting for time on an NVIDIA box. The fencing is the actual design work: the Apple path only activates via hardware.apple_fast_mode: true in the input (or --apple-fast-mode on the CLI), activation emits a logger warning, and the validation suite refuses to run under the flag — a Metal build without it falls back to CPU fp64. The narrowing from complex128 to complex64 happens at exactly one place, the device boundary, and Apple’s unified memory is exposed through the same gpu_malloc interface via shared-mode MTLBuffers, so physics code can’t tell the backends apart. Useful and untrusted, and the type system plus the input schema make sure those two properties never blur.

Numbers I can defend#

Fabricating confidence is easy in scientific code; the repo’s answer is a reproducible validation matrix, with inputs, tolerances, and scope limits committed in VALIDATION.md. The CPU fp64 matrix is 12 paired calculations against Quantum ESPRESSO 7.6 — silicon in several configurations, aluminum, copper, water, MgO, graphene, and spin-polarized BCC iron. All 12 pass the 2 meV/atom energy threshold the DFT reproducibility community uses; the required nonzero forces agree within 3.23e-5 Ry/bohr, stress within 1.12e-6 Ry/bohr³, and constrained-magnetic iron within 0.5455 meV/atom. On the gamma-only silicon LDA benchmark, the disagreement with QE is 0.07 meV/atom. Analytic Hellmann-Feynman forces are separately checked against finite-difference energy derivatives — the test that catches a wrong force expression even when the energy looks right — alongside physics-invariant tests (Hamiltonian hermiticity, sum rules, Parseval) and frozen regression baselines, 450+ GoogleTest cases in all, run in CI.

# Reproduce the QE validation calculation yourself
./build/src/kronos examples/si_qe_validation.yaml

I want to be precise about what those numbers do and don’t claim. They say: within the validated scope — norm-conserving pseudopotentials, LDA/GGA, these systems, CPU fp64 — KRONOS’s numerics match the reference implementation to well under the community threshold. They do not say KRONOS is fast, complete, or a replacement for anything. It’s a v0.1 research engine, and the validation file states its own scope limits on purpose.

Where it goes from here#

The roadmap is in the README and I’ll only repeat what’s written there: v0.2 takes the CUDA/HIP backends through the external QE validation matrix and adds device CI; v0.5 brings finite-displacement phonons and the foundations for noncollinear spin; v0.8 adds PAW pseudopotential support; v1.0 is the production cut with a Python package (pybind11 bindings already exist behind a CMake flag) and full CI/CD; v2.0 targets hybrid functionals and non-collinear magnetism.

The code is GPL-3.0 at github.com/kdpisda/kronos, the docs — including physics notes that build Kohn-Sham DFT up from scratch for engineers — live at kdpisda.github.io/kronos, and the short summary is on the project page. My payments post ends with the observation that the browser is hostile territory and the only thing you believe is what you verify yourself. Floating point is the same kind of territory. Different stack, same discipline: don’t trust the number until you’ve checked it against ground truth you didn’t produce.