§1
Why do this
Scientific software usually earns trust in three ways: comparison with analytic solutions, with experiment, and with other programs. All three work, but they only cover the cases that were tested, and when a result deviates it is hard to say which layer the deviation came from.
The Sleipner A platform sinks
The concrete base of the Sleipner A platform in the North Sea failed during a ballast test and sank. The investigation found that poorly shaped elements in the linear-elastic finite-element analysis underestimated the shear stress in the tricell walls by about 47%, so the reinforcement was insufficient. The program ran normally and reported no error.
Lesson: mesh quality and discretisation error should be checkable preconditions, not a matter of engineering judgement.
The Δ-gauge and reproducibility
Lejaeghere et al. compared the PBE equations of state of 71 elemental crystals across 15 solid-state codes and 40 pseudopotential or basis-set schemes. Recent methods agree well; older schemes deviate noticeably, and the source of each difference could only be located by item-by-item comparison.
Lesson: results depend on many implicit numerical settings. Those settings belong in the specification, tied to an error estimate.
Parallelism makes things harder. An MD code decomposes space over thousands of MPI ranks; if an atom is lost or duplicated as it crosses a subdomain boundary, the energy curve may show only a barely visible jump. The order of floating-point reductions changes with the number of ranks, so the same input gives bitwise-different trajectories on machines of different size. These errors have nothing to do with physics. They are errors in concurrent protocols, and concurrent protocols are exactly what TLA+ is good at.
§2
Core claim: a result is only as trustworthy as the weakest link
We split the path from equation to floating-point number into six layers and write an explicit proof obligation between each pair. The first three obligations belong to PDE theory and numerical analysis, the fourth to software specification, the fifth to floating-point analysis.
- L1
Physical model
The equations themselves: Kohn–Sham, Newton or Langevin equations of motion, linear elasticity.
- Maths Well-posedness: a solution exists, is unique and depends continuously on the data
- L2
Continuous problem
A variational form in a function space, e.g. find u ∈ H¹ with a(u, v) = ℓ(v) for all v.
- Maths Discretisation error estimate: uh → u, with a rate
- L3
Discrete problem
A finite-dimensional algebraic system: stiffness matrices, plane-wave coefficients, particle coordinates.
- Maths Algorithmic convergence: when the iteration stops at step k, the algebraic error is controlled
- L4
Abstract algorithm
SCF iteration, velocity Verlet, conjugate gradients — described mathematically, in exact arithmetic.
- TLA+ Refinement: every step of the parallel program corresponds to one step of the abstract algorithm, or leaves the abstract state unchanged
- L5
Parallel program
MPI ranks, messages, shared-memory threads, GPU streams, checkpoint and restart.
- Floating point Rounding-error bound: the distance between the floating-point result and the exact-arithmetic result
- L6
Floating-point execution
The number actually computed on IEEE 754 hardware.
- Maths Discretisation and algebraic error: from PDE theory and the theory of the algorithm
- Floating point Rounding error: from floating-point analysis
- TLA+ Precondition: the program really computes ũh(k) — guaranteed by refinement
This inequality is the backbone of the project. Each term must be backed by a theorem, and the premise that makes it hold — that the parallel program faithfully executes the abstract algorithm — must be guaranteed by the software-layer specification. The gap between the model and real physics is not in this inequality: that is validation, not verification (§7).
§3
Two layers, one bridge, one foundation
PDE theory and numerical analysis
- Well-posedness: Lax–Milgram, inf-sup conditions, existence of minimisers of the Kohn–Sham energy
- Discrete convergence: a priori estimates (Céa’s lemma, interpolation, plane-wave truncation) and a posteriori estimates
- Structure preservation: symplecticity, time reversibility, conservation of momentum and charge, crystal symmetry
- Iterative convergence: conditions and stopping criteria for SCF fixed points, Newton and Krylov methods
Tool: Lean 4 with Mathlib, which already has Hilbert spaces, measure theory and Lax–Milgram. Strategy: state theorems first as axioms with a literature source so that upper layers can depend on them at once, then replace the axioms with machine proofs. Every unproved axiom is registered in the trust ledger.
TLA+ specification and model checking
- Safety: atoms are never lost or duplicated, every degree of freedom has exactly one owner, parallel assembly has no write conflicts, no deadlock
- Liveness: messages are eventually delivered; the SCF loop either converges or reports failure in finitely many steps
- Protocols: domain decomposition and halo exchange, distributed FFT transposes, checkpoint and restart, adaptive refinement and load rebalancing
Tools: PlusCal for algorithms; TLC for exhaustive checking of small instances; Apalache for symbolic checking; TLAPS for proofs at any scale. Link to the code: the implementation emits an event log at run time, and trace validation checks that the log is a behaviour the specification allows.
Refinement carries the theorems to the parallel program
The abstract algorithm is a TLA+ specification A whose state is the global density or phase-space point, and whose steps are one mathematical iteration each. The parallel implementation is a specification P whose state is spread over processes. A refinement mapping projects the state of P onto that of A — for example by stitching local densities into the global one. Once P ⇒ A is proved, the convergence theorems about A hold for the parallel program.
Performance work is handled the same way: every optimisation is a new refinement, and must be proved to still implement the level above.
Floating point is a first-class citizen
IEEE 754 addition is not associative. Floating-point results never enter a conclusion directly: the output of a large solver goes to a proved checker that uses either exact rational arithmetic or outward-rounded interval arithmetic. The only floating-point fact we then have to trust is “round-to-nearest returns the nearest double”.
Global reductions use a deterministic summation order, so results do not depend on the number of processes or threads. The rounding model is part of the specification and sits in the same error budget as the discretisation and algebraic errors.
§4
The proof-obligation matrix
What each of the three domains must prove at each layer. This table is the project’s work list: every cell must eventually correspond to a specification, a theorem or a registered assumption.
| Layer | DFT · electronic structure | MD · atomic trajectories | FEM · continuum mechanics |
|---|---|---|---|
| Model | Kohn–Sham equations, a nonlinear eigenvalue problem H[ρ]ψi = εiψi, ρ = Σfi|ψi|² | Hamiltonian dynamics; for thermostatted runs, Langevin dynamics and its Fokker–Planck equation | Elliptic, parabolic and hyperbolic PDEs, e.g. linear elasticity −∇·σ(u) = f and its weak form |
| Continuum | Lieb variational principle; existence of minimisers for LDA-type Kohn–Sham models; self-adjointness of the Hamiltonian | Existence and uniqueness for Lipschitz forces; conservation of energy, momentum and angular momentum; Liouville; invariant measure and ergodicity of Langevin dynamics | Well-posedness via Lax–Milgram; inf-sup conditions for mixed problems (Stokes, nearly incompressible elasticity); regularity sets the attainable rate |
| Discrete | A priori error of the plane-wave cutoff; Brillouin-zone k-point integration error; pseudopotential model error registered separately | Symplecticity and time reversibility of velocity Verlet; backward error analysis: a modified Hamiltonian keeps the energy error O(Δt²) for exponentially long times | Céa’s lemma plus interpolation: ‖u − uh‖H¹ ≤ C hk|u|Hk+1; mesh shape regularity as an explicit premise; reliability and efficiency of a posteriori estimators |
| Algorithm | Local convergence of SCF with Anderson/Pulay mixing; orthogonality in LOBPCG and Davidson; the SCF state machine: converge, oscillate, restart | Neighbour lists: no pair is missed while the maximum displacement since the last rebuild is below half the skin; convergence of SHAKE/RATTLE | Convergence and stopping criteria of CG/GMRES with preconditioners; convergence of the SOLVE → ESTIMATE → MARK → REFINE loop |
| Parallel | k-point / band / plane-wave data distribution; distributed 3D FFT transposes lose and duplicate nothing; every process sees the same ρ after reduction | Atom migration is conserved under domain decomposition; ghost exchange is consistent; no deadlock; state is equivalent across checkpoint and restart | Parallel assembly is write-conflict free; every degree of freedom has exactly one owner; the mesh stays conforming after distributed refinement; rebalancing loses no cells |
| Floating point | Accumulated rounding in orthogonalisation; a bound on the numerical deviation of the charge ∫ρ = N | Rounding accumulation over long integrations; bitwise irreproducibility from reduction order | The limit that a condition number κ ∼ h−2 puts on attainable accuracy; quadrature error |
§5
Trust tiers
Not everything can be proved at once. Every component states honestly which tier it has reached and lists the assumptions it depends on in the trust ledger. Tiers only go up, and are never overstated.
Tested
Regression tests, analytic solutions, comparison with established codes. Holds only for the inputs tested.
Paper proof
A proof on paper with a literature source and explicit assumptions, not yet machine-checked.
Model-checked
TLC or Apalache has exhausted every behaviour of a finite instance.
Machine proof
A Lean, Rocq or TLAPS proof valid at any scale, for the abstract model, exact arithmetic or a line-by-line model of the code.
End to end
Machine proof all the way down to executable code and floating-point rounding. Precedent: Boldo et al.’s proof of a C program for the 1D wave equation.
§6
Principles
Specification before code
Every module has a mathematical statement and a TLA+ specification before it has an implementation.
Assumptions are visible
Regularity, mesh conditions, step-size limits and convergence thresholds are written into the specification, not hidden in input-file defaults.
Honest tiers
Better to mark T0 than to overstate T3. The trust ledger is the most important file in the project.
A small trusted kernel
Only a few kernels need T4; upper layers inherit their guarantees through composition theorems.
Floating point enters the specification
The rounding model and its error bounds join the discretisation error in a single error budget.
Reproducible by default
The same input gives bitwise-identical results at any process count, unless the user explicitly gives this up.
Correct first, fast later
First obtain a provable version, then optimise without breaking the proof. Every optimisation is a new refinement.
Verify the checker, not the solver
Large numerical libraries such as LAPACK or FFT libraries are not proved; their outputs are checked by proved checkers. Solvers are replaceable; the guarantee comes from the checker.
Every run comes with a report
Each line of the report cites a claim in the trust ledger, from which its tier and assumptions can be looked up; the report ends with the assumptions the run used.
Every output has a theorem
Every number the program prints must correspond to a machine-proved theorem. Where the result is an approximation — iteration, discretisation, k-point sampling, finite temperature, floating point — what is proved is an error bound. Tests, model checking and sanitizers are supporting evidence only. An output without such a theorem is unfinished and is listed as a gap.
Open
Specifications, proofs, benchmarks and the trust ledger are to be public, so that anyone can re-check them.
§7
What we do not do
Judge the physics
The gap between exchange-correlation functionals, force fields or constitutive laws and the real world is validation. This project does verification: solving the given equations correctly.
Replace existing software at the outset
Quantum ESPRESSO, VASP, LAMMPS, GROMACS, deal.II and FEniCS are benchmarks and references, not targets for replacement.
Prove all the code
Input parsing, I/O and visualisation can stay at T0, as long as they do not enter the error budget.
§8
Roadmap
Current goal G1 — feature parity with DFTK set 29 Sep 2026
Bring mini-DFT to the main feature set of DFTK.jl. A feature is done only when it is implemented, agrees with DFTK on the same discrete problem, and its output carries a Lean certificate. Order: complete everything in 1D first, then generalise the Lean model to d dimensions for 3D, and finally infrastructure and external interfaces. FEM moved to the parallel goal G2; for MD (phase 2), a basic version, mini-MD v0, was built on 1 October 2026. Status →
Parallel goal G2 — certified finite elements against scikit-fem set 30 Sep 2026
Cover the main features of scikit-fem with our own C implementation, with scikit-fem as the independent reference on the same discrete problem. Same completion standard as G1, and wherever possible the certificate speaks about the true solution of the continuous problem, not just the discrete one. 1D first, then 2D triangular meshes. Status →
Parallel goal G3 — certified CALPHAD against pycalphad set 1 Oct 2026
Read TDB thermodynamic databases and compute phase equilibria with our own C implementation, with pycalphad as the independent reference on the same database and the same phases. Same completion standard as G1 and G2, and the certificate speaks about the true global equilibrium of the database model, not just the local solution the driver found. Binary systems first (point equilibria, phase diagrams, invariant reactions, magnetism), then multicomponent and sublattice models. Status →
- PHASE 0Minimal DFTIn progress
A one-dimensional periodic reduced Hartree model, plane waves, serial. Walk through L1–L6 once, and build the infrastructure the whole project reuses: the trust ledger, certificate checkers, the run-report format, trace validation.
Exitevery obligation reaches its target tier; four benchmark groups pass; every run prints a trust report.
- PHASE 1DFT extensionsPartly ahead of plan
1D LDA exchange-correlation (convexity is lost; uniqueness becomes a local result); k-point sampling and its MPI parallelisation — the project’s first real parallel protocol; 3D plane waves with local pseudopotentials, cross-checked against DFTK.jl.
Exitk-point parallel runs are bitwise identical on 1 to 64 processes and every run log passes trace validation.
- PHASE 2MD and FEMExit criteria met
Reuse the phase-0 framework for two minimal versions: a 1D Lennard-Jones chain with velocity Verlet, and 1D Poisson with P1 elements. TLA+ modules for halo exchange, atom migration, parallel assembly and checkpointing.
Exitboth minimal versions print trust reports; Céa’s lemma and Verlet symplecticity reach T3.
- PHASE 3Trusted kernelsNot started
Bring a few kernels to T4: the plane-wave kinetic operator, FEM element-stiffness assembly, one Verlet step — with machine-checked floating-point error bounds connected to the code.
Exitthe floating-point error bounds of these kernels are machine-checked and enter the global error budget.
- PHASE 4BenchmarkingNot started
Standard benchmarks against established codes: the Δ-gauge for DFT, NVE energy drift and radial distribution functions for MD, the method of manufactured solutions for FEM convergence rates.
Exitpublished benchmark results, each with its trust ledger.
Detailed, item-by-item status is on the progress page and on one page per domain: DFT, FEM, MD, CALPHAD.
§9
First milestones: mini-DFT and mini-FEM
Before extending to MD and FEM, we walk through L1 to L6 with the smallest DFT program we could design. It is small enough that every layer’s obligation can be written down, yet it keeps the full skeleton of a Kohn–Sham calculation: a nonlinear eigenvalue problem, self-consistent iteration, plane-wave discretisation.
| Model | One-dimensional periodic reduced Hartree–Fock (rHF), no exchange-correlation. Periodic interval [0, L), N electrons, spin-degenerate. Energy: kinetic + external + ½ Hartree, with the Hartree kernel 4π/q² on nonzero frequencies. |
| Potential | Smooth periodic: constant, cosine, periodised Gaussian double well (optionally repeated). The 1D Coulomb potential is not integrable at the origin and is not used. |
| Discretisation | Plane waves with |k| ≤ K; Γ only, or Nk equally weighted k-points (Γ-centred or Monkhorst–Pack). A real-space grid of at least 4K + 1 points guarantees that density and Hamiltonian action are alias-free. |
| Occupation | Zero temperature with a gap checked in interval arithmetic at every iteration; or Fermi–Dirac, Gaussian, Methfessel–Paxton or Marzari–Vanderbilt smearing with a common chemical potential. |
| Algorithm | SCF with linear or Anderson mixing and optional Kerker preconditioning; dense diagonalisation; shared-memory thread parallelism over k-points with a fixed summation order. |
Why rHF
Its energy is convex in the density matrix, and the Hartree term is strictly convex on nonzero frequencies, so the ground-state density is unique. Within the DFT family it has the cleanest mathematics — a good first target for machine proof. Adding LDA exchange-correlation loses convexity; that step belongs to phase 1.
LAPACK solves, the checker guarantees
LAPACK stays at T0. After convergence the solver writes a certificate; a checker written and compiled in Lean re-checks it in exact rational arithmetic and returns rigorous intervals for the energy, free energy, truncation limit, density, eigenvalues, band structure and density of states. Its soundness is a Lean theorem: it trusts no number from the solver, so a solver bug can only make the check fail or the interval wider.
mini-FEM follows the same pattern for goal G2: 1D Poisson with P1 elements (Dirichlet, Neumann and Robin boundaries) and the beam equation with cubic Hermite elements (clamped, pinned, free and sliding ends), with right-hand sides built from polynomials, exponentials, sines and cosines. The Lean checker returns intervals for the true solution of the continuous problem at every node, and a rigorous bound on the maximum error over the whole interval. Adaptive refinement uses these proved element-wise bounds as its marking criterion.
mini-MD (1 October 2026) does the same for a 1D Lennard-Jones chain with velocity Verlet: the checker bounds the distance from the true solution of Newton’s equations at every step, and from the initial state over about one vibration period. mini-CALPHAD (goal G3) certifies binary phase equilibria from TDB databases — magnetic phases and whole phase diagrams included — with the chemical potentials, driving forces and which phases can appear in any equilibrium. See the MD and CALPHAD pages.
§10
Decisions and open questions
Provers: Lean 4 for the mathematics and the floating-point model; the code-level link is still open
- Division
- Lean 4 + Mathlib for L1–L4 (well-posedness, discretisation error, structure preservation, iterative convergence) and for two parts of L6: the certificate checkers (exact rationals, compiled to executables) and a line-by-line binary64 model of the interval-arithmetic kernel.
- Frama-C/WP for the absence of run-time errors in the trusted C kernel.
- Connecting the Lean floating-point results to the actual C code — candidates are Rocq with Flocq, VST and CompCert, or Frama-C’s floating-point model — is not started.
- Why
- The mathematics layer is the bulk of the work; Mathlib’s analysis and measure theory form one unified library, and most AI-assisted proving tools support Lean first. “Verify the checker, not the solver” shrinks the floating-point obligations enough that a small binary64 model in Lean suffices. For code-level verification, Rocq’s toolchain and the T4 precedent remain strong, so the tool for that step is not yet chosen.
- Cost
- The correspondence between the Lean model and the C code rests on line-by-line review and bitwise differential tests, not on a machine proof; interval-arithmetic correctness therefore stops at T3.
Implementation language: C11, split into a proved kernel and an untrusted driver
- Division
- Trusted kernel — roots of unity, grid density and charge intervals, interval-arithmetic residual bounds and inertia counting. Restricted C with formal annotations; Frama-C/WP proves no run-time errors (T4).
- Certificate checkers in Lean, exact rationals, soundness a Lean theorem.
- Driver — SCF loop, LAPACK, threads, I/O, reports, event logs, certificate output. Stays at T0: every number it computes must pass the kernel or a checker before it appears as a guaranteed result.
- Python and Julia appear only in benchmarks and trace-validation tooling, outside the trust chain.
- Why
- The mature code-level verification routes all start from C. LAPACK and FFTW have C interfaces. The main risk, memory safety, is handled by proof in the kernel and by sanitizers in the driver. Compilation forbids fused multiply-add contraction and any fast-math optimisation, so the compiled arithmetic matches the proofs.
- Not chosen
- Rust (its verifiers do not connect to the code-level tools above), Julia or Python (no route to code-level proof), Fortran (few verification tools), Lean itself (its floating point has no formal semantics).
Open questions
- Continuum theory of 1D rHF. The literature treats the 3D Coulomb interaction; for the 1D periodic Hartree kernel we need to find a reference or supply a proof.
- How deep to formalise Sobolev spaces — what to axiomatise first, and which theorems to prove first.
- Randomness. Langevin dynamics and random initial guesses involve probabilistic properties; TLA+ only describes nondeterminism.
- GPUs. How asynchronous streams and weak memory models enter a TLA+ specification.
§11
References
- B. Jakobsen, F. Rosendahl. The Sleipner platform accident. Structural Engineering International 4(3), 1994.
- K. Lejaeghere et al. Reproducibility in density functional theory calculations of solids. Science 351, aad3000, 2016.
- S. Boldo, F. Clément, J.-C. Filliâtre, M. Mayero, G. Melquiond, P. Weis. Wave equation numerical resolution: a comprehensive mechanized proof of a C program. Journal of Automated Reasoning 50(4), 2013.
- S. Boldo, F. Clément, F. Faissole, V. Martin, M. Mayero. A Coq formal proof of the Lax–Milgram theorem. CPP 2017.
- J. P. Solovej. Proof of the ionization conjecture in a reduced Hartree–Fock model. Inventiones Mathematicae 104, 1991.
- E. Cancès, A. Deleurme, M. Lewin. A new approach to the modeling of local defects in crystals: the reduced Hartree–Fock case. Communications in Mathematical Physics 281, 2008.
- E. Cancès, R. Chakir, Y. Maday. Numerical analysis of the planewave discretization of some orbital-free and Kohn–Sham models. ESAIM: M2AN 46(2), 2012.
- M. F. Herbst, A. Levitt, E. Cancès. DFTK: A Julian approach for simulating electrons in solids. JuliaCon Proceedings, 2021.
- E. Hairer, C. Lubich, G. Wanner. Geometric Numerical Integration, 2nd ed. Springer, 2006.
- B. Leimkuhler, C. Matthews. Rational construction of stochastic numerical methods for molecular sampling. AMRX, 2013.
- P. Binev, W. Dahmen, R. DeVore. Adaptive finite element methods with convergence rates. Numerische Mathematik 97, 2004.
- L. Lamport. Specifying Systems: The TLA+ Language and Tools for Hardware and Software Engineers. Addison-Wesley, 2002.
- H. Cirstea, M. A. Kuppe, B. Loillier, S. Merz. Validating traces of distributed programs against TLA+ specifications, 2024.