Quadrature cap warning: report the real change, name the integral, and test both

[susilehtola]

The adaptive quadrature's cap warning reported a relative change that was identically zero, by construction. This fixes that, makes what it reports testable, and makes the warning say which integral it is about — which turns out to matter, because the obvious reading of the warning pointed at the wrong integrand.

1. The reported change was always zero

In converge_block, the cap branch computed its relative change as (cur - prev) after prev = cur a few lines above — the difference of a matrix with itself. Every warning read:

  block magnitude 7.899e+00, relative change still 0.000e+00

which says "converged perfectly" directly under a sentence saying the opposite. It now reports prevdiff / scale, the last difference the convergence test itself saw. Where the cap is reached on the very first probe no comparison was ever made, and that now says so rather than also printing 0 — the old code conflated "no comparison" with "the change was small".

2. The warning was firing spuriously

With the number honest, the warning turned out to fire on diatomic-H2-hf-r, a CI case, reporting 1.2e-01. The source is kinetic()'s radial_integral(-1,0): the integrand 1/sinh(μ) diverges logarithmically at μ = 0, so the block magnitude grows with quadrature order. The divergence is confined to the row and column of the one basis function non-zero at μ = 0 — exactly the function remove_boundaries drops from every m ≠ 0 shell, the only consumer. The block is now judged without those entries (nskip), and returned whole.

H₂ is unchanged at −1.1336295702, the HF limit to all printed digits, and the warning no longer fires.

3. The label pointed at the wrong integrand

Every weight reported as a bare radial_integral, so a cap couldn't be traced to its integrand. That invites a specific misreading. quadrupole_zz() needs radial_integral(1,4) = sinh(μ)cosh⁴(μ), and a cap seen while building it naturally gets blamed on (1,4). But (1,4) converges to machine precision on the first comparison across the entire focal-separation range, even where the integrand spans eighteen orders of magnitude:

Rhalf	μ_max	(cosh μ_max)⁵	(1,4) block scale	rel at n=60
0.7	4.05	1.9e+07	1.8e+06	7.7e-16
0.005	8.99	1.0e+18	9.3e+16	3.4e-16

It's entire, and Gauss-Lobatto resolves it. The warning that actually fires is kinetic()'s (-1,0), whose magnitude (7.9 on H₂, 8.1 in the diatomic CAS test) is the same class as a cap of magnitude ~7 reported against a quadrupole build — and any Hamiltonian calls kinetic(). The warning now names the (m, n).

4. Making it testable

converge_block sat in an anonymous namespace in basis.cpp, and the numbers its warning reports were computed inside the one-shot print guard — unobservable after the first cap in a process, and unreachable from outside basis.cpp entirely. A test that merely checks the warning fires passes on the original bug.

It moves to converge_block.h (a template over the probe, so a header is its natural home) and gains an optional CapReport, filled on every call whether or not anything prints. converge_block_test drives it with probes known in closed form:

case	probe	asserts
never converges	every entry c·n	rel == (512−320)/512 = 0.375 exactly, and > 0
capped on first probe	seeded past 512	rel == −1: no comparison was made
converges	constant block	not capped, stops at n = 10
nskip	row/col 0 grow as log n	caps judged whole, converges with nskip = 1

It was verified to discriminate. With the old (cur − prev) computation reinstated it prints relative change still 0.000e+00 and three assertions fail.

Verification

converge_block_test   PASSED
diatomic-H2-hf-r      no warning, E = −1.1336295702 (unchanged)
ctest -L unit         10/10
ctest (full enabled)  60/60

🤖 Generated with Claude Code