Potential factory — fitting model surfaces to target curves¶
Location: projects/potential_factory/ (the fitter: target, extract,
tracker, fit, report), qscat.model.flexible (the promoted ansatz),
validation/factory/ (the round trips, the O₂ campaign and its results).
Origin: docs/superpowers/specs/2026-08-24-potential-factory-design.md.
Status: the fitter is toy-stage; the ansatz and the fitted qscat.model.O2
(+ O2_SO12/O2_SO32) are promoted.
Relates to: docs/physics/potential-factory-options.md (why the factory
exists, the target-tier survey it is built from, and the decisions taken on
it), docs/physics/nonlocal-resonance-model.md (the T3 forward model
— the same AsymptoticDiscreteState/v_dk_plus/gamma_from_coupling chain
this package fits against), docs/physics/n2-resonance.md (the published
E_res(R_0)/\(\Gamma(R_0)\) values this package’s pole solver reproduces).
Units: atomic units.
Key result¶
The factory fits a FlexibleDiatomicModel — EMO neutral curve, Gaussian well
with \(\lambda(R)\) and \(\alpha(R)\), optional repulsive shell — to a tiered
Target (T0 neutral curve, T1 pole curves, asymptote, T3 energy-dependent
width) and stops-and-reports per tier. It round-trips the three published
models’ own curves to 1e-3–1e-4 relative (N₂ 9.9e-4, NO 2.1e-3, F₂ 1.5e-3),
and the observable-level check on the emoscat decks reproduces the
published models’ VE/DA cross sections to 1e-9. Its first real-molecule
target, O₂ from the vector-extracted figures of Alt & Houfek (2021), is met
at the extraction floor over the full 1.85–6 bohr range (T1 rms 20 meV,
\(\Gamma\) 8 %/14 %, asymptote to −0.08 mHa at 14 bohr; 87 s on a laptop) once
\(\lambda(R)\) takes the long-range-correct TailR form, and the exact 2-D VE
cross section of the fitted potential — spin–orbit components summed at ⅓
each, 3343 energies × six channels each with MUMPS — sits on the paper’s own
nonlocal-model doublets in all six 0 → v’ panels (positions within 1–8 meV,
heights 0.9–1.1 on the main peaks, doublet separation 19.0–19.3 meV). The
metric that predicted that figure was not the curve rms but the anion’s
quasi-bound levels (±7 meV over v = 0…24); the open item is T3, which is not
discrete-state-consistent. Nothing is fitted to or compared with experiment.
What it does¶
Every model in qscat.model (N₂, NO, F₂) is a hand-tuned testbed: a Morse-like
neutral curve plus a Gaussian electron–molecule well \(V_\mathrm{int}(r,R) = -\lambda(R)\,
e^{-\alpha_c r^2}\) with a single free function of R (\(\lambda\)) and a constant \(\alpha_c\).
projects/potential_factory/ builds a strictly richer ansatz,
FlexibleDiatomicModel (ansatz.py), and a tiered Target format
(target.py) plus an extractor (extract.py) and a staged fitter
(fit.py/report.py) that fits the ansatz to a Target.
The ansatz has three pieces:
Neutral curve. An Expanded Morse Oscillator (EMO, Le Roy’s form) in the dimensionless radial variable
y_p(R) = (R^p − R_e^p)/(R^p + R_e^p): \(v_0(R) = D_e[(1 - e^{-\beta(R)(R - R_e)})^2 - 1]\) with \(\beta(R) = \sum_i \beta_i\, y_p(R)^i\). One \(\beta\) reduces this exactly to a Morse curve.Interaction. A Gaussian well whose depth and range are both smooth functions of
R: \(V_\mathrm{int}(r,R) = -\lambda(R)\,e^{-\alpha(R) r^2}\), where \(\lambda(R)\) and \(\alpha(R)\) are each aSmoothR— a sigmoidf_inf + f_0/(1 + e^{f_1(R − R_f)})times an optional polynomial correction iny_p(R). Withcoeffs == ()this is exactly Houfek’s own sigmoid form for \(\lambda(R)\); the existing models fix \(\alpha\) constant, which is the special case \(\alpha\).f_0\(= 0\).Optional shell. A repulsive Gaussian barrier
shell(R)\(e^{-\alpha_b(r - r_b)^2}\) added to the well, used only when the T3 (coupling) tier needs to reshape the off-resonance background that \(\lambda(R)\)/\(\alpha(R)\) alone cannot reach.
from_diatomic embeds an existing DiatomicResonanceModel (N₂/NO/F₂) into
this family exactly — \(\lambda(R)\)’s sigmoid constants and \(\alpha(R)\)’s constant term
are read off the published model’s own lambda_inf/lambda_1/R_lambda/
alpha_c, so every published model is a point in FlexibleDiatomicModel’s
parameter space with zero embedding error. That is what makes the round-trip
oracle below meaningful: fitting a perturbed seed back to the exact
published constants is a real test of the fitter, not a tautology.
Target (target.py) carries three independent tiers, each optional and
each fitted in order:
T0 — neutral curve.
NeutralTarget: either a table (Curve.from_table, a cubic spline) or the constantsR_e,D_e, \(\omega_e\), \(\omega_e x_e\) directly.T1 — resonance curve.
ResonanceTarget:V_ion(R) = V_0(R) + E_res(R)and \(\Gamma(R)\) asCurves over the range where a pole was actually found, plus the electron affinityEA(theR → ∞asymptote).T3 — energy-dependent width.
CouplingTarget: the nonlocal model’s own \(\tilde\Gamma(\varepsilon,R) = 2\pi\, \varepsilon^\alpha\, |V_{dk}^+(\varepsilon,R)|^2\) — either Alt & Houfek’s closed-form Eq. (25)–(27) (from_alt_houfek) or a table (from_table, log-space interpolated — see Limitations).
extract_target builds a Target from an existing DiatomicResonanceModel’s
own computed curves — T0 from v0, T1 from a per-node gated pole walk
(walk_t1, below), T3 from model_gamma_tilde sampled on a grid — which is
what turns round-tripping into a differential test: fit the factory back to
data generated by the very model it should reproduce.
The forward models and the gradient¶
ElectronicPair.pole (tracker.py) is the gated fixed-R pole finder.
It holds two electronic FEM-DVR-ECS grids that differ only in ECS angle
(default 35°/44°), each with its kinetic-energy matrix built once; .pole(v_fn, window) calls qscat.ecs.find_resonance_pole on the two grids’ Hamiltonians
and gates the match: a bound candidate is accepted only if Re E < 0 and
|Im E| is at round-off; a resonant candidate only if its two-angle residual
clears both an absolute tolerance and residual < gate_frac · \(\Gamma\) (default 5%
of the width). This rejects angle-stable matches that are numerical
coincidences rather than physical poles — a fake ℓ = 1 well’s spike matches
across the two angles with residual 3.9e-3, but that is 2.7× its own
\(0.05\cdot\Gamma\) gate and is correctly rejected, not just caught as an exception. The
package’s convention throughout is the bare well: every window and every
target the fitters build is on the electronic-shift energy scale (`v_int(r,R)
ℓ(ℓ+1)/2r²
, withoutv0(R)), sincev0(R)` is a separate, additive term that the pole solver never needs to see.
pole_sensitivity(H, E) \(= \psi_i^2\) is the exact gradient of a pole energy
with respect to every potential value at once. For a complex-symmetric
Hamiltonian the Hellmann–Feynman theorem holds under the bilinear c-product
(qscat.linalg.c_product, \(\psi^T\psi\) not \(\psi^\dagger\psi\)) rather than the usual Hermitian
inner product: \(\partial E/\partial V_i = \psi_i^2\) for the c-normalized eigenvector \(\psi\) (\(\psi^T\psi =
1\)). This is verified against a finite difference referenced to the
Hamiltonian’s own matched eigenvalue (not the two-angle midpoint
find_resonance_pole returns, which sits ~residual/2 away from either
grid’s true eigenvalue and would otherwise leak a spurious, non-vanishing
term into the difference quotient): at h = 1e-5 the absolute error is
2.22e-7, relative 2.24e-6 — genuine round-off, not a sign or conjugation
defect.
solve_pole_params is damped Newton in \((\lambda, \alpha)\) on a target pole energy,
built entirely from that one gradient: \(\partial E/\partial\lambda = \sum_i \psi_i^2\,(-e^{-\alpha r_i^2})\),
\(\partial E/\partial\alpha = \sum_i \psi_i^2\,(\lambda r_i^2 e^{-\alpha r_i^2})\). A resonant target solves the full
complex 2×2 system; a bound target has one real equation in one unknown
(the true energy has no width to constrain a second parameter), so only \(\lambda\)
moves and \(\alpha\) is held at whatever the seed carried — a single real bound-state
energy genuinely cannot pin two parameters. Recovering N₂’s true well
parameters converges cleanly from a perturbed seed in both regimes: a
resonant target seeded at \((\lambda, \alpha) = (4.5, 0.5)\) against the true (4.9624, 0.4) (λ off by ~9.3%, α by 25%) converges in 8 Newton iterations with no
damping ever triggered; a bound target seeded ~8.7% off in \(\lambda\) (\(\alpha\) fixed,
per the single-equation reduction above) converges in 4. Two structural
difficulties surface only
once continuation is layered on top of this solver (below): a classification
mismatch when the seed’s pole and the target’s class disagree (the seed is
bound but the target is resonant, or vice versa — “bridging”), handled by a
reduced real-part-only step and acceptance test while the mismatch persists;
and the branch point itself, where the bound window’s search range is
capped at Re E < 0 and structurally cannot see a pole that has already
crossed into the resonant regime — bridging therefore searches a window
spanning the union of both regimes. Every candidate acceptance additionally
requires continuity: |E_candidate − E_current| ≤ MAX_STEP (0.15 Ha), so a
more angle-stable but unrelated pole elsewhere in the wide bridging window can
never be silently substituted for the one actually being tracked.
track_curve continues solve_pole_params over an R grid, seeding each
node from the previous node’s converged \((\lambda, \alpha)\). A node whose Newton solve
fails or whose target is undefined (NaN) is dropped, not treated as an
error: its \((\lambda, \alpha)\) are carried forward unchanged so the walk can resume past
it. On N₂’s R = 3.0 → 1.6 continuation (15 nodes crossing the bound/resonant
threshold), 14/15 converge to a maximum relative error of 1.9e-7; the one
dropped node sits inside the ~0.02-bohr-wide slice around R ≈ 2.4 where the
true pole’s real part passes through (0, e_floor) — too positive to gate as
bound, too small to clear the resonant energy floor — a genuine dead zone of
the angle-stability gate itself, not a tracker defect.
fit_neutral’s T0 stage fits the EMO curve at the target table’s own
nodes (a CubicSpline passes exactly through its input data, so an off-node
probe grid would inherit its \(O(h^4)\) interpolation error as a systematic
residual floor), then checks the vibrational ladder the fitted curve
implies: \(\varepsilon_1 - \varepsilon_0\) from qscat.core.vibrational.vibrational_states on
qscat.core.grids.nuclear_grid(), against the target’s own. "met" requires
\(\lvert\Delta\omega_e\rvert / \omega_e \le\) tol.omega_e_rel as well as the curve residual, and the
comparison always appears in TierResult.detail.
Two details make that comparison honest. First, a spectroscopic table quotes
the harmonic constant \(\omega_e\), which is not the 0 → 1 spacing; for the
Morse curve this branch builds from those constants the exact spacing is
\(G(1) - G(0) = \omega_e - 2\,\omega_e x_e\) with \(\omega_e x_e = \omega_e^2 / 4D_e\). Comparing
against the bare \(\omega_e\) would report a fixed anharmonic offset as fit error —
0.83 % (N₂), 1.96 % (NO), 3.33 % (F₂), two of the three over the 1 % default
for a curve reproduced to 1e-13. Second, a table-backed target is
NaN outside its own span and real-argument only, while the nuclear grid runs
to 40 bohr with a complex ECS tail; the reference solve therefore fills the
tail and any out-of-table real point with the fitted curve, and that fill is
admissible only if the two levels never reach it. That is checked, not
assumed: if more than 1e-8 of either level’s probability sits outside the
table, T0 reports omega_e: not checked with the leak rather than quoting a
meaningless number (measured leak on the round-trip tables: 6e-16 N₂,
2e-15 NO, 7e-15 F₂).
fit_resonance’s T1 pipeline chains these into a staged fit. Its R
nodes are the target’s own table nodes when it has them (see Limitations).
Each node is first classified bound or resonant by a Γ floor
(\(\Gamma_{\mathrm{target}} <\) tol.gamma_floor, the same threshold used to pick the
resonant-node mask below) rather than an exact-zero test, matching what the
Newton bridging logic above already treats as “genuinely bound.” The nodes
are tracked once; \(\alpha(R)\) is then fit from the resonant nodes only,
Γ-weighted (weight ∝ \(\Gamma_{\mathrm{target}}\), normalized to max 1) — a bound node’s single
real eigenvalue cannot determine \(\alpha\) at all (\(\partial E/\partial\alpha\) never enters a
lam-only Newton step there), and the near-degeneracy \(\partial E/\partial\alpha \approx -6.06\cdot\partial E/\partial\lambda\)
(quantified below) means even resonant nodes need down-weighting toward the
more open, more informative ones. The nodes are then re-tracked with the
just-fitted \(\alpha(R)\) supplied as the per-node seed (not merely the first
node’s value), so the corrected \(\alpha\) also governs the bound branch before
\(\lambda(R)\) is fit against this second track — both smoothing fits run in the
log-amplitude reparametrized coordinates described under Limitations. A
joint polish (_joint_polish) follows: a single least-squares over
\(\lambda\)’s and \(\alpha\)’s free constants together, minimizing the pole residual (the
real part at every node, the Γ-weighted imaginary part at resonant nodes
only) using the analytic c-product Jacobian — pole_sensitivity chained
through the reparametrization, not a finite difference — which is what
takes F₂ from a several-percent curve error after smoothing alone to
round-off (see the measured table below). Its per-node pole searches recentre
on the same node’s pole from the previous residual evaluation and, once that
pole exists, accept a candidate only within one step (MAX_STEP) of it — the
same continuity guard solve_pole_params and walk_t1 carry, and for the
same reason: find_resonance_pole is a global residual-argmin over its
window with no preference for the state being tracked, so an unguarded wide
fallback window can substitute a different angle-stable state and then anchor
every later iteration on it. A node with no admissible candidate is scored as
a miss (the penalty above), not latched. An electron-affinity (EA)
asymptotic constraint is applied if the node range reaches far enough
toward R → ∞ to bracket it. Finally the fitted model is re-verified by a
fresh gated walk (walk_t1); "met" requires the resulting E_res/\(\Gamma\)
residuals to clear tol.e_res_rms/tol.gamma_rel and at least 75% of
the requested nodes to have survived tracking, re-tracking, and re-walking
— a curve that reads accurate only on a handful of surviving nodes does not
count as fit.
ECS analyticity and the boundedness check¶
Every term in the ansatz is entire in r and in R (aside from y_p’s
poles at |R| = R_e, which an ECS tail pivoted at R_0 > R_e never
reaches), so each has an analytic bound on how far it can be rotated
before its ECS tail stops decaying. Under \(r \to r\, e^{i\theta}\) the Gaussian well’s
exponent becomes \(-\alpha r^2 e^{2i\theta}\), which decays iff \(\cos 2\theta > 0\), i.e.
\(\theta < 45^\circ\); the EMO neutral curve’s exponentials pick up only \(e^{i\theta}\) and
stay bounded iff \(\cos\theta > 0\), i.e. \(\theta < 90^\circ\). These caps
(electronic_max_deg = 45, nuclear_max_deg = 90) hold for every member of
the ansatz family, independent of any grid — they are a property of the
functional form, not a measurement.
report.ecs_bounded(model, pair, R_tail) is the complementary, measured
check: it evaluates |V| on the two tails a given fit actually used —
pair.grid_a’s ECS-rotated electronic tail and the nuclear probe tail
R_tail — and compares each to |V| over the real region, raising
ValueError if either tail’s potential has grown by more than 10× (not
absorbing). The angles actually probed (probed_electronic_deg,
probed_nuclear_deg) are reported alongside the analytic caps so a caller
can see the margin, not just a pass/fail; tail_growth is the one number in
the returned dict that is a genuine measurement rather than a derived
constant. FitReport.ecs_bounds_deg carries all four fields (the two
analytic caps and the two probed angles) plus tail_growth, so a saved
report records both what the ansatz can tolerate and what a specific fit
used.
Round-trip oracle — measured¶
The acceptance test (test_roundtrip.py) round-trips each of N₂, NO, F₂:
seed a perturbed FlexibleDiatomicModel (D_e ×1.2, R_e ×1.05, \(\beta_0 \times 0.9\),
\(\lambda\).f_inf \(\times 1.1\), \(\alpha\).f_inf \(\times 1.3\) from the published model), fit it against
extract_target’s own tiered Target built from the unperturbed published
model, and check the recovered curves and constants. All three currently pass
all three tiers (“met”) through fit()’s T0 → T1 → T3 pipeline:
N₂ |
NO |
F₂ |
|
|---|---|---|---|
T0 (neutral curve) |
met |
met |
met |
T1 (resonance curve) |
met |
met |
met |
T3 (coupling / width) |
met |
met |
met |
T0 recovers D_e, R_e, \(\beta_0\) to round-off on every molecule (rms
~7e-17, max ~2e-16 — the least-squares problem is well-posed independent
of the T1 seed).
T1, per molecule:
N₂ (9/9 target nodes, 5 resonant): \(\lambda(R)\)/\(\alpha(R)\) curve relative error
1.8e-13/1.2e-13;E_resrms2.10e-14Ha, max3.49e-14Ha; \(\Gamma\) relative rms and max0.000; the fitted crossing radius lands within 0.1 bohr of the publishedR_c = 2.405.NO (10/10, 4 resonant): curve relative error
4.3e-10/3.6e-10;E_resrms1.59e-11Ha, max3.58e-11Ha; \(\Gamma\) relative rms and max0.000; fittedR_cwithin 0.1 bohr of the publishedR_c = 2.285.F₂ (8/8, 3 resonant): curve relative error
1.2e-9/8.8e-10;E_resrms1.01e-11Ha, max2.24e-11Ha; \(\Gamma\) relative rms and max0.000. The recoveredSmoothRconstants match the published ones to 8–9 significant figures:constant
recovered
published
\(\lambda\).
f_inf18.848999984001
18.849
\(\lambda\).
f_0−8.874597817145
−8.874597806914
\(\lambda\).
f_13.212999999418
3.213
\(\lambda\).
R_f1.831999999437
1.832
\(\alpha\).
f_inf2.999999997362
3.0
F₂’s clean result depends on one condition its window was chosen to satisfy,
explained below under Limitations: enough of the surviving nodes must be
genuinely resonant (not bound, where \(\alpha\) is frozen) to break a near-degeneracy
between \(\lambda\) and \(\alpha\). It no longer depends on the fit’s R grid happening to
coincide with extract_target’s survivors — T1 now reads the target on
those nodes.
T0 also reports the v=0 → v=1 vibrational spacing of the fitted curve
against the target’s own, and gates "met" on omega_e_rel. Measured on the
round trip: N₂ 9.9e-4, NO 2.1e-3, F₂ 1.5e-3 — all bounded by the
extraction table’s own cubic-spline interpolation error (10 nodes over the
window), not by the fit, which reproduces the curve to ~1e-16. With the
curve replaced by exact constants (NeutralTarget.constants) the same check
reads 1.5e-13.
T3 is met on all three by the same fit_coupling least-squares over
\(\log\tilde\Gamma(\varepsilon,R)\) described below.
Wall time for one full round trip (T0+T1+T3, ElectronicPair()’s default
grid): N₂ ≈ 4.1 s, NO ≈ 4.5 s, F₂ ≈ 3.7 s. The full
projects/potential_factory -m slow suite — the three-molecule round trip
plus the fit_neutral/fit_resonance/fit_coupling slow tests — measures
31.6 s on the same default grid.
Base experiments — converged curves and the observable-level round trip¶
validation/factory/base_experiments.py (python -m validation.factory.base_experiments --molecule N2|NO|F2 --stage curves|fit|xs|all) runs the factory on each published model at production
quality and asks the question the toy-grid round trip cannot: what does
curve-level agreement buy at the level of the exact 2-D cross sections?
Nothing here is compared with experiment. The runs below were made on
2026-08-25 in the MUMPS container (single-threaded BLAS — on tiny 113×113
electronic eigenproblems a 32-thread OpenBLAS is ~400× slower than one
thread, so the ladder below takes seconds, not minutes); the CSVs,
convergence summaries and FitReports are committed under
validation/factory/results/.
Converged resonant curves. walk_t1 extracts E_res(R), \(\Gamma(R)\) and
V_ion = V_0 + E_res on 45 nuclear nodes over the whole crossing region for a
ladder of four electronic grids (r_max 16→40 bohr, DVR order 7→13, 6→10 tail
elements; 113→347 points). The grid-to-grid differences on the finest step
(r32_o11_c8 → r40_o13_c10) are the convergence statement:
N₂ |
NO |
F₂ |
|
|---|---|---|---|
scan |
3.2 → 1.5 |
3.4 → 1.6 |
4.2 → 1.9 |
max \(\Delta E_\mathrm{res}\), finest step (Ha) |
5.5e-9 |
3.6e-7 |
6.2e-7 |
max \(\Delta\Gamma\), finest step (Ha) |
4.1e-9 |
1.2e-6 |
2.9e-5 |
nodes with a gated pole |
44/45 |
44/45 |
44/45 |
\(\Gamma\) at the innermost node (eV) |
4.5 |
— |
14.1 |
The one node without a gated pole is, for every molecule, the one inside the
crossing slice where the pole sits below the gate’s threshold floor. Away from
that slice the curves converge to 1e-12–1e-14 Ha on the bound side; the
residual differences on the resonant side come from the very short-R end,
where the width is several eV and the real-region extent still matters
(F₂’s \(\Gamma\) reaches 14 eV at 1.9 bohr). The crossing found by the factory agrees
with the published R_c to 3e-4–2e-3 bohr (N₂ 2.4048 vs 2.405, NO 2.2845 vs
2.285, F₂ 2.5970 vs 2.595). Figures:
docs/physics/figures/{n2,no,f2}-factory-curves.png.
Refit at production quality. The factory refits each model from its own
converged curves (24 nodes, n_eps = 8, the finest grid) from the perturbed
seed of the round-trip test. All three meet all three tiers; the T1 E_res
rms is 1.7e-12 Ha (N₂), 2.3e-12 Ha (NO) and 6.7e-7 Ha (F₂ — the largest,
from the two crossing-slice nodes the tracker drops and re-verifies by
interpolation), the T3 log-width rms 1e-9–7e-9. The DA-threshold signs come
out + for N₂ and NO (endothermic) and − for F₂ (exothermic).
Cross sections, published vs refitted. Exact 2-D \(\sigma_{0\to v'}(E)\) for
v' = 0, 1, 2 on each molecule’s emoscat production deck over its preset
window, and \(\sigma_\mathrm{DA}\) where the channel exists — F₂ over the same window, NO on a
window above its DA threshold (E = 0.175–0.30 Ha; the channel opens at
+0.172 Ha, above the whole VE sweep):
N₂ |
NO |
F₂ |
|
|---|---|---|---|
deck (unknowns) |
26,857 |
78,804 |
128,568 |
VE energies (sweep, s) |
40 (2) |
30 (10) |
25 (15) |
max |σ_refit/σ_pub − 1|, VE |
1.9e-10 / 3.3e-10 / 5.2e-10 |
5.2e-10 / 1.5e-10 / 1.4e-9 |
9.0e-9 / 8.9e-9 / 8.7e-9 |
max |σ_refit/σ_pub − 1|, DA |
— |
4.3e-7 (σ ≤ 7.4e-10 bohr², 4.7–8.2 eV) |
4.3e-9 (σ up to 5.6 bohr²) |
So a model rebuilt from curves that agree to 1e-12 Ha reproduces the exact
observables to 1e-9 — including N₂’s boomerang structure point by point and
F₂’s DA over three orders of magnitude — and the 1e-7 on NO’s DA is the
exponentially small (1e-10–1e-19 bohr²) tail of a closed-at-threshold channel
where the cross section itself is at the solver’s floor. Curve-level agreement
at this level is sufficient; how much curve-level disagreement the
observables tolerate is the sensitivity budget, still to be measured. Figures:
docs/physics/figures/{n2,no,f2}-factory-ve.png,
{no,f2}-factory-da.png.
O₂ — the first real-molecule target, from the figure to its VE cross section¶
validation/factory/fit_o2.py fits the factory to O₂ as published by Alt &
Houfek, Phys. Rev. A 103, 032829 (2021) (reference/literature/ alt-houfek-2021-pra103-032829.md). The target is not digitised by eye: the
paper’s Fig. 2 is embedded as vector graphics, and validation/factory/ extract_fig2.py recovers V_0(R), the anion curve V_ion(R) (the dashed real
part of the resonance energy below the crossing, the bound curve above it)
and \(\Gamma(R)\) from the PDF’s paths — axes calibrated from the tick marks
(residuals 3e-4 bohr, 1.4e-3 eV), each curve’s centreline as the median of
its filled outline’s two edges per bin — to a vertical precision of ~0.02 eV
(validation/factory/data/o2/README.md). The extracted V_0 has its minimum
at −5.259 eV, Table I’s calculated D₀ plus the zero-point energy, and its
ladder G(1)−G(0) = 1607 cm⁻¹ against the spectroscopic 1551 (the paper’s
MRCI curve is ~3 % stiffer than experiment; the T0 check therefore uses the
curve’s own ladder). The energy-dependent width comes from Table II’s
\(\tilde\Gamma(\varepsilon,R) = 2\pi\,\varepsilon^{5/2} A(R)\,e^{-B(R)\varepsilon}\); the electron-affinity asymptote from
Table I’s EA(O) = 1.4611 eV. Because the data are a figure, the tolerances of
this phase are the extraction floor (IMAGE_TOL: 20 meV on V_0, 40 meV on
E_res, 20 % on \(\Gamma\)), not the sensitivity budget; fit() runs with
continue_on_miss so every tier reports.
The asymptote is theory, not the figure. Both curves end at energies
that are atomic constants, and the fitter treats them as such rather than
reading them off a table that stops at 6 bohr: V_0(R→∞) = 0 (free atoms; the
EMO form has it built in) and V_ion(R→∞) = −EA(O) (O + O⁻; Table I),
approached through the ion–atom polarisation \(-\alpha_d/(2R^4)\) with \(\alpha_d(\mathrm{O}) =
5.3(2)\) a.u. (reference/literature/schwerdtfeger-nagle-2019-molphys117-1200.md).
Two things in that sentence are the operator’s per-molecule judgement,
declared on the ResonanceTarget (R_inf, tail) and not library
constants: which tail applies, and from where the curve follows it. For
O₂, Fig. 2 still shows 0.15 eV of binding at 6 bohr beyond the polarisation
term — the O⁻ ↔ O charge resonance, decaying as exp(−√(2 EA) R) with
√(2 EA) = 0.33/bohr — so it is under the 20 meV extraction floor from ~12
bohr and R_inf = 14 is the choice recorded in targets/o2.py, with its
reasoning. Nothing is placed between the table’s end and R_inf: there the
curve is neither tabulated nor asymptotic, and the model interpolates. The
polish pins the anion’s bound state to −EA + tail(R) at R_inf, 2.2 R_inf
and 5 R_inf, and the report carries a fourth tier, asymptote, that
re-measures those residuals and V_0(R_inf) on the final model. The far node
exists because one node at 10 bohr proved worthless: the sigmoid’s inflection
ran off to large R, held −EA at 10 and sat 0.2 eV above it at 20 (measured
on the first O₂ fit).
The \(\lambda(R)\) form. With the asymptote pinned, the sigmoid × polynomial
SmoothR could hold the table or the tail but not both (pinning it cost
T1 20 meV and 30 % of \(\Gamma\) and still missed by 3.5 mHa), and over the full
1.85–6.0 bohr range it fell into a wrong basin (\(\lambda\) flat at 3.55, \(\alpha\) 0.245)
regardless. TailR replaces it for O₂: \(\lambda(R) = \lambda_\infty + (1 - y_q(R))\,P(y_p(R))\)
with q = 4, so every term dies as R^{−4} — the polarisation power — and
\(\lambda(\infty) = \lambda_\infty\) exactly; P is a plain polynomial in Le Roy’s bounded y_p, so
the fit is linear in its coefficients and has no inflection to run away. The
published models keep their sigmoid (from_diatomic embeds them exactly),
and a floor f_1 ≥ 1/(R_max − R_min) now stops a sigmoid from degenerating
into a constant under the polynomial (the O₂ fits had reached \(f_1 = 10^{-6}\),
R_f = 25, coefficients of order 10²).
Result (2026-08-25, range 1.85–6.0 bohr, 40 nodes, 5-term EMO, TailR
\(\lambda(R)\) with nine and \(\alpha(R)\) with three polynomial terms, r32/order-11
electronic pair, R_inf = 14, 87 s on the laptop;
validation/factory/results/o2-fit-report.json,
docs/physics/figures/o2-factory-fit.png):
tier |
status |
measured |
|---|---|---|
T0 neutral curve |
met |
rms 14 meV, max 86 meV (the wall); ladder to 0.3 % of the curve’s own |
T1 pole curves |
met |
|
asymptote |
met |
\(V_\mathrm{ion} - (-\mathrm{EA} - \alpha_d/2R^4)\) = −0.08 / −0.02 / −0.02 mHa at 14 / 31 / 70 bohr; |
T3 energy-dependent width |
not met |
log-rms 1.1 against Table II; the shell fit cannot move it (see below) |
Beyond the table the model is monotone: V_ion + EA = −89, −44, −15, −4,
−1.4 meV at 7, 8, 10, 14, 20 bohr, the charge-resonance-plus-polarisation
approach one expects. The fitted \(\lambda(R)\) runs from ~4.3 at 1.85 bohr through a
maximum of 5.23 near 3.5 to \(\lambda_\infty = 5.06\) with \(\alpha \approx 0.36\) — a d-wave well of
N₂’s kind, as expected for the same \(^2\Pi_g\) symmetry. The whole run is 13×
faster than the sigmoid fits (85–87 s vs 430–1164 s) because the polish now
converges in under twenty evaluations.
The spectral check — the metric that predicts the VE figure. O₂’s
vibrational-excitation cross section is a comb of narrow peaks at the anion’s
quasi-bound levels, so the number that matters is not an rms over the curve
but the levels: validation/factory/o2_levels.py solves the 1-D nuclear
problem in \(V_\mathrm{ion}(R) - i\Gamma(R)/2\) once in the extracted curves and once in the
fitted model’s (a seeded walk_t1), and compares E_v − E_0(neutral) — where
each peak sits — and \(\Gamma_v\) — how wide it is
(validation/factory/results/o2-anion-levels.csv). On the committed model
the peak positions agree to ±7 meV over v = 0…24 (0–2.3 eV, the paper’s
whole Fig. 5 window) and the widths to ~10 % from v ≈ 7 up; the sub-0.1 meV
widths of v = 2–5 are 30–50 % over, set by \(\Gamma(R)\) right at the crossing where
the extracted curve is at its precision floor. Two earlier fits make the point
that this check is not redundant with T1: the 4.5-bohr sigmoid fit (T1 rms 74
meV) put every level −20 meV low with widths 10–25 % broad, and the first
full-range TailR fit with five coefficients (T1 rms 72 meV — the same rms)
scattered them from +73 meV at v = 0 to −3 at v = 8 and +62 at v = 24, with
widths 30–50 % broad, because its error sat in the anion well rather than on
the wall. The rms did not distinguish them; the levels did.
What O₂ forced on the fitter (all now in fit.py/extract.py): the shell
is bounded repulsive — free, it dug a second well at large R that bound a
state of its own and the verification walk followed that state instead; every
verification walk is seeded on the target’s own pole, because the default
window’s global residual-argmin picks the deepest angle-stable state in a
well that holds more than one; the electron-affinity asymptote is a set of
nodes of the polish (bound poles at R_inf and beyond, see above) rather
than a shift of lam.f_inf applied afterwards, which moved \(\lambda(R)\) at every
R and undid a polish that had reached 1.3 mHa; T3 re-walks the T1 nodes after the shell is
fitted and fails if the shell broke them; a crossing-slice node (E_res > 0
with \(\Gamma\) under the floor) is “no target” rather than an impossible bound
state; and the EMO’s \(\beta(R)\) is kept positive by bounds and a penalty, since
a 5-term EMO free to go negative at small R produced a curve binding no
levels at all.
One open finding. The T3 comparison is not discrete-state-consistent: the factory’s \(\tilde\Gamma(\varepsilon,R)\) is built on the R-independent discrete state (choice B), while Table II is a Breit–Wigner fit of R-matrix eigenphase sums, i.e. the width of the resonance pole as a function of energy — their energy dependences differ (the factory’s is too flat: ~ε^{1.5} over 0.05–0.5 eV where Table II goes as ε^{2.5}e^{−Bε}), and the repulsive shell cannot reconcile them. A real-molecule T3 target should be compared through the physical, R-dependent discrete state (choice A) or through the pole width \(\Gamma(R) = \tilde\Gamma(E_\mathrm{res}(R), R)\) alone, which T1 already gates. It belongs to the O₂ plan’s next phase, together with the authors’ tables; it does not stand between this model and its VE cross section, which the exact 2-D solver computes from the potential surface alone.
From the fit to the run surface. The ansatz is promoted to
qscat.model.flexible (FlexibleDiatomicModel, SmoothR, TailR,
from_diatomic; projects/potential_factory/ansatz.py re-exports it), and
the fitted model is the registry entry qscat.model.O2 — the first entry that
is a fit rather than a published parameter set, its constants copied from the
committed report key for key and locked to it by
validation/factory/test_o2_report.py. That lock caught a real gap: the
y_p reference radius inside \(\lambda(R)\)/\(\alpha(R)\) is a frame constant the report
does not record (it stays at the seed’s 2.2819 bohr, not the fitted R_e of
2.268), and a registry built from the report alone was 1e-2 Ha off in the
well. Its decks come from the discretisation tuner (validation/factory/ o2_grids.py): the electronic grid as proposed (324 points, real region to
20 bohr, 35° tail to 138); the nuclear grid’s mesh as proposed but its real
region cut at 8 bohr with the tuner’s tail re-attached — the VE path’s fixed
18-bohr extent put 62 elements where nothing lives (DA is closed until 3.7 eV;
the anion’s outer turning point at 2.3 eV is 4.0 bohr) — 289 points. All
three 1-D probes pass at the window’s top (0.10 Ha): channel representation
4e-5, the anion bound state to 2e-10, the 12 lowest neutral levels to 1.2e-3
(the last level marginal). The 2-D spot check overruled them, exactly as
it did for F₂: one nuclear h-refinement moved \(\sigma_{0\to 1}\) at 1.36 eV by
69 %, after which both further refinements move it under 2 %. The reason
is the observable, not the mesh: the cross section at an energy sits on a
comb of 1–8 meV peaks, so a level shift of 3 meV — the 1e-3 the probe
tolerates — is most of a width. The carried nuclear deck is therefore the
tuner’s mesh refined once (549 points, 324 × 549 = 178k unknowns — MUMPS
territory, 440 s for the whole refine loop on the laptop with SuperLU). Two
small things the tuner needed for a fitted model: max_stable_angle read round-off as growth on
O₂’s fast-decaying EMO (|V_0| is 1e-17 at the pivot; β(∞) = 3.9 against
N₂’s 1.15) and returned 0°, so growth below an absolute 1e-12 Ha floor no
longer counts. The apps/qscat-run preset O2:tuner carries the decks
verbatim (test_o2_grids.py locks them), VE only, and the example
o2-ve.yaml is generated by validation/factory/o2_ve_energies.py with a
level-aware mesh — a 0.5 mHa background grid plus 121 points across ±6 widths
of every anion level from the spectral check’s own table (3343 energies) —
because O₂’s peaks are 0.01–8 meV wide and any uniform sweep one can afford
walks past them.
The comparison is theory against theory. The paper’s Fig. 5 (p. 032829-6)
— its own nonlocal-model and LCP cross sections for 0 → v’, v’ = 0…5 — is
vector graphics like Fig. 2, and validation/factory/extract_fig5.py pulls
both curves of every panel out of it (validation/factory/data/o2/ fig5_ve_0{v}_{nrm,lcp}.csv). Two things differ from Fig. 2: the tick labels
are glyph outlines rather than text, so each panel’s ranges are read off the
rendered page once and fixed in the extractor while the tick rectangles do
the calibration; and the centreline is the outline’s upper envelope per
0.15-pt bin minus half the stroke width (0.25 pt, 0.3–0.7 % of a panel’s
range; energy resolution 0.9–1.4 meV), because a comb of meV-wide peaks is
exactly where the median-of-edges centreline fails — an x-bin on a peak holds
both walls. One trap is worth recording: gnuplot strokes the legend key
samples inside the curve paths, so unmasked they appear as tall spikes at
the top right of every panel, and every panel’s “maximum” but the first was
the key line at 0.94 of the range; they are masked by their measured position
(a real peak crossing the two thin bands keeps its top to 0.0002 of the
range). Fig. 5’s curves are spin–orbit resolved — every peak is a doublet,
the sum of a ²Π_{1/2} and a ²Π_{3/2} component of statistical weight 1/3
each — so the only comparison that means anything is with the model treated
the same way; validation/factory/o2_ve_figure.py overlays the sum of the
two components’ exact 2-D cross sections (next section) on those curves. The
overlay answers two separate questions: whether the exact 2-D comb of the
fitted potential lands where the paper’s nonlocal comb does (positions are
the fit’s, to the spectral check’s ±7 meV; heights are where the
discrete-state choice shows), and whether the paper’s LCP fails the way this
repository’s own LCP does. The paper’s measured traces (its Figs. 7–9) are
not extracted; nothing is compared with experiment.
The energy mesh has to draw the peaks. The sweep on the 178k-unknown deck runs at 0.38 s per energy in the MUMPS container on a 32-core host (the committed 3343-energy, six-channel sweeps in ~1290 s each); the same 80-energy slice took 61 min on the laptop with SuperLU, 46 s per energy — what the sparse-LU hot path costs without MUMPS. Cheap energies are what make the mesh affordable, and the mesh decides the heights: at 15 points across ±5Γ per level (spacing Γ/1.5) every peak height read 0.6–0.9× its converged value — a Γ/3 miss on a Lorentzian is exactly ×0.69 — at 41 points (Γ/4) the heights came within ~10 %, and at the committed 121 points across ±6Γ (Γ/10) the peaks are drawn in full; a 0.02-meV-wide level (the elastic peak at 0.199 eV) is resolved, a 0.004-meV one (0.069 eV) is not.
The spin–orbit splitting. The paper’s Fig. 5 doublets come from the
²Π_g curve being two curves: the ab-initio splitting Δ_SO(R) (their Fig. 1 —
19.7 meV at 2.1 bohr, a 12.2 meV atomic plateau beyond 9 bohr; a Gaussian
fit the paper used directly) puts ²Π_{1/2} and ²Π_{3/2} symmetrically at
∓Δ_SO(R)/2, each is run independently with the same width, and the cross
sections are summed with the statistical factor 1/3 each (§III A, p.
032829-4). The factory does the same thing to the potential: Fig. 1 is
vector-extracted (extract_fig1.py, ~0.05 meV), o2_target(so=±1) moves the
anion curve and its atom + ion asymptote by ∓Δ_SO/2, and fit_o2_so.py
produces qscat.model.O2_SO12/O2_SO32 from O2 by polishing only
(refine_resonance: every node’s pole seeded from the model’s own gated
walk, the asymptote nodes attached, the joint polish moved from where the
coefficients are). That last point was measured, not assumed: re-running the
full T1 pipeline — tracking from a seed and re-fitting the smooth forms —
against a curve shifted by 10 meV landed ²Π_{1/2} in a wrong basin (E_res
rms 31 mHa, asymptote 28 mHa off), while the polish reproduces the parent’s
quality on both components to the digit (T1 rms 0.739 mHa, asymptote 0.05
mHa, 77 s and 50 s; the spectral check ±7 meV on each against its own
shifted curve). The components share O2’s deck (their curves differ by far
less than the discretisation resolves) and get their own level-aware meshes;
o2_ve_figure.py --so12 --so32 sums the two sweeps at ⅓ each and overlays
Fig. 5 directly, doublets included, and reports the separation of the 0 → 1
doublet near 1.04 eV — Allan’s v’ = 9 pair, which the paper’s model puts at
17.8 meV against a measured 19.6 ± 1.0.
Result (2026-08-26; validation/factory/results/o2-so{12,32}-ve/,
docs/physics/figures/o2-2d-ti-ve-spin-orbit-vs-alt-houfek.png). Two
sweeps of 3343 energies each (six channels, Γ/10 mesh), ~1290 s apiece with
MUMPS. Summed at ⅓ each, the exact 2-D doublets sit on Fig. 5’s in all six
panels: both members of every doublet in 0 → 1…3 within 1–8 meV of the
paper’s (0.3166/0.4604/0.5821/0.6822/0.7018 eV against 0.3186/0.4666/
0.5888/0.6858/0.7039 in 0 → 1 — the spectral check’s ±7 meV made visible),
with heights 0.9–1.1 of the paper’s on the main peaks; the weak 0 → 5 channel
(peaks of 0.005–0.017 a₀²) reads 0.8–0.96 on its stronger members; and the
paper’s LCP is offset from both by 10–40 meV in the inelastic channels,
exactly as this repository’s own LCP is on N₂ and F₂. The weakest peak of
each inelastic channel, at its threshold, is where the fitted model runs
high (1.4–1.6× at 0.57 eV in 0 → 2, 2× at 0.81 eV in 0 → 3) — the honest
exact-vs-nonlocal and fitted-vs-published residual to chase with the
authors’ tables; everything else is the paper’s comb, from a potential
fitted to a figure. The doublet separation is 19.0–19.3 meV
across the 0 → 1 comb (19.0 meV at the v’ = 9 pair near 1.04 eV), against
the paper’s 17.8 and Allan’s measured 19.6 ± 1.0 — our components are
separated by Δ_SO(R) at the levels’ own ⟨R⟩ ≈ 2.3–2.6 bohr, where Fig. 1
reads 18–19.5 meV. That the number lands
nearer the measurement than the paper’s own is noted, not claimed: it is a
theory-vs-theory comparison of two nuclear-dynamics treatments of one
splitting curve, and nothing here was tuned to it.
Limitations¶
The sigmoid-constant degeneracy. N₂’s own published \(\lambda(R)\) sigmoid has
extreme raw constants (\(R_f = -27.98\), \(f_0 \approx -7\times 10^{13}\)) — the sigmoid’s
midpoint sits far outside the physical R range, so the curve is shaped by
the ratio of huge, nearly-cancelling terms. Fitting those constants directly
with an unscaled least squares is catastrophically ill-conditioned: it can
converge to a materially different point (f_0 doubled, R_f shifted ~5%)
that still leaves a multi-percent residual against the very tracked data
being fit. The fit therefore reparametrizes in log-amplitude coordinates
(u = log|f_0|, sign frozen from the seed) with Jacobian-scaled steps before
optimizing, and — the more fundamental point — the fitter is judged, and
"met" is decided, by the curve \(\lambda(R_{\mathrm{desc}})\)/\(\alpha(R_{\mathrm{desc}})\) the constants
reproduce, never by how close the raw constants land to the published ones.
Two sigmoid parameterizations can differ substantially in every constant and
still trace the same curve to round-off.
The \(\lambda\)/\(\alpha\) near-degeneracy. On the single-Gaussian well, the two Hellmann-Feynman gradients \(\partial E/\partial\lambda\) and \(\partial E/\partial\alpha\) are nearly parallel: measured on F₂ at \((\lambda, \alpha) = (17.46, 3.0)\), \(\partial E/\partial\alpha \approx -6.06 \cdot \partial E/\partial\lambda\). A single pole therefore barely constrains \(\lambda\) and \(\alpha\) independently — the per-node 2×2 Newton Jacobian is close to singular along that direction — and a bound node cannot help at all (\(\alpha\) is frozen there, contributing no equation in \(\alpha\)). Enough genuinely resonant nodes are needed to break the degeneracy from several different points along the near-null direction; two were not enough for F₂, three were.
T1 evaluates the target on the target’s own nodes. walk_t1 gates out
any node it cannot pin a pole at and reports only the survivors, and the
target curve is a spline through those survivor nodes — exact there, but
interpolated elsewhere. E_res(R) has a branch point at the resonance
crossing, and no polynomial interpolant is accurate across it: reading the
target even one node off the survivor set cost, measured on F₂, 1.39e-3 Ha
and 1.67e-3 Ha of spurious target error at the two fit nodes straddling the
threshold — 140–170 % of the e_res_rms = 1e-3 tolerance, before the fit
itself had done anything. fit_resonance therefore builds its R grid from
target.v_ion.x (restricted to R_range, evenly subsampled to at most
n_nodes, descending) whenever the target is table-backed, falling back to a
linspace only for a callable-backed target that has no nodes and no
interpolation error to avoid; the same grid is used for the re-verification
walk, and TierResult.detail reports which was used (nodes: table 8/8 /
nodes: linspace 24). This is the convention T0 (fit_neutral, at
target.curve.x) and T3 (_coupling_eval_grid, subsampling
eps_nodes/R_nodes) already followed. Measured effect on the F₂ round trip
at the (3.6, 1.8) window, which drops one extraction node: with an
independent linspace grid T1 was not met (E_res rms 4.6e-3 Ha, \(\Gamma\)
relative max 0.51, \(\lambda\) 65 % off, T3 never attempted, 12.1 s); on the
target’s own nodes it is met at E_res rms 1.01e-11 Ha with \(\lambda\)/\(\alpha\)
recovered to ~1e-9, in 3.8 s.
Tolerances (report.py) are placeholders — v0_rms, omega_e_rel,
e_res_rms, gamma_rel, gamma_floor, coupling_log_rms are round numbers
chosen to make the round-trip oracle discriminating, not the output of a
measured sensitivity-to-observable budget. A future budget pass should derive
each from how much a given tier residual actually moves a downstream cross
section.
The fake near-threshold pole and the gate. A bare Gaussian well has
no ℓ = 0 bound or resonant state at all — an s-wave target needs the
optional shell (a barrier) to have anything to fit, and without one
find_resonance_pole can still report an angle-stable “pole” near threshold
that is a numerical artifact rather than a physical state (the known weakness
this package’s residual/gamma gate exists to reject — see
docs/physics/h2plus-resonance-states.md). The gate catches the specific
spike case measured during development (residual 2.7× over its own
tolerance), but is not a proof that every spurious near-threshold match is
caught.
Polar targets are out of scope, and T3 says so. A polar molecule’s
threshold exponent \(\varepsilon^{\sqrt{d+\frac14}}\) (a dipole-coupling law) cannot be produced by
a centrifugal barrier of any integer ℓ; matching it needs an explicit
long-range term the current ansatz does not have. fit_coupling therefore
compares CouplingTarget.alpha_exponent against the ansatz’s own ℓ + ½
first and, on a mismatch, returns "not met" naming both exponents without
fitting — rather than fitting a shell against a threshold law the model cannot
produce and reporting whatever residual that leaves.
Untested code paths. SmoothR with coeffs != () (the polynomial
correction beyond a pure sigmoid) and FlexibleDiatomicModel.params()’s
shell-present branch are implemented but not exercised by any test on N₂/NO/F₂
(none of the three published models needs either). fit_resonance’s EA
(electron-affinity) constraint has a bracket-search fallback for when the
naive bracket fails; it is likewise implemented but untested, since none of
the three round-trip windows reaches the R → ∞ asymptote where it would
fire. There is no guard against a Target table with fewer nodes than the
ansatz has free parameters at a given tier — an under-determined fit would
currently proceed rather than raise. A node whose target is NaN (dropped by
track_curve) carries forward its raw, un-refit \((\lambda, \alpha)[j]\) rather than an
interpolated value; this is a deliberate, documented choice (continuation
needs some seed to resume from) but means a caller inspecting TrackResult
directly should not read a frozen node’s parameters as a fit result.
fit.py has grown to ~1000 lines holding three fitting stages (T0/T1/T3)
plus the smoothing reparametrization and the joint-polish refinement in one
module. The natural next seam before promoting any of this out of
projects/ is splitting the smoothing/polish machinery into its own module
(e.g. smooth.py), leaving fit.py as orchestration.
Wall time. The full projects/potential_factory -m slow suite measures
34.2–34.6 s on ElectronicPair()’s small default grid (r_max=16, order=7, n_complex=6); a production-sized grid, or a molecule needing a larger ℓ or
a finer R continuation, will cost more per pole search and per Newton
step, and the factory has not yet been run at production grid sizes.