Test whether the SCF reference is the ground state, and act on the answer

[calvinp0]

⚠️ Changes a physical quantity. Reviewed by a quantum-chemistry pass as well as an adversarial code pass, three times each.

Teaches ARC to test whether the SCF reference it chose is actually the ground state, and to act on the answer. Gaussian and ORCA.

Off by default ('stability': False in settings.py) — a run that does not ask for it is unchanged.

What it adds

An opt-in stability job type, spawned once per species when its optimisation converges — for every TS, and for any other species whose opt actually ran restricted, which is the only reference the analysis can inform. The analysis is a single point, so it needs no frequencies.

The optimisation's tail — freq, sp, IRC, rotors and the rest — is held until the verdict is in, because each of those inherits the optimisation's reference or its geometry:

opt → stability → freq / sp / IRC / rotors
opt → stability → opt (unrestricted) → freq / sp / IRC / rotors

An instability means the geometry is wrong too: it is a stationary point of the restricted surface only. So an adopted verdict re-optimises rather than merely re-referencing — a Hessian taken at the restricted geometry on an unrestricted reference sits at a non-stationary point and can produce spurious imaginary modes. At most one re-optimisation per species, guarded by stability_reoptimized on ARCSpecies and carried in restart.yml, and a resumed run releases any species whose analysis finished while it was down.

The re-optimisation's starting orbitals are per-ESS, because the two codes leave different things behind. ORCA follows the instability and writes the relaxed broken-symmetry orbitals, which seed the re-optimisation. Gaussian's stable=(rext,noopt) does not follow, so its checkfile still holds the restricted orbitals — and a restricted determinant is a stationary point of the unrestricted equations, so reading it back converges to the solution the analysis rejected. That case drops the checkfile instead and takes Gaussian's guess=mix.

• Gaussian: stable=(rext,noopt) at the freq level on the freq geometry.
• ORCA: a single point with STABPerform true / STABRestartUHFifUnstable true, plus .gbw orbital tracking.

A parser reads the verdict, the negative stability-matrix roots and their eigenvalues, and derives whether the analytic frequencies are invalidated. An unreadable verdict reports unknown, never stable.

Alongside it, the <S**2> spin-contamination diagnostic is re-homed from an arcbench-based branch onto main (Gaussian, ORCA, Q-Chem), and arc/checks/spin.py gains the Yamaguchi approximate spin-projection arithmetic.

The adoption contract

1. A user-declared number_of_radicals always wins and is never overwritten by a calculation.
2. ARC still runs and still assesses the check regardless.
3. Disagreement is a warning, never a crash; both pictures are recorded.
4. Only when nothing is declared does an external (R→U) instability of a restricted reference get adopted for subsequent jobs, with provenance recorded.

An internal instability never flips the reference — it is a lower solution inside the reference's own spin symmetry, which is a different problem and not evidence of broken-symmetry character.

The derived verdict lives in its own ARCSpecies.derived_stability_verdict, deliberately not in number_of_radicals, which feeds molecular-graph perception in eight places plus xTB's UHF count. A measured SCF property must not decide which molecule ARC thinks it has.

Measurement widened to every restricted species; adoption stays TS-only. A well is never re-optimised on its verdict, so adopting for one would guarantee the E_elect(unrestricted) + ZPE(restricted) splice a TS only risks. Because adoption is TS-only, a TS whose restricted reference was unstable runs unrestricted (broken-symmetry, contaminated, biased low) while its reactants and products stay restricted — so the barrier is systematically underestimated, and nothing projects. That direction is now stated in the docstring and in advanced.rst rather than left implicit.

Provenance — source, declared count, verdict, per-job-type SCF references, mismatch flag and the log it was read from — is recorded in output.yml and survives a restart. A TS switch carries an adopted external instability forward, dropping everything geometry-specific.

The ORCA path, and what had to be measured to write it

ORCA's manual documents the input keywords but not the output, so the parser could not be written from documentation. Thirteen ORCA 6.0.0 jobs were run (five fixtures, four diagnostics, four functional-matched controls) on the four geometries the Gaussian fixtures already use.

Verdict wording, which ORCA's manual does not publish:

The stability analysis shows that the wavefunction is stable
The stability analysis indicates that the wavefunction is unstable

STABRestartUHFifUnstable true is mandatory, not a preference. With it false, ORCA 6.0.0 terminates in LEANSCF with a BLAS incompatible-matrices error (exit 62) after printing a complete verdict — reproduced at nprocs 8 and nprocs 1, at STABNRoots 6 and 3; LeanSCF false fails earlier still. The three stable jobs terminate normally. ORCA therefore cannot do Gaussian's noopt "report but do not follow", and ARC always tells it to follow.

A followed log holds two analysis blocks with opposite verdicts. verdict is taken from the first — the wavefunction under test. The second describes what ORCA relaxed into and is exposed separately as followed_to_stable. Gaussian's noopt logs never exercise this.

On the unstable singlet TS, following the instability gained −0.02022 Eh (−12.7 kcal/mol) and took <S**2> from 0.000000 to 0.864742 — a genuine open-shell singlet whose restricted energy was wrong by 12.7 kcal/mol.

.gbw tracking

ORCA's analysis is an SCF post-step: it converges an SCF first, so without reading the tested orbitals it may converge to a different solution — the hazard Gaussian's checkfile requirement exists to prevent. ARC had no .gbw or MORead handling.

checkfile was already ESS-generic (psi_4 → check.chk, terachem → teracheck.chk), so this rides existing plumbing via check_file_name / guess_file_name class attributes on JobAdapter. ORCA names its own output after the input file, so it cannot read and write one .gbw the way Gaussian reuses one .chk: the guess is uploaded as guess.gbw and read with !MORead / %moinp, while the job's own input.gbw is what returns.

Every ORCA job that runs an SCF on a single structure reads the guess when one exists — opt, conf_opt, optfreq, scan, freq, sp, conf_sp, stability — mirroring Gaussian's guess=read, which sits outside any job-type branch. Job types for which the adapter writes no calculation keyword, monatomic species, and job arrays (whose members share one remote path) read none. The emission set and the upload set are the same predicate, OrcaAdapter.reads_orbital_guess(), and a test asserts they agree across all 15 job types.

This is what keeps the chain consistent. Measured on the C₅H₁₀ TS — identical !UKS B3LYP def2-TZVP input, same geometry, only the guess differing:

	E (Eh)	<S**2>
ORCA's own initial guess	−196.344572	0.000000
!MORead	−196.364789	0.864739

12.7 kcal/mol. The fresh guess collapses onto the closed-shell solution, which the stability analysis reports unstable at that geometry (lowest root −0.0647); reading the optimisation's orbitals holds the stable broken-symmetry solution. A basis change between jobs is handled by ORCA itself — a def2-SVP job reading a def2-TZVP .gbw logs N(Shells)= 6 and 11 - projection required and converges — so no level tracking is required.

The sector is measured, not assumed

ORCA prints one unlabelled stability matrix, which for a restricted reference spans both the internal (singlet) and external (R→U triplet) sectors. Assuming external would drive reference adoption and suppress the analytic-frequency warning on evidence the log does not contain.

Because ARC is forced to run the follow anyway, the post-restart <S**2> is free and decides it: a nominal singlet reaching a stable solution above SPIN_SYMMETRY_BREAKING_S_SQUARED = 0.01 proves the spin symmetry broke. Where the follow never converged, a new verdict unattributed_instability is reported with both flags None — never a fabricated external, and it deliberately does not trigger reference adoption.

Cross-code validation

ORCA analyses RHF/RKS in UHF/UKS space and UHF/UKS in UHF/UKS space. Gaussian's Stable=RExt uses the same Ms-conserving block — all four Gaussian fixtures print Stability analysis using <AA,BB:AA,BB> singles matrix: — so the two tests span the same space and neither reaches the GHF sector.

The two codes agree on the verdict for all four systems. Comparing eigenvalues requires a matched functional, since ORCA's B3LYP uses VWN-5 and Gaussian's uses VWN3; ORCA's matching keyword is B3LYP/G:

system	Gaussian E(0)	ORCA B3LYP/G E(0)	ratio
stable restricted singlet	0.0241461	0.0242519	1.004
stable unrestricted doublet	0.0024619	0.0024662	1.002
RHF→UHF unstable singlet	−0.0642219	−0.0641958	0.9996
spin-contaminated doublet	0.0005803	0.0657159	113×

Total energies agree to 0.0002 Eh, so the functional accounted for the entire raw offset. The outlier is not a code difference: at matched functional its energies still differ by 0.025 Eh with <S**2> 1.7488 vs 1.700055, because the two codes converged to different UHF solutions of a near-dissociated O(³P)···CH₃ pair at r(O–C) = 3.78 Å. No cross-code conclusion is drawn from that system.

What review changed

A shipping blocker. arc/job/adapters/orca.py emitted %moinp "guess.gbw" and uploaded the file, but no ORCA template in arc/settings/submit.py copied it into the scratch working directory — Gaussian has cp "$SubmitDir/check.chk" .; every ORCA template copied only input.in. With repo defaults any ORCA job reading a guess would abort with "Cannot open file guess.gbw", be classified errored/['Unknown'], and be retried forever. It passed development testing only because that cluster's ~/.arc/submit.py overlay happens to glob *.gbw. All six templates now copy guess.gbw in and input.gbw out; all 36 templates were re-checked to still .format().

Two silent-corruption paths. A failed .gbw download leaves a 0-byte file (paramiko opens the local file before the remote), which passed os.path.isfile() and was adopted as species.checkfile; size is now checked. And an ORCA .gbw could be handed to Gaussian as check.chk with guess=read, reachable whenever opt_level routes to ORCA and a later job to Gaussian — adapters now refuse a checkfile that is not theirs, via JobAdapter.readable_checkfile().

A test-ordering bug of the kind that has been destabilising this suite. Five tests asserted record.levelname == 'WARNING', but arc.common.initialize_log calls addLevelName(logging.WARNING, 'Warning: '). They pass alone and fail whenever anything initialises ARC's log first — which -n 6 --dist worksteal can arrange. They compare levelno now.

Earlier rounds: the Gaussian pass refuted the premise the work started from (g16 does not ignore guess=mix on a restricted reference — Mixing orbitals ... Coef= 7.07106781D-01 on a closed-shell singlet, systematic across 24 fixtures). The chemistry pass caught the Yamaguchi arithmetic being silently singlet-only, hard-coding <S²>_LS = 0 — a 16.8 kcal/mol error for a BS-doublet/HS-quartet pair. The adversarial pass measured the separation guard returning −101 Ha from a −1.0 Ha reference; the floor is now physical rather than a divide-by-zero guard, and NaN/inf are rejected. before-annihilation is the <S**2> fed to the projection, because that is the value the SCF energy belongs to — the annihilated one biases it +21.2 kcal/mol.

Discoverability

run_stability_job refuses an unsupported ESS with a warning naming the species and the ESS, emitted once per ESS per run. It is a warning rather than an info because it means the job type the user explicitly switched on will never run for any species that ESS handles — a configuration mismatch, not a situational skip — and once per ESS so a 200-species project gets one line rather than 200.

docs/source/input_reference.rst lists stability among the job-type keys and in the job_types example, and states which keys default true and which false.

Verification

11 commits, no file touched by more than one. Commit order is bisect-safe: readable_checkfile sits in arc/job/adapter.py so the adapter commits can use it without a forward dependency.

Full suite, serial: 2992 passed, 36 skipped, 5 failed. The 5 are arc/job/adapters/torch_ani_test.py, failing on an unconfigured TORCHANI_PYTHON (arc_env has no torchani; it lives in ts_gcn); nothing here touches that adapter.

🤖 Generated with Claude Code

[Copilot]

Copilot was unable to review this pull request because the user who requested the review has reached their quota limit.

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 65.02%. Comparing base (45d73a0) to head (afc47f9).
⚠️ Report is 1 commits behind head on main.

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #1014      +/-   ##
==========================================
+ Coverage   64.60%   65.02%   +0.41%     
==========================================
  Files         119      120       +1     
  Lines       39785    40462     +677     
  Branches    10307    10473     +166     
==========================================
+ Hits        25703    26310     +607     
- Misses      11105    11124      +19     
- Partials     2977     3028      +51     

Flag	Coverage Δ
functionaltests	65.02% <ø> (+0.41%)	⬆️
unittests	65.02% <ø> (+0.41%)	⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[alongd]

Reviewed at 32f60644aaa535eae16b01636c99ccf0fe595d17 against merge-base d285b1ed004bf3b7b5140b9da11986236ea4a8b1 — 33 files, +7378/−29, 16 source files across 9 packages.

First, the honest framing: this is careful work. The instrumentation is good, the adoption contract is stated precisely enough to be tested against, the fixtures are real ESS output rather than hand-made, and the PR description is unusually self-critical. Most of what follows is downstream of a single line that predates this branch, and which this PR's own docstring already explains. I've measured that one, so it isn't a hypothesis.

How to read the confidence numbers — they record how a claim was established, not how strongly it reads:

Score	Established by
10	Executed. Output pasted in.
8–9	Read the specific lines and traced the path.
6–7	Read and reasoned; nobody ran the failure.
≤5	Pattern match; check before acting.

---

Fix first

1. The adopted unrestricted reference collapses back to restricted — the feature is currently a no-op

arc/job/adapters/gaussian.py:405 · confidence 10 — code path executed, SCF consequence measured

guess=read is emitted whenever a checkfile exists and guess=mix otherwise, chosen with no reference to restricted/unrestricted. The u prefix is decided ~30 lines later. So a post-adoption job reads a restricted checkpoint into a UKS calculation.

Composed inputs, with an adopted verdict on a singlet TS:

sp    chk=True  adopted=True  restricted_used=False
        #P integral=(grid=ultrafine, Acc2E=12) guess=read ub3lyp/def2tzvp ...
opt   chk=True  adopted=True  restricted_used=False
        #P opt=(calcfc,maxcycle=100,noeigentest,ts)  guess=read ub3lyp/def2tzvp ...
sp    chk=False adopted=True  restricted_used=False
        #P integral=(grid=ultrafine, Acc2E=12) guess=mix ub3lyp/def2tzvp ...

The line is pre-existing at the merge-base — this PR doesn't introduce it, it collides with it:

$ git show d285b1ed0:arc/job/adapters/gaussian.py | grep -n "guess=read' if self.checkfile"
396:  input_dict['job_type_1'] += ' guess=read' if self.checkfile is not None and os.path.isfile(self.checkfile) \

And the collision is guaranteed on exactly the path that needs the other branch: run_stability_job:1812-1816 requires species.checkfile to exist and match the freq job's before a verdict is recorded at all. species.checkfile is set on the first opt/optfreq/composite (scheduler.py:1275), cleared only by ESS troubleshooting (:4314), and not cleared by switch_ts. So at every adoption a checkfile is present and guess=mix is unreachable.

You already state the mechanism, at scheduler.py:1747-1750:

without a checkfile the route would fall back to guess=mix, whose deliberately symmetry-broken SCF is a different wavefunction than the one under test

That reasoning is right, and it's applied to protect the measurement. It isn't applied to the jobs the verdict is adopted for.

The measurement

PySCF 2.14, B3LYP/def2-TZVP, grids level 4, conv_tol=1e-9, on the 15-atom C5H10 TS geometry taken from the Input orientation block of arc/testing/stability/rhf_uhf_instability_singlet_ts.out. 215 basis functions.

--- RKS reference ---
E(RKS)              = -196.490056185  converged=True

--- UKS seeded from the converged RKS density  [guess=read analogue] ---
E(UKS | from RKS)   = -196.490056186  converged=True  <S^2>=-0.000000

--- UKS from a HOMO/LUMO-mixed start  [guess=mix analogue] ---
E(UKS | from mix)   = -196.510081567  converged=True  <S^2>=0.863364

================ VERDICT ================
E(UKS from RKS) - E(RKS) = -4.724e-10 Ha  (-0.0000 kcal/mol)
E(UKS from mix) - E(RKS) = -2.003e-02 Ha  (-12.5661 kcal/mol)
UKS-from-RKS collapsed back to the restricted solution : True
UKS-from-mix found a lower, broken-symmetry solution   : True
energy the adoption path leaves on the table           : -12.5661 kcal/mol

Three things this establishes:

1. The setup is your fixture's. E(RKS) = −196.490056185 Ha against the fixture's own SCF Done: E(RB3LYP) = -196.490056259 (line 523) — agreement to 7.4e-8 Ha across two different programs and grids.
2. A UKS calculation seeded from a converged restricted density returns the restricted energy exactly — 4.7e-10 Ha, ⟨S²⟩ = 0.000000. It doesn't resemble the restricted solution; it is the restricted solution, carrying an unrestricted label.
3. A symmetry-broken start finds a genuine lower solution — 12.57 kcal/mol below, ⟨S²⟩ = 0.863. The RKS reference at this geometry really is externally unstable, exactly as the fixture's Gaussian verdict (Triplet-A at −0.0642 Ha) says.

The internal/external split agrees across both programs: PySCF's RKS.stability() returns (mo_i, mo_e, stable_i, stable_e) = (ndarray, None, True, None) — internal stability True, matching your stable Singlet-A, with the external test simply not run by default (that slot is None, not False; reading it as a negative is a trap I hit on the first pass).

Consequence. On this fixture, the adoption path produces −196.4901 Ha and records scf_reference: unrestricted, reference_mismatch: false, when the broken-symmetry solution it was reaching for is at −196.5101 Ha. The barrier is wrong by 12.57 kcal/mol, in the direction the PR set out to correct, with every record saying the correction was applied.

Caveat on provenance, stated plainly: this was measured in PySCF, not Gaussian. The mechanism is program-independent — the RKS solution is a stationary point of the UKS equations and α=β is preserved by the SCF — and the reference energy cross-validates against your own fixture. But PySCF's guess=read analogue is a hand-built density seed, not Gaussian's checkpoint reader. A g16 run would close that last gap.

Reproduce it

import numpy as np
from pyscf import dft, gto

GEOM = """
C   0.42653500   1.47789200  -0.00007700
C   0.49427800   0.01191500  -0.00000300
C   1.79225900  -0.71265700   0.00003600
C  -0.76952500  -0.70318200   0.00003500
C  -2.09552400  -0.11756100   0.00001100
H   1.39854400   1.97008900  -0.00010900
H  -0.16267200   1.83051900  -0.86690400
H   2.40600800  -0.45644800  -0.87511100
H   2.40600000  -0.45636800   0.87516500
H   1.65286800  -1.79552300   0.00008500
H  -0.70100100  -1.41388300   0.84851300
H  -0.70100600  -1.41396700  -0.84837200
H  -2.27706200   0.94604900  -0.00003500
H  -0.16265800   1.83060800   0.86672200
H  -2.94716600  -0.77952000   0.00003900
"""   # arc/testing/stability/rhf_uhf_instability_singlet_ts.out, first Input orientation block
mol = gto.M(atom=GEOM, basis='def2-tzvp', charge=0, spin=0, verbose=0, max_memory=16000)

mf_r = dft.RKS(mol); mf_r.xc = 'b3lyp'; mf_r.grids.level = 4
mf_r.conv_tol = 1e-9; mf_r.max_cycle = 200
e_r = mf_r.kernel()

def make_uks():
    mf = dft.UKS(mol); mf.xc = 'b3lyp'; mf.grids.level = 4
    mf.conv_tol = 1e-9; mf.max_cycle = 200
    return mf

# guess=read analogue: seed UKS with the converged restricted density
dm_r = mf_r.make_rdm1()
mf_u1 = make_uks()
e_u1 = mf_u1.kernel(dm0=np.array([dm_r / 2, dm_r / 2]))

# guess=mix analogue: rotate HOMO/LUMO oppositely in the two spin channels
nocc = mol.nelectron // 2
h, l = nocc - 1, nocc
c = 1.0 / np.sqrt(2.0)
homo, lumo = mf_r.mo_coeff[:, h].copy(), mf_r.mo_coeff[:, l].copy()
mo_a, mo_b = mf_r.mo_coeff.copy(), mf_r.mo_coeff.copy()
mo_a[:, h], mo_a[:, l] = c * (homo + lumo), c * (homo - lumo)
mo_b[:, h], mo_b[:, l] = c * (homo - lumo), c * (homo + lumo)
occ = np.zeros(mf_r.mo_coeff.shape[1]); occ[:nocc] = 1.0
mf_u2 = make_uks()
e_u2 = mf_u2.kernel(dm0=mf_u2.make_rdm1((mo_a, mo_b), (occ, occ)))

print(e_r, e_u1, mf_u1.spin_square()[0], e_u2, mf_u2.spin_square()[0])

Fix direction

On adoption emit guess=mix (or guess=(read,mix)), or add stable=opt to the first post-adoption job. Then assert s_squared > BROKEN_SYMMETRY_S2_THRESHOLD on that job's output and warn or re-run if it collapsed — otherwise the same silent failure returns the moment the guess logic is touched again.

Worth noting this is not reliably a uniform no-op: an unstable RKS solution is a saddle in orbital space, so DFT grid noise can seed the collapse in some species and not others. One campaign can yield a mix of restricted and broken-symmetry TSs, all labelled unrestricted. And it is invisible to the instrument this PR ships — _parse_spin_diagnostic (arc/output.py:334) reports s_squared: 0.0 against s_squared_expected: 0.0 for a multiplicity-1 species, i.e. perfectly healthy, for exactly this failure.

---

2. Five of the new tests are order-dependent and will flake CI red

arc/scheduler_test.py:737,768,781,796,809 · confidence 10

arc/common.py:186-189, inside initialize_log():

logging.addLevelName(logging.WARNING, 'Warning: ')   # process-global, permanent

The five new tests assert r.levelname == 'WARNING'. Once any earlier test in the same process constructs ARC(...) — e.g. arc/reaction/reaction_test.py:1054,1075 — level 30 is permanently named 'Warning: ' for the rest of the run.

$ pytest arc/scheduler_test.py -n0 -q
89 passed in 4.61s

$ pytest arc/reaction/reaction_test.py arc/scheduler_test.py -n0 -q
AssertionError: False is not true : no warning for an unstable TS: ['Warning: ', 'Warning: ']
5 failed, 115 passed in 16.18s

Makefile:56 runs -n auto, and xdist worker assignment isn't stable across runs, so this is a flaky red that won't reproduce in a scoped local run. grep -rn "levelname == 'WARNING'" arc/*_test.py matches only these five lines — the pattern is new to this PR.

This is also why the "zero failures" claim doesn't hold: full tree, single process, clean HEAD is 13 failed, 2974 passed, 10 skipped.

Fix: assert r.levelno == logging.WARNING (an immutable int).

---

3. Adoption is inert for the species it measures

arc/scheduler.py:3182 vs :3186 · confidence 9

if not switch_ts and species_has_sp(...):
    self.check_rxn_e0_by_spc(label)                     # 3182
if not switch_ts and (self.species_dict[label].is_ts
                      or job_scf_reference_is_restricted(job) is True):
    self.run_stability_job(label=label, freq_job=job)   # 3186

The barrier is computed and E0-checked before the stability job is spawned. The verdict lands asynchronously, and check_stability_job never re-runs or re-checks the TS. So the only channel through which adoption changes a number is a TS switch — and finding 1 says that channel is defeated.

Two independent passes reached this from opposite ends: the ordering above, and the observation that a TS guess's sp is typically already composed by the time its verdict lands.

Worth deciding explicitly whether adoption is meant to be retroactive. If yes, it needs a re-run path. If no, the description's "adopted for subsequent jobs" should name which subsequent jobs actually exist.

---

Then

4. A DFT-level verdict flips the reference at every other level and ESS

arc/job/adapters/common.py:479 · confidence 8

The stability job runs at freq_level (scheduler.py:1789) — B3LYP-class DFT, notoriously prone to RHF→UHF instability. The published energy comes from sp_level, typically CCSD(T)-F12 in Molpro or DLPNO-CCSD(T) in Orca. The gate consults nothing but is_ts and reference-agnostic method types, and every adapter reads it: orca.py:310 ('r' if is_restricted else 'u'), molpro.py:233, qchem.py:225, terachem.py:275.

Neither Orca nor Molpro has symmetry-breaking machinery in ARC's templates (Orca needs %scf BrokenSym/FlipSpin), so the flag is likely a no-op there too, while changing the algorithm and losing closed-shell efficiency. Magnitude unquantified — would need real output from those ESSs.

5. The mixed-reference check is structurally blind in the most common configuration

arc/scheduler.py:1544-1569 · confidence 8

When sp_level == opt_level, run_sp_job calls post_sp_actions directly and check_sp_job is never reached. record_scf_reference is called only from check_sp_job and check_freq_job, so scf_references['sp'] is permanently absent, check_scf_reference_consistency returns at its first guard forever, and reference_mismatch is permanently null. A single-level campaign cannot report a mixed reference under any circumstances.

Fix: record the reference in post_sp_actions, from the opt job's memo.

6. invalidates_analytic_freq has no consumer

arc/parser/adapters/gaussian.py:295-301 · confidence 10 (call-site grep: only tests, docstrings, and arc/output.py:248)

An internal instability — which check_stability_job itself describes as putting the analytic frequencies "outside the range in which they are defined" — produces a logger.warning, and then the ZPE and imaginary frequency go to Arkane unchanged.

Asymmetry worth closing: the mixed-reference case writes MIXED_SCF_REFERENCE_MESSAGE into output[label]['warnings']; the invalid-Hessian case writes nothing there, so it never reaches output.yml's warnings field.

7. Spin contamination is measured and not acted on

confidence 10

arc/testing/stability/stable_spin_contaminated_doublet_ts.out
stability: verdict 'stable', restricted False, invalidates_analytic_freq False
s2       : {'s_squared': 1.7488, 's_squared_expected': 0.75, 's_squared_annihilated': 0.774}

A doublet TS at 133% contamination passes as stable and its E0 goes to Arkane. Nothing compares s_squared against s_squared_expected outside output.yml. In practice this condition is far more common than an RHF→UHF instability, and the diagnostic you've already built could gate on it with two lines in check_freq_job.

8. arc/checks/spin.py ships three API traps, and no consumers

confidence 10 for each measurement

• s2_ls defaults to 0.0 (:52, :120) — silently wrong for every non-singlet target:

  BS doublet (S2=1.0) / HS quartet (3.80), target doublet:
    doublet target   : -100.00892857142857
    default (singlet): -100.03571428571428
    error kcal/mol   : -16.81
  triplet target: -50.0556 vs -50.3333  ->  17.4 kcal/mol
  

  s_squared_expected_from_multiplicity (arc/parser/parser.py:282) computes exactly this value and isn't reachable from the API that needs it. Drop the default; take the target multiplicity.
• MIN_S2_SEPARATION is a cliff (:103-106):

  sep=0.0999 -> None       (logger.debug only)
  sep=0.1001 -> -109.99000999000998    # a -9.99 Ha "correction" from E_BS = -100.0
  sep=0.11   -> -109.09090909090908
  sep=1.0    -> -101.0
  

  Refusing below the threshold is right and 0.1 is a defensible floor. The problem is the discontinuity, and that the first value past the guard is nonsense returned without a warning. Cap the amplification ratio rather than the denominator, and warn rather than debug when refusing.
• Nothing binds e_bs and e_hs to the same geometry, basis and functional (:48-53, :115-121); the returned record (:152-163) carries no level, geometry or path. The obvious future mistake — pulling E_HS from the triplet's own optimized geometry, which ARC has lying around — is silently wrong and unguarded. Take a Level and an xyz, and refuse mismatches.
• Two boundary mutants survived the suite: :103 <→<= and :150 >→>= both leave pytest arc/checks/spin_test.py at 26 passed. test_a_separation_at_the_floor_is_projected uses 0.2, not 0.1.

9. Three docstrings teach the wrong error direction

arc/job/adapters/common.py:394-396, repeated at scheduler.py:1859-1862 and in check_stability_job · confidence 10

They state the true low-spin energy lies between E_BS and E_RKS. Your own function refutes it:

E_restricted -100.0   E_BS -100.01   E_HS -99.99
E_projected  -100.03000000000002
projected below BS? True
is truth between E_BS and E_R? False
BS-R kcal/mol -6.28    proj-BS kcal/mol -12.55

Yamaguchi AP with E_BS < E_HS always gives E_proj < E_BS — below both, here by twice the R→BS gap. Adoption is an under-correction in a monotone direction, which is a better story than the one currently told; but whoever wires spin.py up will size the residual error at half its magnitude with the wrong sign.

10. stable=(rext,noopt) on an unrestricted reference runs only the internal test

arc/parser/adapters/gaussian.py:299 · confidence 10

All four shipped fixtures contain exactly one Stability analysis using <AA,BB:AA,BB> singles matrix and zero <AB,BA blocks — including both UB3LYP ones. So for an unrestricted reference the external (UHF→GHF) test never runs, yet check_stability_job's docstring and invalidates_analytic_freq treat a 'stable' verdict on an unrestricted reference as covering "any instability, internal or external". The verdict is narrower than the record claims.

For the restricted case your claim is confirmed: one matrix carries both Singlet-A (internal) and Triplet-A (external) roots, and the fixture's Triplet-A at −0.0642 Ha produces RHF -> UHF.

11. Two records of the same run disagree after a TS switch

arc/scheduler.py:4425-4433 · confidence 8

delete_all_species_jobs resets output[label]['paths'] and ['job_types'], so arc/output.py:657 reports wavefunction_stability: null for the final TS. But output[label]['wavefunction_stability'] (set at :1916) and output[label]['info'] (:1917) are top-level keys nothing clears, so arc/main.py:757 prints the abandoned guess's verdict against the surviving geometry. carry_stability_verdict_across_ts_switch was written to stop exactly this on the species object, and misses the output record.

A related consequence: the carried verdict drops the log path, so output.yml can show scf_reference.source: derived with log: null.

12. Order dependency in switch_ts is untested

arc/scheduler.py ~:3401-3416 · confidence 10 (mutation survived)

Swapping the call order of carry_stability_verdict_across_ts_switch and determine_most_likely_ts_conformer leaves pytest arc/scheduler_test.py -n0 -q at 89 passed. test_switch_ts_reduces_the_stability_verdict mocks out determine_most_likely_ts_conformer entirely, so chosen_ts never changes and the ordering is invisible. A future refactor that swaps them would record measured_on_ts_guess against the new guess instead of the abandoned one, with nothing catching it.

---

Unresolved — two passes reached opposite conclusions

Does the once-per-species spawn guard survive a restart?

arc/scheduler.py:1785 · confidence 6 — read-only on both sides, nobody ran a restart

if not self.job_types.get('stability', False) or 'stability' in self.job_dict[label].keys():
    return

• One reading: job_dict is rebuilt empty at Scheduler init (:467) and repopulated by restore_running_jobs only from restart_dict['running_jobs'] — a snapshot of running jobs. A completed stability job isn't in it, so after a restart the key is gone and a subsequent freq re-run spawns a duplicate.
• The other: no re-fire path exists, because a freq job already completed in a prior session is never resubmitted, and post_freq_actions only fires on a freq-completion event.

The disagreement is over whether a freq re-run can occur after the stability job has already completed, across a restart boundary. Impact either way is duplicate HPC spend, not a wrong number.

To settle: exercise a restart with stability: true and a freq re-run, or add self.output[label]['job_types'].get('stability') to the guard — that flag is already persisted in restart_dict['output']. Note the interaction you'd be choosing: delete_all_species_jobs resets it on a TS switch, so an output-based guard re-enables one spawn per TS guess. That may be what you want, but it's a change from the documented "once per species".

---

Housekeeping

Item	Location	Confidence
The six new fixtures carry the submitting account, hostname, compute-node name and scratch-path layout in their headers and archive blocks	arc/testing/{stability,spin}/*.out	10
Unguarded int() on a log-derived digit run; qchem.py already guards the same thing	gaussian.py:361, orca.py:226	10
parse_wavefunction_stability has no base-class default on ESSAdapter, unlike parse_s_squared/parse_ess_version; safe only via the getattr in make_parser	arc/parser/adapter.py	9
Returns: <type> inline instead of Google style	spin.py:23,90,143; parser.py:294; job/adapters/common.py:363; output.py:267,296	9
No Returns: section at all	run_stability_job, check_stability_job, log_open_shell_character_sources	9
Two import tempfile inside test bodies	parser_test.py:1258,1475	10
parse_s_squared implemented three times with a common skeleton	gaussian/orca/qchem	8

On the fixtures: this is the repo's existing norm rather than anything new — 13 tracked files at the merge-base carry the same host identifier and 6 the same account, and /gtmp is hardcoded in arc/settings/settings.py. No keys, tokens, emails or licence numbers are present. Scrubbing the six new ones is cheap and none of the new parsers key off those substrings; whether to backfill the pre-existing files is a separate question.

---

Checked and sound — please don't spend time re-verifying these

Chemistry

• The Yamaguchi formula is correct. Denominator is ⟨S²⟩_HS − ⟨S²⟩_BS, derived from E_BS = (1−w)E_LS + wE_HS with w = (S²_BS−S²_LS)/(S²_HS−S²_LS). Verified algebraically and numerically. The s2_ls plumbing is complete — it reaches the formula, the ordering guard (:107) and broken_symmetry (:150).
• Before-annihilation ⟨S²⟩ is right and consistently fed. s_squared: 0.7536 from the <Sx>= line vs s_squared_annihilated: 0.75. The Initial guess exclusion works (uhf_died_before_scf_septet.out → None, not 12.0) and the first-Multiplicity rule survives guess=fragment (→ expected 0.75, not the fragment's 1).
• AEC/BAC is not a new problem. arc/statmech/arkane.py:373-398 applies corrections uniformly; a TS's energy enters only via E0(TS) − E0(reactants), where atom counts are identical and AEC cancels exactly. Wells are never adopted. Your stated reason for excluding wells is sound, and the AEC half doesn't additionally condemn the TS case.
• The composite/semiempirical early return should be closed, not carried as a follow-up. uCBS-QB3 isn't a Gaussian keyword — CBS-QB3 selects its open-shell variant from the multiplicity itself, as do AM1/PM6 — and run_stability_job:1802-1806 refuses any non-DFT/HF level, so a derived verdict can't exist for a composite species anyway.
• Your guess=mix re-derivation holds. arc/testing/restart/1_restart_thermo/calcs/freq_a19031.out: route line 84 #P guess=mix wb97xd/def2tzvp freq, line 105 Multiplicity = 1, line 223 Mixing orbitals, IMix= 1 ... Coef= 7.07106781D-01 7.07106781D-01, line 422 SCF Done: E(RwB97XD). g16 does apply mixing on a restricted reference. Note this is also the fact finding 1 turns on.
• The fixtures are genuine — g16 RevC.02, real l1.exe paths, real archive blocks, real timestamps.

Persistence

• The provenance does survive a restart — via restart.yml / ARCSpecies.as_dict() (species.py:810-813), read back at :947-949. Not via output.yml, which is confirmed write-only/export. The claim holds; the description's stated mechanism is the part to correct.
• Old pre-PR restart fixtures resume cleanly. All 4 species of 1_restart_thermo/restart.yml constructed; the three new attributes defaulted to None/{}/None. No KeyError/TypeError/AttributeError.
• No stability respawn on resumed old projects — default_job_types['stability'] is False and initialize_job_types fills the missing key from defaults_to_false.
• The number_of_radicals provenance question is moot. Only two assignments exist in arc/: the constructor default (species.py:382) and the restart read (:946). ARC never writes it programmatically, so a restored value is exactly as user-declared as a fresh one.
• No tuple→list degradation. In real usage relaxations/negative_eigenvectors are always lists.
• ts_checks semantics untouched by this diff.
• restricted_used round-trips safely — popped in restore_running_jobs before job_factory.

Contracts

• ESSAdapter.parse_s_squared has a safe return None default (arc/parser/adapter.py:210), not @abstractmethod. All 9 subclasses probed; the 6 without an override inherit cleanly, no AttributeError.
• No JobAdapter abstract method added, so no adapter is forced to implement anything new. job_type='stability' is only ever constructed by run_stability_job, which gates on job_adapter != 'gaussian'.
• The default_job_types in-place-mutation trap does not bite — arc/common.py:92 lists stability in defaults_to_false; zero import-time default_job_types[key] subscripts anywhere.
• specific_job_type: stability runs nothing, confirmed by execution, matching your own admission — and it does not also hang, since check_all_done (:3843-3845) skips stability. It's the only such trap this PR introduces.
• check_negative_freq's return-shape change is fully absorbed — both call sites unpack, all 9 returns in the body are 2-tuples.

Security / performance

• No ReDoS in any of the nine new regexes at n = 200 / 2,000 / 20,000 / 200,000 — nothing over 50 ms. Both parsers gate on cheap substring checks before invoking re.search.
• No eval / exec / pickle / yaml.load / subprocess / shell=True in the new source.
• No log-derived string reaches an input deck. The stability route section (gaussian.py:333-336) is composed exclusively of module constants; run_stability_job passes only xyz and level.
• No path traversal — the only path written to output.yml is ARC-constructed, never parsed out of a log.
• Checkfile warm-start confirmed — each stability job is an SCF restart plus a matrix diagonalization, not a re-optimization. (Ironically, the same mechanism as finding 1.)
• No completed job is invalidated or re-queued on adoption.
• No per-tick re-parse — check_stability_job is called once, from the single job-completion branch at :869.
• Parser cost at scale: a synthetic 210 MB log gives 1.73 s / 392 MB peak for parse_wavefunction_stability and 1.51 s / 394 MB for parse_s_squared — the same readlines() pattern every existing ARC parser uses, and real logs are orders smaller.
• Parser false positives swept: across all tracked .out/.log fixtures, only the four new stability fixtures parsed — parsed_count 4, errors_count 0.
• stability cannot block convergence and is correctly reset by delete_all_species_jobs.

Mutants the suite killed — i.e. the tests work:

• 'unknown' → 'stable' on an unreadable verdict: killed by test_unparsed_verdict_is_not_reported_as_stable.
• number_of_radicals is None → truthiness: killed, confirming declared-0 is distinguished from undeclared-None.
• number_of_radicals > 1 → >= 1: killed, confirming a declared 1 isn't credited as an open-shell source.
• Q-Chem len(tokens) >= 2 → > 2: killed by test_parse_s_squared_qchem.

---

Where this leaves the PR

Two chains account for nearly everything above.

Finding 1 is the review. guess=read from a restricted checkpoint returns E(RKS) exactly, so the adoption path changes no number today. Findings 3, 10 and 11 are downstream detail of the same defect: the feature measures correctly, records correctly, and then applies a reference flip that the SCF quietly undoes. Fix finding 1 and the PR does what its description says.

arc/checks/spin.py has no consumers, which is why findings 8 and 9 are latent rather than live. Every trap in it is harmless today and becomes a wrong published number the day someone wires it up. Fixing the API now costs a fraction of debugging it later. The same holds for finding 6 — invalidates_analytic_freq is computed and discarded.

A smaller third chain: findings 5, 6, 7 and 11 are all the same shape — a diagnostic that is computed, recorded, and never gated on. Acting narrowly by design is defensible, but four separate measurements currently reach no decision, and finding 7 is the one most likely to matter in a real campaign.

Method, for what it's worth

Eight parallel review passes over the same pinned diff — testing/mutation, security, API-contract, maintainability, performance, data-migration, an adversarial pass, and a cross-model adversarial pass — none of which saw each other's conclusions, plus the PySCF measurement above. Findings were deduplicated by file:line and by mechanism; cross-model findings were re-derived against the source before being accepted.

Three conclusions were reached independently by two passes each, which is why I'd weight them: adoption being inert for the measured species (finding 3), the provenance riding restart.yml rather than output.yml, and the TS-switch record losing its log path (finding 11). The restart spawn guard is recorded as unresolved precisely because two passes read the same lines and disagreed.

All test gates ran in arc_env with an empty HOME and pytest -n0. Known-failing at any SHA on my box, not attributable to this branch: checks/ts_test.py::test_check_rxn_e0, ::test_compute_rxn_e0, processor_test.py::test_compare_rates, statmech/arkane_test.py::test_run_statmech_using_molecular_properties (missing arkane in rmg_env), plus molecule_test.py::TestConnectTheDots.

No ESS job was submitted at any point.

[calvinp0]

Reviewed at 32f60644aaa535eae16b01636c99ccf0fe595d17 against merge-base d285b1ed004bf3b7b5140b9da11986236ea4a8b1 — 33 files, +7378/−29, 16 source files across 9 packages.

First, the honest framing: this is careful work. The instrumentation is good, the adoption contract is stated precisely enough to be tested against, the fixtures are real ESS output rather than hand-made, and the PR description is unusually self-critical. Most of what follows is downstream of a single line that predates this branch, and which this PR's own docstring already explains. I've measured that one, so it isn't a hypothesis.

How to read the confidence numbers — they record how a claim was established, not how strongly it reads:
Score Established by
10 Executed. Output pasted in.
8–9 Read the specific lines and traced the path.
6–7 Read and reasoned; nobody ran the failure.
≤5 Pattern match; check before acting.

Fix first

1. The adopted unrestricted reference collapses back to restricted — the feature is currently a no-op

arc/job/adapters/gaussian.py:405 · confidence 10 — code path executed, SCF consequence measured

guess=read is emitted whenever a checkfile exists and guess=mix otherwise, chosen with no reference to restricted/unrestricted. The u prefix is decided ~30 lines later. So a post-adoption job reads a restricted checkpoint into a UKS calculation.

Composed inputs, with an adopted verdict on a singlet TS:

sp    chk=True  adopted=True  restricted_used=False
        #P integral=(grid=ultrafine, Acc2E=12) guess=read ub3lyp/def2tzvp ...
opt   chk=True  adopted=True  restricted_used=False
        #P opt=(calcfc,maxcycle=100,noeigentest,ts)  guess=read ub3lyp/def2tzvp ...
sp    chk=False adopted=True  restricted_used=False
        #P integral=(grid=ultrafine, Acc2E=12) guess=mix ub3lyp/def2tzvp ...

The line is pre-existing at the merge-base — this PR doesn't introduce it, it collides with it:

$ git show d285b1ed0:arc/job/adapters/gaussian.py | grep -n "guess=read' if self.checkfile"
396:  input_dict['job_type_1'] += ' guess=read' if self.checkfile is not None and os.path.isfile(self.checkfile) \

And the collision is guaranteed on exactly the path that needs the other branch: run_stability_job:1812-1816 requires species.checkfile to exist and match the freq job's before a verdict is recorded at all. species.checkfile is set on the first opt/optfreq/composite (scheduler.py:1275), cleared only by ESS troubleshooting (:4314), and not cleared by switch_ts. So at every adoption a checkfile is present and guess=mix is unreachable.

You already state the mechanism, at scheduler.py:1747-1750:

without a checkfile the route would fall back to guess=mix, whose deliberately symmetry-broken SCF is a different wavefunction than the one under test

That reasoning is right, and it's applied to protect the measurement. It isn't applied to the jobs the verdict is adopted for.

The measurement

PySCF 2.14, B3LYP/def2-TZVP, grids level 4, conv_tol=1e-9, on the 15-atom C5H10 TS geometry taken from the Input orientation block of arc/testing/stability/rhf_uhf_instability_singlet_ts.out. 215 basis functions.

--- RKS reference ---
E(RKS)              = -196.490056185  converged=True

--- UKS seeded from the converged RKS density  [guess=read analogue] ---
E(UKS | from RKS)   = -196.490056186  converged=True  <S^2>=-0.000000

--- UKS from a HOMO/LUMO-mixed start  [guess=mix analogue] ---
E(UKS | from mix)   = -196.510081567  converged=True  <S^2>=0.863364

================ VERDICT ================
E(UKS from RKS) - E(RKS) = -4.724e-10 Ha  (-0.0000 kcal/mol)
E(UKS from mix) - E(RKS) = -2.003e-02 Ha  (-12.5661 kcal/mol)
UKS-from-RKS collapsed back to the restricted solution : True
UKS-from-mix found a lower, broken-symmetry solution   : True
energy the adoption path leaves on the table           : -12.5661 kcal/mol

Three things this establishes:

1. **The setup is your fixture's.** E(RKS) = −196.490056185 Ha against the fixture's own `SCF Done: E(RB3LYP) = -196.490056259` (line 523) — agreement to 7.4e-8 Ha across two different programs and grids.

2. **A UKS calculation seeded from a converged restricted density returns the restricted energy exactly** — 4.7e-10 Ha, ⟨S²⟩ = 0.000000. It doesn't resemble the restricted solution; it _is_ the restricted solution, carrying an unrestricted label.

3. **A symmetry-broken start finds a genuine lower solution** — 12.57 kcal/mol below, ⟨S²⟩ = 0.863. The RKS reference at this geometry really is externally unstable, exactly as the fixture's Gaussian verdict (Triplet-A at −0.0642 Ha) says.

The internal/external split agrees across both programs: PySCF's RKS.stability() returns (mo_i, mo_e, stable_i, stable_e) = (ndarray, None, True, None) — internal stability True, matching your stable Singlet-A, with the external test simply not run by default (that slot is None, not False; reading it as a negative is a trap I hit on the first pass).

Consequence. On this fixture, the adoption path produces −196.4901 Ha and records scf_reference: unrestricted, reference_mismatch: false, when the broken-symmetry solution it was reaching for is at −196.5101 Ha. The barrier is wrong by 12.57 kcal/mol, in the direction the PR set out to correct, with every record saying the correction was applied.

Caveat on provenance, stated plainly: this was measured in PySCF, not Gaussian. The mechanism is program-independent — the RKS solution is a stationary point of the UKS equations and α=β is preserved by the SCF — and the reference energy cross-validates against your own fixture. But PySCF's guess=read analogue is a hand-built density seed, not Gaussian's checkpoint reader. A g16 run would close that last gap.

Reproduce it

import numpy as np
from pyscf import dft, gto

GEOM = """
C   0.42653500   1.47789200  -0.00007700
C   0.49427800   0.01191500  -0.00000300
C   1.79225900  -0.71265700   0.00003600
C  -0.76952500  -0.70318200   0.00003500
C  -2.09552400  -0.11756100   0.00001100
H   1.39854400   1.97008900  -0.00010900
H  -0.16267200   1.83051900  -0.86690400
H   2.40600800  -0.45644800  -0.87511100
H   2.40600000  -0.45636800   0.87516500
H   1.65286800  -1.79552300   0.00008500
H  -0.70100100  -1.41388300   0.84851300
H  -0.70100600  -1.41396700  -0.84837200
H  -2.27706200   0.94604900  -0.00003500
H  -0.16265800   1.83060800   0.86672200
H  -2.94716600  -0.77952000   0.00003900
"""   # arc/testing/stability/rhf_uhf_instability_singlet_ts.out, first Input orientation block
mol = gto.M(atom=GEOM, basis='def2-tzvp', charge=0, spin=0, verbose=0, max_memory=16000)

mf_r = dft.RKS(mol); mf_r.xc = 'b3lyp'; mf_r.grids.level = 4
mf_r.conv_tol = 1e-9; mf_r.max_cycle = 200
e_r = mf_r.kernel()

def make_uks():
    mf = dft.UKS(mol); mf.xc = 'b3lyp'; mf.grids.level = 4
    mf.conv_tol = 1e-9; mf.max_cycle = 200
    return mf

# guess=read analogue: seed UKS with the converged restricted density
dm_r = mf_r.make_rdm1()
mf_u1 = make_uks()
e_u1 = mf_u1.kernel(dm0=np.array([dm_r / 2, dm_r / 2]))

# guess=mix analogue: rotate HOMO/LUMO oppositely in the two spin channels
nocc = mol.nelectron // 2
h, l = nocc - 1, nocc
c = 1.0 / np.sqrt(2.0)
homo, lumo = mf_r.mo_coeff[:, h].copy(), mf_r.mo_coeff[:, l].copy()
mo_a, mo_b = mf_r.mo_coeff.copy(), mf_r.mo_coeff.copy()
mo_a[:, h], mo_a[:, l] = c * (homo + lumo), c * (homo - lumo)
mo_b[:, h], mo_b[:, l] = c * (homo - lumo), c * (homo + lumo)
occ = np.zeros(mf_r.mo_coeff.shape[1]); occ[:nocc] = 1.0
mf_u2 = make_uks()
e_u2 = mf_u2.kernel(dm0=mf_u2.make_rdm1((mo_a, mo_b), (occ, occ)))

print(e_r, e_u1, mf_u1.spin_square()[0], e_u2, mf_u2.spin_square()[0])

Fix direction

On adoption emit guess=mix (or guess=(read,mix)), or add stable=opt to the first post-adoption job. Then assert s_squared > BROKEN_SYMMETRY_S2_THRESHOLD on that job's output and warn or re-run if it collapsed — otherwise the same silent failure returns the moment the guess logic is touched again.

Worth noting this is not reliably a uniform no-op: an unstable RKS solution is a saddle in orbital space, so DFT grid noise can seed the collapse in some species and not others. One campaign can yield a mix of restricted and broken-symmetry TSs, all labelled unrestricted. And it is invisible to the instrument this PR ships — _parse_spin_diagnostic (arc/output.py:334) reports s_squared: 0.0 against s_squared_expected: 0.0 for a multiplicity-1 species, i.e. perfectly healthy, for exactly this failure.

2. Five of the new tests are order-dependent and will flake CI red

arc/scheduler_test.py:737,768,781,796,809 · confidence 10

arc/common.py:186-189, inside initialize_log():

logging.addLevelName(logging.WARNING, 'Warning: ')   # process-global, permanent

The five new tests assert r.levelname == 'WARNING'. Once any earlier test in the same process constructs ARC(...) — e.g. arc/reaction/reaction_test.py:1054,1075 — level 30 is permanently named 'Warning: ' for the rest of the run.

$ pytest arc/scheduler_test.py -n0 -q
89 passed in 4.61s

$ pytest arc/reaction/reaction_test.py arc/scheduler_test.py -n0 -q
AssertionError: False is not true : no warning for an unstable TS: ['Warning: ', 'Warning: ']
5 failed, 115 passed in 16.18s

Makefile:56 runs -n auto, and xdist worker assignment isn't stable across runs, so this is a flaky red that won't reproduce in a scoped local run. grep -rn "levelname == 'WARNING'" arc/*_test.py matches only these five lines — the pattern is new to this PR.

This is also why the "zero failures" claim doesn't hold: full tree, single process, clean HEAD is 13 failed, 2974 passed, 10 skipped.

Fix: assert r.levelno == logging.WARNING (an immutable int).

3. Adoption is inert for the species it measures

arc/scheduler.py:3182 vs :3186 · confidence 9

if not switch_ts and species_has_sp(...):
    self.check_rxn_e0_by_spc(label)                     # 3182
if not switch_ts and (self.species_dict[label].is_ts
                      or job_scf_reference_is_restricted(job) is True):
    self.run_stability_job(label=label, freq_job=job)   # 3186

The barrier is computed and E0-checked before the stability job is spawned. The verdict lands asynchronously, and check_stability_job never re-runs or re-checks the TS. So the only channel through which adoption changes a number is a TS switch — and finding 1 says that channel is defeated.

Two independent passes reached this from opposite ends: the ordering above, and the observation that a TS guess's sp is typically already composed by the time its verdict lands.

Worth deciding explicitly whether adoption is meant to be retroactive. If yes, it needs a re-run path. If no, the description's "adopted for subsequent jobs" should name which subsequent jobs actually exist.

Then

4. A DFT-level verdict flips the reference at every other level and ESS

arc/job/adapters/common.py:479 · confidence 8

The stability job runs at freq_level (scheduler.py:1789) — B3LYP-class DFT, notoriously prone to RHF→UHF instability. The published energy comes from sp_level, typically CCSD(T)-F12 in Molpro or DLPNO-CCSD(T) in Orca. The gate consults nothing but is_ts and reference-agnostic method types, and every adapter reads it: orca.py:310 ('r' if is_restricted else 'u'), molpro.py:233, qchem.py:225, terachem.py:275.

Neither Orca nor Molpro has symmetry-breaking machinery in ARC's templates (Orca needs %scf BrokenSym/FlipSpin), so the flag is likely a no-op there too, while changing the algorithm and losing closed-shell efficiency. Magnitude unquantified — would need real output from those ESSs.

5. The mixed-reference check is structurally blind in the most common configuration

arc/scheduler.py:1544-1569 · confidence 8

When sp_level == opt_level, run_sp_job calls post_sp_actions directly and check_sp_job is never reached. record_scf_reference is called only from check_sp_job and check_freq_job, so scf_references['sp'] is permanently absent, check_scf_reference_consistency returns at its first guard forever, and reference_mismatch is permanently null. A single-level campaign cannot report a mixed reference under any circumstances.

Fix: record the reference in post_sp_actions, from the opt job's memo.

6. invalidates_analytic_freq has no consumer

arc/parser/adapters/gaussian.py:295-301 · confidence 10 (call-site grep: only tests, docstrings, and arc/output.py:248)

An internal instability — which check_stability_job itself describes as putting the analytic frequencies "outside the range in which they are defined" — produces a logger.warning, and then the ZPE and imaginary frequency go to Arkane unchanged.

Asymmetry worth closing: the mixed-reference case writes MIXED_SCF_REFERENCE_MESSAGE into output[label]['warnings']; the invalid-Hessian case writes nothing there, so it never reaches output.yml's warnings field.

7. Spin contamination is measured and not acted on

confidence 10

arc/testing/stability/stable_spin_contaminated_doublet_ts.out
stability: verdict 'stable', restricted False, invalidates_analytic_freq False
s2       : {'s_squared': 1.7488, 's_squared_expected': 0.75, 's_squared_annihilated': 0.774}

A doublet TS at 133% contamination passes as stable and its E0 goes to Arkane. Nothing compares s_squared against s_squared_expected outside output.yml. In practice this condition is far more common than an RHF→UHF instability, and the diagnostic you've already built could gate on it with two lines in check_freq_job.

8. arc/checks/spin.py ships three API traps, and no consumers

confidence 10 for each measurement

* **`s2_ls` defaults to `0.0`** (`:52`, `:120`) — silently wrong for every non-singlet target:
  ```
  BS doublet (S2=1.0) / HS quartet (3.80), target doublet:
    doublet target   : -100.00892857142857
    default (singlet): -100.03571428571428
    error kcal/mol   : -16.81
  triplet target: -50.0556 vs -50.3333  ->  17.4 kcal/mol
  ```
  
  
      
        
      
  
        
      
  
      
    
  `s_squared_expected_from_multiplicity` (`arc/parser/parser.py:282`) computes exactly this value and isn't reachable from the API that needs it. Drop the default; take the target multiplicity.

* **`MIN_S2_SEPARATION` is a cliff** (`:103-106`):
  ```
  sep=0.0999 -> None       (logger.debug only)
  sep=0.1001 -> -109.99000999000998    # a -9.99 Ha "correction" from E_BS = -100.0
  sep=0.11   -> -109.09090909090908
  sep=1.0    -> -101.0
  ```
  
  
      
        
      
  
        
      
  
      
    
  Refusing below the threshold is right and 0.1 is a defensible floor. The problem is the discontinuity, and that the first value past the guard is nonsense returned without a warning. Cap the amplification ratio rather than the denominator, and warn rather than debug when refusing.

* **Nothing binds `e_bs` and `e_hs` to the same geometry, basis and functional** (`:48-53`, `:115-121`); the returned record (`:152-163`) carries no level, geometry or path. The obvious future mistake — pulling E_HS from the triplet's own optimized geometry, which ARC has lying around — is silently wrong and unguarded. Take a `Level` and an xyz, and refuse mismatches.

* **Two boundary mutants survived the suite**: `:103` `<`→`<=` and `:150` `>`→`>=` both leave `pytest arc/checks/spin_test.py` at `26 passed`. `test_a_separation_at_the_floor_is_projected` uses 0.2, not 0.1.

9. Three docstrings teach the wrong error direction

arc/job/adapters/common.py:394-396, repeated at scheduler.py:1859-1862 and in check_stability_job · confidence 10

They state the true low-spin energy lies between E_BS and E_RKS. Your own function refutes it:

E_restricted -100.0   E_BS -100.01   E_HS -99.99
E_projected  -100.03000000000002
projected below BS? True
is truth between E_BS and E_R? False
BS-R kcal/mol -6.28    proj-BS kcal/mol -12.55

Yamaguchi AP with E_BS < E_HS always gives E_proj < E_BS — below both, here by twice the R→BS gap. Adoption is an under-correction in a monotone direction, which is a better story than the one currently told; but whoever wires spin.py up will size the residual error at half its magnitude with the wrong sign.

10. stable=(rext,noopt) on an unrestricted reference runs only the internal test

arc/parser/adapters/gaussian.py:299 · confidence 10

All four shipped fixtures contain exactly one Stability analysis using <AA,BB:AA,BB> singles matrix and zero <AB,BA blocks — including both UB3LYP ones. So for an unrestricted reference the external (UHF→GHF) test never runs, yet check_stability_job's docstring and invalidates_analytic_freq treat a 'stable' verdict on an unrestricted reference as covering "any instability, internal or external". The verdict is narrower than the record claims.

For the restricted case your claim is confirmed: one matrix carries both Singlet-A (internal) and Triplet-A (external) roots, and the fixture's Triplet-A at −0.0642 Ha produces RHF -> UHF.

11. Two records of the same run disagree after a TS switch

arc/scheduler.py:4425-4433 · confidence 8

delete_all_species_jobs resets output[label]['paths'] and ['job_types'], so arc/output.py:657 reports wavefunction_stability: null for the final TS. But output[label]['wavefunction_stability'] (set at :1916) and output[label]['info'] (:1917) are top-level keys nothing clears, so arc/main.py:757 prints the abandoned guess's verdict against the surviving geometry. carry_stability_verdict_across_ts_switch was written to stop exactly this on the species object, and misses the output record.

A related consequence: the carried verdict drops the log path, so output.yml can show scf_reference.source: derived with log: null.

12. Order dependency in switch_ts is untested

arc/scheduler.py ~:3401-3416 · confidence 10 (mutation survived)

Swapping the call order of carry_stability_verdict_across_ts_switch and determine_most_likely_ts_conformer leaves pytest arc/scheduler_test.py -n0 -q at 89 passed. test_switch_ts_reduces_the_stability_verdict mocks out determine_most_likely_ts_conformer entirely, so chosen_ts never changes and the ordering is invisible. A future refactor that swaps them would record measured_on_ts_guess against the new guess instead of the abandoned one, with nothing catching it.

Unresolved — two passes reached opposite conclusions

Does the once-per-species spawn guard survive a restart?

arc/scheduler.py:1785 · confidence 6 — read-only on both sides, nobody ran a restart

if not self.job_types.get('stability', False) or 'stability' in self.job_dict[label].keys():
    return

* **One reading:** `job_dict` is rebuilt empty at Scheduler init (`:467`) and repopulated by `restore_running_jobs` only from `restart_dict['running_jobs']` — a snapshot of _running_ jobs. A completed stability job isn't in it, so after a restart the key is gone and a subsequent freq re-run spawns a duplicate.

* **The other:** no re-fire path exists, because a freq job already completed in a prior session is never resubmitted, and `post_freq_actions` only fires on a freq-completion event.

The disagreement is over whether a freq re-run can occur after the stability job has already completed, across a restart boundary. Impact either way is duplicate HPC spend, not a wrong number.

To settle: exercise a restart with stability: true and a freq re-run, or add self.output[label]['job_types'].get('stability') to the guard — that flag is already persisted in restart_dict['output']. Note the interaction you'd be choosing: delete_all_species_jobs resets it on a TS switch, so an output-based guard re-enables one spawn per TS guess. That may be what you want, but it's a change from the documented "once per species".

Housekeeping

Item Location Confidence
The six new fixtures carry the submitting account, hostname, compute-node name and scratch-path layout in their headers and archive blocks arc/testing/{stability,spin}/*.out 10
Unguarded int() on a log-derived digit run; qchem.py already guards the same thing gaussian.py:361, orca.py:226 10
parse_wavefunction_stability has no base-class default on ESSAdapter, unlike parse_s_squared/parse_ess_version; safe only via the getattr in make_parser arc/parser/adapter.py 9
Returns: <type> inline instead of Google style spin.py:23,90,143; parser.py:294; job/adapters/common.py:363; output.py:267,296 9
No Returns: section at all run_stability_job, check_stability_job, log_open_shell_character_sources 9
Two import tempfile inside test bodies parser_test.py:1258,1475 10
parse_s_squared implemented three times with a common skeleton gaussian/orca/qchem 8

On the fixtures: this is the repo's existing norm rather than anything new — 13 tracked files at the merge-base carry the same host identifier and 6 the same account, and /gtmp is hardcoded in arc/settings/settings.py. No keys, tokens, emails or licence numbers are present. Scrubbing the six new ones is cheap and none of the new parsers key off those substrings; whether to backfill the pre-existing files is a separate question.

Checked and sound — please don't spend time re-verifying these

Chemistry

* **The Yamaguchi formula is correct.** Denominator is ⟨S²⟩_HS − ⟨S²⟩_BS, derived from E_BS = (1−w)E_LS + wE_HS with w = (S²_BS−S²_LS)/(S²_HS−S²_LS). Verified algebraically and numerically. The `s2_ls` plumbing is complete — it reaches the formula, the ordering guard (`:107`) and `broken_symmetry` (`:150`).

* **Before-annihilation ⟨S²⟩ is right and consistently fed.** `s_squared: 0.7536` from the `<Sx>=` line vs `s_squared_annihilated: 0.75`. The `Initial guess` exclusion works (`uhf_died_before_scf_septet.out` → `None`, not 12.0) and the first-`Multiplicity` rule survives `guess=fragment` (→ expected 0.75, not the fragment's 1).

* **AEC/BAC is not a new problem.** `arc/statmech/arkane.py:373-398` applies corrections uniformly; a TS's energy enters only via E0(TS) − E0(reactants), where atom counts are identical and AEC cancels exactly. Wells are never adopted. Your stated reason for excluding wells is sound, and the AEC half doesn't additionally condemn the TS case.

* **The composite/semiempirical early return should be closed, not carried as a follow-up.** `uCBS-QB3` isn't a Gaussian keyword — `CBS-QB3` selects its open-shell variant from the multiplicity itself, as do AM1/PM6 — and `run_stability_job:1802-1806` refuses any non-DFT/HF level, so a derived verdict can't exist for a composite species anyway.

* **Your `guess=mix` re-derivation holds.** `arc/testing/restart/1_restart_thermo/calcs/freq_a19031.out`: route line 84 `#P guess=mix wb97xd/def2tzvp freq`, line 105 `Multiplicity = 1`, line 223 `Mixing orbitals, IMix= 1 ... Coef= 7.07106781D-01 7.07106781D-01`, line 422 `SCF Done: E(RwB97XD)`. g16 does apply mixing on a restricted reference. Note this is also the fact finding 1 turns on.

* **The fixtures are genuine** — g16 RevC.02, real `l1.exe` paths, real archive blocks, real timestamps.

Persistence

* **The provenance does survive a restart** — via `restart.yml` / `ARCSpecies.as_dict()` (`species.py:810-813`), read back at `:947-949`. Not via `output.yml`, which is confirmed write-only/export. The claim holds; the description's stated mechanism is the part to correct.

* **Old pre-PR restart fixtures resume cleanly.** All 4 species of `1_restart_thermo/restart.yml` constructed; the three new attributes defaulted to `None`/`{}`/`None`. No `KeyError`/`TypeError`/`AttributeError`.

* **No stability respawn on resumed old projects** — `default_job_types['stability']` is `False` and `initialize_job_types` fills the missing key from `defaults_to_false`.

* **The `number_of_radicals` provenance question is moot.** Only two assignments exist in `arc/`: the constructor default (`species.py:382`) and the restart read (`:946`). ARC never writes it programmatically, so a restored value is exactly as user-declared as a fresh one.

* **No tuple→list degradation.** In real usage `relaxations`/`negative_eigenvectors` are always lists.

* **`ts_checks` semantics untouched** by this diff.

* **`restricted_used` round-trips safely** — popped in `restore_running_jobs` before `job_factory`.

Contracts

* **`ESSAdapter.parse_s_squared`** has a safe `return None` default (`arc/parser/adapter.py:210`), not `@abstractmethod`. All 9 subclasses probed; the 6 without an override inherit cleanly, no `AttributeError`.

* **No `JobAdapter` abstract method added**, so no adapter is forced to implement anything new. `job_type='stability'` is only ever constructed by `run_stability_job`, which gates on `job_adapter != 'gaussian'`.

* **The `default_job_types` in-place-mutation trap does not bite** — `arc/common.py:92` lists `stability` in `defaults_to_false`; zero import-time `default_job_types[key]` subscripts anywhere.

* **`specific_job_type: stability` runs nothing**, confirmed by execution, matching your own admission — and it does _not_ also hang, since `check_all_done` (`:3843-3845`) skips `stability`. It's the only such trap this PR introduces.

* **`check_negative_freq`'s return-shape change is fully absorbed** — both call sites unpack, all 9 returns in the body are 2-tuples.

Security / performance

* **No ReDoS** in any of the nine new regexes at n = 200 / 2,000 / 20,000 / 200,000 — nothing over 50 ms. Both parsers gate on cheap substring checks before invoking `re.search`.

* **No `eval` / `exec` / `pickle` / `yaml.load` / `subprocess` / `shell=True`** in the new source.

* **No log-derived string reaches an input deck.** The stability route section (`gaussian.py:333-336`) is composed exclusively of module constants; `run_stability_job` passes only `xyz` and `level`.

* **No path traversal** — the only path written to `output.yml` is ARC-constructed, never parsed out of a log.

* **Checkfile warm-start confirmed** — each stability job is an SCF restart plus a matrix diagonalization, not a re-optimization. (Ironically, the same mechanism as finding 1.)

* **No completed job is invalidated or re-queued on adoption.**

* **No per-tick re-parse** — `check_stability_job` is called once, from the single job-completion branch at `:869`.

* **Parser cost at scale**: a synthetic 210 MB log gives 1.73 s / 392 MB peak for `parse_wavefunction_stability` and 1.51 s / 394 MB for `parse_s_squared` — the same `readlines()` pattern every existing ARC parser uses, and real logs are orders smaller.

* **Parser false positives swept**: across all tracked `.out`/`.log` fixtures, only the four new stability fixtures parsed — `parsed_count 4, errors_count 0`.

* **`stability` cannot block convergence** and is correctly reset by `delete_all_species_jobs`.

Mutants the suite killed — i.e. the tests work:

* `'unknown'` → `'stable'` on an unreadable verdict: killed by `test_unparsed_verdict_is_not_reported_as_stable`.

* `number_of_radicals is None` → truthiness: killed, confirming declared-`0` is distinguished from undeclared-`None`.

* `number_of_radicals > 1` → `>= 1`: killed, confirming a declared `1` isn't credited as an open-shell source.

* Q-Chem `len(tokens) >= 2` → `> 2`: killed by `test_parse_s_squared_qchem`.

Where this leaves the PR

Two chains account for nearly everything above.

Finding 1 is the review. guess=read from a restricted checkpoint returns E(RKS) exactly, so the adoption path changes no number today. Findings 3, 10 and 11 are downstream detail of the same defect: the feature measures correctly, records correctly, and then applies a reference flip that the SCF quietly undoes. Fix finding 1 and the PR does what its description says.

arc/checks/spin.py has no consumers, which is why findings 8 and 9 are latent rather than live. Every trap in it is harmless today and becomes a wrong published number the day someone wires it up. Fixing the API now costs a fraction of debugging it later. The same holds for finding 6 — invalidates_analytic_freq is computed and discarded.

A smaller third chain: findings 5, 6, 7 and 11 are all the same shape — a diagnostic that is computed, recorded, and never gated on. Acting narrowly by design is defensible, but four separate measurements currently reach no decision, and finding 7 is the one most likely to matter in a real campaign.

Method, for what it's worth

Eight parallel review passes over the same pinned diff — testing/mutation, security, API-contract, maintainability, performance, data-migration, an adversarial pass, and a cross-model adversarial pass — none of which saw each other's conclusions, plus the PySCF measurement above. Findings were deduplicated by file:line and by mechanism; cross-model findings were re-derived against the source before being accepted.

Three conclusions were reached independently by two passes each, which is why I'd weight them: adoption being inert for the measured species (finding 3), the provenance riding restart.yml rather than output.yml, and the TS-switch record losing its log path (finding 11). The restart spawn guard is recorded as unresolved precisely because two passes read the same lines and disagreed.

All test gates ran in arc_env with an empty HOME and pytest -n0. Known-failing at any SHA on my box, not attributable to this branch: checks/ts_test.py::test_check_rxn_e0, ::test_compute_rxn_e0, processor_test.py::test_compare_rates, statmech/arkane_test.py::test_run_statmech_using_molecular_properties (missing arkane in rmg_env), plus molecule_test.py::TestConnectTheDots.

No ESS job was submitted at any point.

Thanks.

Please note that the the branch has moved quite ahead since this review was enacted - I believe even when reviewed it had already moved a few commits ahead. So this will require another thorough review.

• 2 - the order dependent tests. Note there were 7 lines. All seven now compare levelno, since initialize_log calls addLevelName(logging.WARNING, 'Warning: ').
• 3 - Adoption is no longer inert. The analysis spawns from spawn_post_opt_jobs rather than post_freq_actions and freq/sp/IRC are held until the verdict lands. An unstable reference also means the geometry is wrong, so an adopted verdict re-optimises rather than only re-referencing - a Hessian at the restricted geometry on an unrestricted reference sits at a non-stationary point
• Unresolved (restart guard) - resolved. The guard is now ARCSpecies.stability_analysis_ran, persisted through as_dict/from_dict

Fixed: 5, 6, 7, 8, 9, 10, 11, 12, H1, H2, H5, H6

Also, just for clarification:

• 8 - spin.py's silent s2_ls = 0 default is gone; both entry points now require a multiplicity and derive <S**2>_LS fomr it. The 16.8kcal/mol figure for a BS-doublet/HS-quartet pair reproduced exactly. The MIN_S2_SEPARATION cliff is closed by capping amplification (<S**2>_BS - <S**2>_LS)/separation rather than the denominator - that ratio is w/(1-w) in the BS determinant's high spin weight, so an idel fully spin-flipped solution sits at exactly 1 and the cap admits up to w=2/3
• 10 - fixed. The guassian reader no reports external_instability: None for an unrestricted reference matching ORCA. All four Gaussian fixtures hold exactly one <AA,BB:AA,BB> matrix and zero <AB,BA> blocks, including both UB3LYP ones, so Gaussian never tested that sector either. This changes exisitng Gaussian behaviour.