Stage 1: a small verified mathematical core¶
Stage 1 adds 30 Lean theorems for unit conversions, density-fitted Coulomb contractions, orthogonal nearest images and pair-energy accounting. The combined RHF/UHF/DF/MP2/core audit covers 58 theorems. All pass with the pinned Lean 4.35.0 and mathlib v4.35.0 environment. This implements the first milestone of the formal-verification design.
The core is optional development infrastructure. Production numerical code is unchanged. The Lean specifications are mathematical models over real numbers; the Python reference adapters and native implementation checks are separate, tested bridges. They are not extracted verified executables, and passing those checks does not establish correctness of the entire calculation engine.
Units: dimension, scale and frame¶
Units.lean represents dimensions by integer
length and energy exponents. Length is (1,0), energy (0,1) and force (-1,1).
A unit carries a strictly positive scale to a reference unit. A quantity’s type
also carries its Cartesian or fractional frame. Conversion preserves both the
dimension and frame; it cannot perform a lattice transformation. Two Lean
#check_failure checks reject dimension-changing and frame-changing terms.
This is a bounded type vocabulary, not a full SI or reciprocal-space ontology.
If unit a has scale s_a and unit b has scale s_b to the same reference, then
The intermediate scale cancels in a composed conversion. Positivity supplies
the nonzero denominators needed to prove composition, identity and inversion
(conversion_composes, conversion_identity, conversion_inverse). These
claims apply to each coordinate of a vector as a scalar quantity.
For derivative conversion, let x=s_L x’ and E’(x’)=E(s_L x’)/s_E. If E is differentiable at s_L x’, the chain rule gives
derivative_conversion proves the derivative statement from HasDerivAt;
force_conversion proves the sign/scaling identity. Physical unit scales must
be positive even though the derivative identity itself needs no positivity
hypothesis. This does not prove that a particular computed gradient is the
energy derivative. The present theorem is scalar; tensor stress and lattice
variation remain future work.
The numerical bridge checks molecule.from_xyz against an independent scale
adapter, and runs native water RHF plus the ASE calculator at the same geometry.
ASE energies are Hartree times eV/Hartree; forces are minus the gradient times
Hartree/Bohr. The XYZ Bohr-to-angstrom constant is 0.529177210903, while the tested
ASE dependency supplies 0.5291772105638411. Each interface is checked against its
actual convention. A cross-interface round trip must not silently assume these
constants are identical. No constants were changed in this stage.
DF contractions: finite sums and pair packing¶
Contractions.lean uses finite AO and
auxiliary index spaces. Let B have shape (n_aux,n_ao,n_ao) and D shape
(n_ao,n_ao), with real entries. Define the fitted ERI and Coulomb matrix by
Distributing products and interchanging finite sums gives
df_staged_equals_dense proves this full sum identity for arbitrary B and D.
No density symmetry, ERI permutation symmetry, positivity, spin factor or
conjugation is implicit. Complex k-point tensors need a separate specification.
This theorem starts with supplied factors; it does not bound auxiliary-basis
fitting error, verify metric regularization or prove integral accuracy.
For a logical n by m matrix, packPair uses offset i*m+j. The theorems
pack_pair_offset, unpack_pack_pair and pack_unpack_pair prove that offset
and both inverse laws, including empty index spaces. Physical Eigen/NumPy
strides are not logical index conventions. The native bridge tests explicit
packing rather than assuming an array’s raw memory already has that order.
The checked production paths are Python DensityFitting.build_J, resident C++
DensityFitting::build_J and contract_B_with in cpp/src/df.cpp. The independent
Python oracle uses explicit indexed products and math.fsum, with no einsum or
BLAS packing copied from those implementations. Synthetic nonsymmetric factors
expose permutations hidden by physical ERI symmetry. A real STO-3G water basis
with def2-universal-jkfit exercises native integral generation and both fitted
routes. Strided, Fortran-contiguous, symmetric and zero densities are included.
The two routes’ factors need not coincide, so each contraction is compared with
its own factors, followed by a comparison of the resulting Coulomb matrices.
Real arithmetic permits reassociation; floating-point arithmetic does not. The cancellation fixture has residual about 1.19e-7 under its explicit empirical 2e-6 absolute tolerance. This result is a regression check, not a roundoff bound. Streaming, screened, distributed and GPU DF paths are outside this milestone.
Orthogonal images: optimality and the complete tie set¶
Images.lean takes lattice vectors as columns
of A. For active orthogonal axes, assume A^T A=diag(L_i^2), with L_i>0. Decompose
a displacement into A f plus a perpendicular component. The theorem
orthogonal_cartesian_norm expands the Cartesian squared norm using the Gram
hypothesis. For an integer image n, its active part is therefore
The perpendicular squared distance is an additive constant, so it does not
change the minimizer (perpendicular_constant). With zero active dimensions,
molecular_image_sq gives zero for the periodic contribution, not zero for the
whole molecular distance. Orthogonal decomposition of arbitrary input data is
an interface assumption; the numerical fixtures construct it explicitly.
Theorem orthogonal_nearest_iff_axes proves that n globally minimizes q if and
only if each n_i is a nearest integer to f_i. The forward proof changes just one
coordinate; positive L_i lets it cancel that coordinate’s metric factor. The
reverse proof sums the coordinatewise inequalities. Integer spacing then gives
the complete characterization
nearest_axis_bound, integer_nearest_of_half and
orthogonal_nearest_iff_half_cell prove necessity and sufficiency. This includes
both integers at an exact half tie. Independent tied axes yield the familiar
face/edge/corner multiplicities; the cardinality formula itself is not a
separate Lean theorem. orthogonal_round_optimal proves optimality of one
rounded representative, which is insufficient for CCM boundary weights.
Translation f -> f+m maps the entire minimizer set by n -> n+m; reversal maps
(f,n) -> (-f,-n). image_translation, image_reversal, the corresponding
nearest_* implications and equivalences prove these properties.
nearest_value_unique proves equality of all minimizing squared distances,
without requiring a unique image representative.
The rational Python oracle enumerates per-axis floor/ceil candidates and their
Cartesian product. It treats a float as its exact binary rational value.
Checks cover interiors, faces, edges, corners, negative ties, displacements
several cells away, translation, reversal and Cartesian rotation. CCM’s
minimum_image accepts lattice rows; NEB’s _minimum_image_diff accepts columns.
The bridge transposes explicitly and checks NEB distances for 1, 2 and 3 active
axes plus unchanged perpendicular components. Exact sets are compared with
CCM at tol_bohr=0; a single rounded NEB representative is compared by distance.
Default CCM selection is tolerance-based. At displacement (2+1e-8,0,0) in a
cell of lengths (4,6,8), the exact set has one image and the production default
returns two. That distinction is recorded, not described as exact correctness.
A separate tolerance/topology contract remains necessary. The skew-cell
counterexample in the harness demonstrates why fractional component rounding
cannot be generalized to triclinic nearest images. Cutoff-neighbor completeness,
multi-k phases and general CCM topology are also outside these proofs.
Pair energies: multiplicities before normalization¶
Accounting.lean separates the pair domain from its contribution. For a finite directed domain P, assume an explicit bijection between P and R times a two-element direction index. If the two contributions for each representative r agree, then
pair_orbit_accounting proves this result by transporting the sum along the
bijection. Producing that bijection for a production pair list remains an
obligation of the bridge. For molecules, molecular_pair_accounting separately
proves that half the off-diagonal ordered sum equals the i<j sum under e_ij=e_ji.
For Coulomb interactions e_ij=Z_i Z_j/r_ij, distinct centers must have r_ij>0;
the algebraic theorem alone does not enforce physical nonsingularity.
For periodic labels (i,j,n), reversal is (j,i,-n). The reference ledger checks
closure, equality of reverse contributions and absence of fixed points. It
omits only (i,i,0). A synthetic nonzero self-image pair verifies that accounting
does not indiscriminately discard all i=j entries. The finite CCM nearest-image
model is not an infinite electrostatic sum or an Ewald specification.
For a nonempty selected image set S, uniform_boundary_weights proves
sum_{s in S} 1/|S|=1. This establishes weight normalization only. Completeness of
S, geometric multiplicity and tolerance-induced topology are separate issues.
The CCM numerical example uses a two-cell orthogonal supercell and compares
its nuclear energy with a complete independent weighted pair ledger.
Finally, normalization_sum and normalized_components prove distribution of
one common divisor across sums and across electronic/nuclear components.
Physical per-cell or per-particle reporting requires a positive count; Lean’s
totalized real division allows the algebraic identities even at zero, which
must not be interpreted as valid physical normalization. The native RHF check
also verifies E_total=E_electronic+E_nuclear. Embedding/double-counting terms and
energy derivatives across image or cutoff switches remain outside this stage.
No fixed-topology differentiability theorem is claimed for the pair engine.
Reproduction and observed evidence¶
Run the pinned proof audit with lake on PATH:
python formal/check_proofs.py
The audit builds all modules, checks the exact claim manifest and inspects every
theorem’s transitive axioms. All 58 claims use only the permitted standard base
(propext, Classical.choice, Quot.sound); none uses an admitted proof or a
compiler-trust axiom. The 30 new claims comprise 5 unit, 4 contraction, 16 image
and 5 accounting theorems. The manifest is
claims.json.
Supply formal/check_core.py, formal/core_reference.py and
formal/check_rhf.py to a calculation-scheduler job and run:
python check_core.py
Use a qualified vibe-qc Python installation with NumPy, ASE and the resident
native DF bindings. Do not use Python optimization that disables assertions.
The script writes core-checks.json with source-function hashes, native-library
hash, version, residuals, rejected invalid inputs and discriminating controls.
Reports contain environment paths; keep raw evidence outside public source.
The final native run passed on vibe-qc 0.18.6. The six inspected production Python function hashes matched the source functions under review; relevant native DF code was unchanged from that runtime’s source. This is source correspondence evidence, not proof of compiler or binary refinement.
Check |
Largest observed absolute residual |
|---|---|
XYZ coordinate conversion |
0 |
ASE energy / force conversion |
0 / 8.89e-16 |
Synthetic DF, ordinary fixture |
3.56e-15 |
Resident native DF, ordinary fixtures |
1.12e-15 |
Native pair packing |
2.23e-16 |
Orthogonal distances and weights |
0 |
Molecular and CCM nuclear pair ledgers |
0 |
Ten negative controls detect missing/double conversion, force sign/scale, DF permutation/spin factor, lost ties, skew rounding, missing pair half and double cell normalization. Nine invalid-contract checks reject incompatible quantities, nonfinite data, malformed DF shapes, a degenerate cell, a missing reverse pair and coincident distinct centers. These rejections belong to the reference adapters unless the harness explicitly calls a production boundary.
The complete conventional H2O/STO-3G RHF calculation gives -74.96033706904934 Hartree. The matching existing ORCA 6.1.1 result is -74.960337068945 Hartree, a difference of about 1.05e-10 Hartree. That independent reference supports the molecular fixture; it does not independently validate CCM’s physical model. The CCM nuclear ledger value is 2.2888543819998315 Hartree before dividing by two cells.
Provenance, limits and next stage¶
The formulation follows the explicit-assumption approach of Bobbin et al. (2024), sections 2-4. The separation of image rules, energy accounting and computational benchmarking is informed by Ugwuanyi et al. (2025), sections 4-5 and 7. The dimension/scale distinction is informed by Bobbin et al., sections 3-5, and explicit index spaces by Tooby-Smith, sections 2-3. The DF equations use Weigend (2002), equations 1-4; CCM boundary accounting is compared with the 2014 CCM formulation, equations 4-7, 11-13 and 18-24. Relevant primary texts were read for the design and reused here.
The original core follows MPL-2.0. Lean and mathlib remain pinned optional development dependencies with their existing upstream licenses. No external formal-physics project is vendored and no runtime dependency is added.
Next work should prioritize a quantitative floating-point contraction bound, certified general-cell image search, a separate tolerance/boundary policy, and fixed-topology energy/derivative contracts. Incremental adoption of quantity adapters at production interfaces can follow once their compatibility policy is agreed. This stage supplies a maintainable verified mathematical core and meaningful implementation checks; it does not establish a verified production engine or a priority claim over earlier formal chemical-physics work.