Add OrbitalEvaluator and Gaussian cube writer utility

[ConradJohnston]

Summary

Adds OrbitalEvaluator and write_cube / CubeGrid to GauXC. The evaluator wraps the host collocation driver with thread-parallel batched evaluation and per-batch shell screening; the cube writer is self-contained and has no kernel dependency. Either can be used independently. Also adds write_cube_hdf5 behind GAUXC_HAS_HDF5.

API sketch

auto eval = OrbitalEvaluatorFactory::make_orbital_evaluator(
    ExecutionSpace::Host, basis);   // optional 3rd arg: screening tolerance

auto grid = CubeGrid::from_molecule(mol, 80, 80, 80);
eval.eval_orbital(grid, C, out);        // or (npts, points, C, out)
eval.eval_density(grid, D, ldd, out);   // or (npts, points, D, ldd, out)

write_cube(path, mol, grid, field);
write_cube_hdf5(path, mol, grid, field);  // if GAUXC_HAS_HDF5

Performance

Versus PySCF cubegen 2.12.1, cc-pVDZ, 16 threads on an AMD Ryzen AI 7 PRO 350 under WSL2. Single-threaded BLAS on both sides, grid evaluation only — no SCF, no file I/O. The harness parses geometries out of tests/standards.cxx so the two codes cannot drift apart, and basis size is asserted identical for every row.

system	nbf	grid	GauXC orbital	vs cubegen	GauXC density	vs cubegen
water	24	80³	18 ms	3.2×	15 ms	9.2×
benzene	114	80³	45 ms	3.4×	68 ms	10.9×
taxol	1099	40³	34 ms	11×	64 ms	63×
taxol	1099	80³	59 ms	56×	204 ms	155×

The advantage is per-batch shell screening. cubegen evaluates the full basis at every point — verified directly: eval_gto 500 Bohr from the molecule, where every AO is identically zero, costs the same as on top of it. So its cost is O(npts · nbf) for orbitals and O(npts · nbf²) for density, while GauXC's tracks the locally surviving nbe. Concretely: taxol at 80³ has 10× benzene's basis but takes only 1.3× longer.

Grids are walked in spatially compact tiles rather than contiguous index ranges, since a contiguous range spans the full z extent as soon as it crosses a row boundary and screens poorly. Tiling is the only grid traversal: counting the collocation and GEMM work each one actually does, tiling never lost anywhere measured — benzene and taxol, 64³ to 200³, margins of 3 and 6 Bohr, and a single-plane grid where its scatter degenerates to runs of one point — ranging from a 2–3% saving on benzene at 128³ to 53% on taxol at 200³.

Density is compared against a GEMM-contracted cubegen. As shipped it contracts with numpy.einsum('pi,ij,pj->p', ao, dm, ao), which costs a further 15× on taxol and 5.6× on benzene; against that, taxol at 40³ is ~970×.

Accuracy against an unscreened PySCF reference evaluated at identical points (taxol, 40³, 1e-10 screening tolerance): density agrees to 4.9e-16 L2 relative, orbitals to 1.0e-10 absolute on a field spanning [-0.71, 0.32].

On the timing methodology

Wall-clock timing on this host is unreliable at the tens-of-milliseconds scale and simply taking more repeats does not fix it. Over 30 independent runs of each code, individual timings spanned a factor of 20 (cubegen water orbital: 27 ms to 545 ms, median 143 ms), and the running minimum was still falling at N=30 rather than converging — under WSL2 the available vCPU capacity drifts with host activity, so there is no fixed floor for a minimum to converge to. Bootstrapped spread on a 30-sample minimum ranged from 2% to 188% depending on the quantity.

The water and benzene rows are therefore the median of 25 tightly interleaved pairs, each code measured immediately after the other so both see the same machine state and the drift cancels in the ratio. Interquartile ranges: water 2.2–4.1× (orbital) and 7.5–12.6× (density); benzene 3.1–4.1× and 9.9–12.2×. Quoted GauXC times are the fastest observed over those runs. The taxol rows are best of 3 and their margins are orders of magnitude outside any of this.

Tests

77,123 assertions across 22 cases for [orbital_evaluator] and [cube]; 3,640,979 across 56 for the full suite. Covers the evaluator against a direct eval_collocation reference, the screened paths (nbe == 0 and 0 < nbe < nbf), multi-batch grids, the tile bounding boxes and the scatter they require, multi-orbital scatter into a padded output, single-plane grids, accessors and move semantics, every argument guard and empty-work no-op, a dense density matrix through a padded leading dimension, invariance of the results to the thread count, safety of concurrent const invocation, CubeGrid overload consistency, cube and HDF5 round-trips, multi-chunk cube writes, and %13.5E formatter edge cases including NaN, infinities, subnormals and 3-digit exponents.

Known limitations

• The CubeGrid overloads are not bitwise equal to the pointer overloads given the same points. A grid is walked in spatially compact tiles, which screen better, whereas a caller-supplied point array is batched in the order it arrives, so the two drop different shells and agree to the order of the screening tolerance rather than bitwise. Results are otherwise bitwise reproducible, including across thread counts, which is pinned by a test.
• Screening tolerance is a constructor argument applied to the evaluator's private basis copy; it defaults to GauXC's 1e-10. Orbital error tracks it directly, density error goes as its square, so density tolerates a looser setting (1e-6 is ~2.5× faster on density with ~1e-12 error).
• format_e13_5 matches glibc %13.5E except for values within a few ulp of a rounding boundary in the 6th significant digit, where scaling the mantissa can tip it across; the two then differ by one unit in the last printed digit, never more. Fuzzed against snprintf: 3M random bit patterns were byte-identical, and of 3M values deliberately parked within 3 ulp of a 6th-digit tie, 2.7% differed, all by exactly one in the last digit.
• ldo and ldd above the range BLAS is handed are rejected, since BLAS takes leading dimensions as int.

Follow-ups

Device backends via the existing CUDA/HIP collocation drivers; eval_orbital_gradient / eval_density_gradient; native-layout collocation to avoid Gau2Grid's transpose; NUMA-aware first touch.

Closes #128