Kernel cclib GaussGrid

Kernel: cclib GaussGrid

cclib-mojo accelerates cclib's electron-density / wavefunction-on-grid path (cclib.method.volume: Volume.wavefunction() / electrondensity()), whose PyQuante backend evaluates contracted Gaussian basis functions one grid point at a time in pure Python. Package: python/cclib_mojo · Kernel: kernels/gaussgrid · pip install cclib-mojo

The algorithm

A molecular orbital on a grid is a sum over contracted Cartesian Gaussian primitives, ψ(r) = Σ cᵢ·Nᵢ·xᵃyᵇzᶜ·exp(-αᵢr²), evaluated at every grid point. The reference path interprets that formula per point. The kernel restructures it:

1. Flatten (Python, shared by both backends): validate gbasis/geometry, compute normalization constants with PyQuante's exact formulas and operation order (THO eq. 2.2 and 2.12), and flatten everything into typed arrays.
2. Separable-axis precompute (kernel): the Gaussian factorizes exactly — exp(-α·dx²)·exp(-α·dy²)·exp(-α·dz²) — so the per-axis factors are computed once per call, not per point.
3. One SIMD FMA per (point, primitive): the inner loop is a single float64 fused multiply-add in compiled code, with the grid row under update kept cache-resident. One FFI call evaluates the whole 3-D grid for one MO (or accumulates the density over all MOs) — FFI cost per call, not per point.

Autopsy: why the reference is slow (measured, two honest baselines)

The PyQuante-1 path interprets ~5 Python operations per (grid point × basis function × primitive) — hundreds of millions of interpreter steps for a production cube grid. cclib's hot loop on that backend (pyamp() → per-point CGBF.amp(x, y, z)) is Python-2-era code that cannot run on Python 3, so the benchmark drives a verbatim Python-3 transcription of that amplitude path (tests/pyquante1_oracle.py, same formulas, same operation order, same per-point interpreter cost profile) exactly the way cclib drives it.

But the headline is not "7,000× vs ancient code" — there are two measured baselines on the same machine and workload:

• vs the PyQuante1 path: 1,898–7,033× (the table below).
• vs cclib's pyquante2 backend (NumPy-vectorized cgbf.mesh, the fastest backend cclib ships today): cclib's real Volume.wavefunction() takes 0.689 s for the 50³ single-MO workload vs 0.0095 s for cclib-mojo — 72×, max abs diff 8.4e-13. The remaining gap over the PyQuante1 table is pyquante2's NumPy vectorization, which still pays per-(basis-function, grid) temporaries; the kernel pays one SIMD FMA per (point, primitive).

Measured numbers

Benzene, 6-31G* (12 atoms, 102 contracted Cartesian basis functions, 192 primitives — basis data from PyQuante 1.6.5's published basis library), seeded MO coefficients, median of 5, single-threaded. Apple M4 Max, macOS 26.6.2 arm64, Python 3.12.14, NumPy 2.5.3, Mojo 1.1.0. Reproduce: pixi run bench-cclib.

workload	PyQuante1 path (s)	cclib-mojo (s)	speedup
wavefunction 1 MO, 50³	32.23	0.0126	2,567×
wavefunction 1 MO, 100³	345.37	0.0491	7,033×
density 3 MOs, 50³	63.59	0.0335	1,898×
wavefunction 1 MO, 50³ — NumPy fallback	32.23	0.1691	191×

Setup is per-call on both sides (cclib rebuilds PyQuante basis objects each wavefunction() call; the wrapper rebuilds its flat arrays): 2.09 ms vs 1.48 ms.

Parity proof

• Gate: element-wise vs the PyQuante1 reference within 1e-10 relative, asserted before every timing run; measured agreement is 1.5e-12 (vs pyquante2: 8.4e-13) — the values are the same numbers.
• Differential suite (19 tests per backend, 38 per full run — native + forced fallback): s/p/d/f shells (STO-3G H₂O, 6-31G* d on O, cc-pVTZ f on C), multiple atoms, MO amplitude vs summed density, non-cubic anisotropic grids, grids offset from the origin and from all atoms, exact-zero coefficients, grid-point ordering, validation errors — plus end-to-end equivalence with cclib's actual Volume.wavefunction()/electrondensity() on the pyquante2 backend.
• Conventions verified against cclib, not assumed: centers Å→bohr with cclib's convertor constant (×1.8897261245), grid axes ÷0.5291772109 (the two are not exact reciprocals — we mirror cclib exactly), shell expansion order of cclib's sym2powerlist (S/P/D/F), Volume.data C-ordering (x outer, z inner), and the abs(coeff) > 0.0 inclusion rule. Both backends share basis flattening and normalization (_basis.py), so they cannot disagree about constants.

Options matrix

Surface	Status
wavefunction_on_grid(gbasis, atomcoords, mocoeffs_1d, origin=, step=, shape=)	standalone, cclib Volume.data layout
density_on_grid(gbasis, atomcoords, mocoeffs_nmo, …)	Σ|ψ|² over MOs (×2 for closed shell)
cclib_mojo.cclib_integration	drop-in mirrors of cclib.method.volume functions — no monkeypatching; parse with cclib, evaluate here, write cube files with cclib
Shells S/P/D/F (Cartesian)	exactly cclib's sym2powerlist scope
G shells, spherical-harmonic (5d/7f)	not supported — cclib's volume path doesn't handle them either
Threading	single-threaded (determinism first); multithreading is future work
Dependencies	NumPy only — cclib is not required (cclib_integration duck-types ccdata/Volume)
CCLIB_MOJO_NATIVE_LIB, CCLIB_MOJO_DISABLE_NATIVE=1, backend_info(), native_available()	standard controls
CCLIB_MOJO_ALLOW_PURE_WHEEL=1	build a py3-none-any fallback wheel (e.g. for Windows)

Runnable example: quickstart.py at the repo root (water, STO-3G). Citing: please cite both cclib (O'Boyle et al., J. Comput. Chem. 29 (2008) 839–845) and this package — the kernel is a clean-room implementation of the textbook contracted-Gaussian formulas (Taketa, Huzinaga, O-ohata, J. Phys. Soc. Jap. 21, 2313 (1966)).

---

Benchmarks · Kernels · How It Works · Apache-2.0, © 2026 Algenta