vibe-qc day 2: closing out Ewald and laying the symmetry foundation
vibe-qc is an open-source quantum chemistry and solid-state code with a C++17 numerical backend, a Python frontend via pybind11, and an ASE Calculator interface for geometry optimization and structure I/O. It is being written in the open, one coding session at a time, with Claude Opus as the implementation engine and me steering at the architecture level. Day 1 shipped the molecular HF/DFT/MP2 stack (v0.1.0, April 18) and then pushed into periodic territory: one-electron lattice integrals, Gamma-only and multi-k RHF, multi-k KS-DFT, and the first Ewald infrastructure through Madelung-cancellation helpers. Day 2 closed out the v0.2.0 Ewald milestone and planted the first stake for v0.2.5 space-group symmetry. Four deliverables, 26 new tests (720 to 746 total), Sphinx docs building clean.
Phase 12e-c-4: end-to-end EWALD_3D dispatch
The previous session had the Ewald machinery in place but not yet wired to the user-facing entry points. Today that changed. Both run_rhf_periodic_scf(...) and run_rhf_periodic_gamma_scf(...) now branch on options.lattice_opts.coulomb_method: pass DIRECT_TRUNCATED and you get the existing C++ driver unchanged; pass EWALD_3D and you get the new Python driver, with nuclear_repulsion_per_cell automatically rerouted through Ewald summation as well. That last part matters: an Ewald electronic energy paired with a direct-truncation nuclear repulsion would be self-consistent in neither the long-range cancellation nor the ω-dependence, and the resulting total energy would be meaningless. Electronic and nuclear sides now always agree on which summation scheme is in use.
The benchmark suite that shipped alongside covers H₂ crystal, LiH rocksalt, MgO rocksalt, and Ne FCC — 11 tests in total, checking convergence, $\omega$-invariance (energy independence of the Ewald splitting parameter), an equation-of-state scan, and the equivalence of a [1,1,1] multi-k mesh with a Gamma-only calculation. All pass.
One finding worth documenting: atom positions relative to the FFT Poisson grid origin are not numerically irrelevant. Atoms sitting at box corners, rather than centered in the cell, inflate the $\omega$-invariance residual by roughly two orders of magnitude. The physics is correct either way, but the numerical conditioning is not. This is now flagged in the docs and the benchmark suite uses centered geometries as the reference.
Phase 12e-c-4c-iii-c: multi-k Pulay DIIS
The multi-k Ewald driver shipped without DIIS. For loose cells and [1,1,1] meshes that was acceptable — pure density damping converged, slowly. For tight cells with denser k-meshes it was not. An H₂ chain at a = 10 bohr with a [2,2,2] mesh simply plateaued under damping and never converged.
The fix is Pulay DIIS extended to the Brillouin zone. The commutator error vector at each k-point is:
$$\mathbf{e}(\mathbf{k}) = \mathbf{F}(\mathbf{k})\mathbf{D}(\mathbf{k})\mathbf{S}(\mathbf{k}) – \mathbf{S}(\mathbf{k})\mathbf{D}(\mathbf{k})\mathbf{F}(\mathbf{k})$$
and the Pulay B matrix couples iterations i and j with a k-weighted inner product:
$$B_{ij} = \sum_{\mathbf{k}} w_{\mathbf{k}}\, \mathrm{Re}\,\mathrm{tr}\!\left[\mathbf{e}_i(\mathbf{k})^\dagger\, \mathbf{e}_j(\mathbf{k})\right]$$
The k-weights are the standard Monkhorst-Pack weights that already appear in the energy and density assembly, so no new parameters appear. The DIIS coefficients solve the same linear system as in the molecular case; the only change is that the scalar error measure is now a sum over the zone rather than a single matrix norm.
The numbers are clean. H₂ chain, [1,1,1]: 13 iterations under damping, 7 with DIIS. H₂ chain, [2,2,2]: previously non-converging under damping, now 4 iterations with DIIS. This is the same qualitative behavior the molecular DIIS delivered on H₂O back on day 1 (52 iterations to 9), just extended to k-space.
Phase 12f: periodic Becke partition for tight-cell DFT
The molecular Becke fuzzy-cell partition (J. Chem. Phys. 88, 2547, 1988) assigns each grid point a weight that partitions unity across atomic cells using a smooth step function of interatomic distances. In a periodic system with a tight unit cell, image atoms from neighboring cells fall within the smoothing radius of the step function, and ignoring them breaks the partition: the weights no longer sum to the cell volume.
The fix is to extend the partition denominator over home atoms plus image atoms within a user-specified radius. The implementation is wired into run_rks_periodic via two new options: PeriodicKSOptions.use_periodic_becke (boolean, default False) and becke_image_radius_bohr (float). Default behavior is unchanged from v0.1.0; the periodic extension is an explicit opt-in for users who know they are working in a tight cell.
The empirical witness is a 5-bohr cubic H₂ cell. Without periodic Becke, the molecular partition produces a total grid weight of 18,526 — wildly wrong for a cell with volume 125 bohr³. With periodic Becke, the total grid weight comes out as 126, matching V_cell = 125 bohr³ to within the grid quadrature error. The discrepancy between the two modes (18,526 vs 126) is not a subtle numerical issue; it is a categorical correctness failure of the molecular partition in tight periodic geometry. Any DFT calculation on a tight crystal without this correction is integrating a density that does not integrate to the correct electron count per cell.
Phase SYM3a: orbit-reduced storage for lattice integrals
The v0.2.5 symmetry milestone starts here. The one-electron lattice integrals $S(\mathbf{g})$, $T(\mathbf{g})$, $V(\mathbf{g})$ are indexed by a cell offset vector $\mathbf{g} = (g_1, g_2, g_3)$. For a system with space-group symmetry, many (g, atom-pair) combinations are related by point-group operations and carry the same numerical content up to a known rotation. SYM3a exploits this to compress storage.
The machinery builds on SYM2c, which already identifies atom-pair and cell-index orbits for arbitrary (including non-origin-fixed) structures like NaCl. SYM3a takes those orbits and stores one representative sub-block of shape $(n_{\mathrm{AO},a},\, n_{\mathrm{AO},b})$ per orbit, rather than a full $(n_{\mathrm{bf}} \times n_{\mathrm{bf}})$ block at every cell triple. A round-trip compress/reconstruct test on simple cubic Pm-3m He at STO-3G, with a 10 bohr cutoff, partitions 33 cell triples into 5 orbits and achieves 6.6× memory reduction. Reconstruction is exact to 10⁻¹⁶.
What SYM3a does not yet do is reduce compute: the integral kernels still evaluate every shell-pair times cell-triple combination before handing off to the compressor. That is SYM3b, the next item on the list. SYM3b is where the |G|-fold reduction hits wall-clock time rather than just peak memory, because the kernel skips shell-pair/cell-triple combinations that are equivalent under the point group instead of evaluating and then discarding them. The storage handle SYM3a provides is what SYM3b will plug into.
Where v0.2.0 stands
The v0.2.0 contract is “ω-invariant 3D periodic SCF that reduces to molecular HF in the loose-box limit and integrates the unit-cell volume correctly in tight cells.” As of today that contract is end-to-end testable with two option flags:
options.lattice_opts.coulomb_method = EWALD_3D
options.ks_opts.use_periodic_becke = TrueInternal-consistency benchmarks are in place across four crystal systems. The remaining gating item for the v0.2.0 tag is a cross-check pass against published CRYSTAL reference energies on LiH, NaCl, MgO, and Si at published geometries. That is a validation exercise against an external reference rather than new implementation work, and it is next on the list.
What is next
Immediately: SYM3b, the kernel-level compute reduction that makes the |G|-fold saving from space-group symmetry show up in wall-clock time. The substrate — Wigner D-matrices (SYM1), AO permutation matrices and orbit identification (SYM2a/b), atom-pair-resolved orbits for non-origin-fixed structures (SYM2c), and now the compressed storage layer (SYM3a) — is all in place. SYM3b is the payoff.
Further out, the headline feature of the entire project remains the cyclic cluster model. The roadmap runs from v1.0 (feature-complete molecular and periodic HF/DFT/MP2) through v2.0 (HF-CCM in 3D) and on through a sequence of independently publishable steps: MP2-CCM, local-MP2-CCM, CCSD-CCM, and eventually projection-based embedding where a CCM-correlated central cluster sits inside a DFT-treated periodic environment. The capstone is multi-scale defect chemistry at chemical accuracy on systems with hundreds of atoms in the cluster. Every piece shipped in the first two days of this project is load-bearing for that goal.
The code is at vibe-qc.com. Licensed MPL 2.0.