demoAI

Two methods compute water's 6-31G integrals

Adam Zettel · September 24, 2026 · written with an agent

(μν ∣ λσ) = \iint \frac{χ _{μ} ( 1 ) χ _{ν} ( 1 ) χ _{λ} ( 2 ) χ _{σ} ( 2 )}{r _{12}} d 1 d 2

McMurchie-Davidson and Head-Gordon-Pople compute the same 4186 electron repulsion integrals for water in the 6-31G basis by different recurrences. They agree to rounding, each matches PySCF, and each one alone gives PySCF's Hartree-Fock energy. With a third implementation in the comparison, the cross-check names the implementation that disagrees and the integral where it does.

Open in the workbench Download .vivOutputs below are from the site build’s run; Run all runs it again in your browser.

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

The electron repulsion integrals are most of the cost of a Hartree-Fock calculation and the easiest place for a bug that does not announce itself. lib/hf computes them more than one way and compares the answers. This post takes the two methods the library cross-checks by default and runs each one alone on water in 6-31G: 13 basis functions and 4186 integrals after the eightfold permutation symmetry.

McMurchie-Davidson writes each product of two Gaussians as a sum of Hermite Gaussians. The integral becomes two matrix products over those expansions and the Hermite Coulomb integrals $R_{t uv}$ at the distance between the two product centers.

Head-Gordon-Pople starts from Obara and Saika's vertical recurrence, which builds up the angular momentum on one center of each pair, primitive by primitive. A horizontal recurrence then moves it onto the other center once per contracted shell quartet.

They share no recurrence. Both take the Boys function from lib/hf/boys, which is checked on its own against mpmath values (tools/boys_reference.py).

About 60 percent of the integrals come out as the same double from both methods, bit for bit: 2483 of the 4186. The exact count depends on the last bit of exp, so a different math library would change it. Every Viv engine, native or in the browser, uses the same one, so the count is the same in all of them, and a test checks it. The figure sorts the integrals that differ by how far apart they are, measured as $∣ a - b ∣/ (1 + ∣ a ∣)$ , the measure the library's cross-check uses. None is as far as $1 0^{- 15}$ . For the integrals larger than 0.01 the gap is at most 5.75 times the double-precision epsilon relative to the integral itself. The screened McMurchie-Davidson code skips any shell quartet whose Schwarz bound is under $1 0^{- 14}$ , so a few integrals no larger than $3 \times 1 0^{- 20}$ are exactly zero in its set.

Each method is also checked against PySCF's stored integrals within $1 0^{- 12}$ , and each one alone, fed to the self-consistent field with the default implementation of every other step, gives PySCF's energy $- 75.98394849810573$ within $1 0^{- 9}$ .

Finding which implementation is wrong. The cross-check takes the median of the implementations' answers, entry by entry, and blames the implementation farthest from it, at a named integral. For an even count the program takes the upper of the two middle values as the median. With two implementations that is the larger value, so a disagreement is detected but the blame falls on whichever one gave the lower value. With three, the median is a correct value. The last part of the program adds a third implementation that returns the Head-Gordon-Pople integrals with integral 2000, which is $- 0.117$ , multiplied by $1 - 1 0^{- 9}$ . That makes it larger, so with only two implementations the median is the wrong value and the correct one takes the blame. With three, the check names the altered implementation and integral 2000.

Printed

E_mcmurchie_davidson: -75.9839

E_head_gordon_pople: -75.9839

Integrals where McMurchie-Davidson and Head-Gordon-Pople differ, by size of the difference

the integrals that are not bit-identical

Tests 8 of 8 pass

each method alone matches PySCF's 4186 integrals within 1e-12, and the two agree to 1e-15 (passes)
2483 integrals are bit-identical; the rest are within 5.75 epsilon or screened below 3e-20 (passes)
each method alone gives PySCF's RHF energy -75.98394849810573 within 1e-9 (passes)
with three implementations the check names the altered one and integral 2000; with two it detects the error but blames the correct one (passes)
water in 6-31G has 13 basis functions and 4186 integrals after the eightfold symmetry (passes)
about 60 percent of the integrals are bit-identical from the two methods (passes)
the screened McMurchie-Davidson code skips quartets whose Schwarz bound is under 1e-14 (passes)
the integral the wrong implementation alters, integral 2000, is -0.117 (passes)