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 carries a broken-symmetry geometry and ZPE while its reactants and products do not. Yamaguchi projection gives E_projected < E_BS < E_restricted, so a broken-symmetry energy sits above the spin-pure one it approximates and the barrier the run reports is systematically overestimated by the residual contamination of the TS alone — less so than the all-restricted barrier it replaces. That direction is stated in the docstring and in advanced.rst.

Adoption governs DFT except double hybrids, and Hartree-Fock including HF-3c. A verdict measured at the optimisation level decides the reference for the geometry and the ZPE, and does not reach a correlated single point. Measured on the same TS in Molpro at cc-pVDZ:

reference	SCF	CCSD(T)-F12	T1	D1
RHF	−194.974539	−196.064366	0.0410	0.2172
broken-symmetry UHF	−195.077240	(withdrawn — see below)	—	—
triplet ROHF	−195.066148	−196.053302	0.0118	0.0369

The argument rests on the triplet ROHF row: spin relaxation drops T1 from 0.0410 to 0.0118, below the 0.015 threshold at which ARC reports a multireference species, because the orbitals absorb the static correlation the diagnostic detects — while the total energy moves under 7 kcal/mol. The diagnostic that exists to catch a multireference species stops catching it, at no cost in the energy.

A broken-symmetry coupled-cluster row was withdrawn after review, and no BS-UHF correlated energy is quoted because none exists. Molpro's ccsd(t)-f12 is the closed-shell program and uccsd(t)-f12 the ROHF-orbital open-shell one — the "U" names spin-unrestricted amplitudes, not orbitals — so neither takes a genuinely spin-broken UHF determinant as its reference, and the difference that was reported subtracted two different references' energies. The arithmetic gave it away before the mechanism did: the row implied the BS reference recovers 0.22 Eh less correlation than the triplet ROHF row, another open-shell F12 treatment of the same molecule in the same basis, while being variationally 64 kcal/mol below RHF at the SCF level. Its T1/D1 are withdrawn with it, since a diagnostic read off a calculation that is not an energy of the state is worth no more than the energy.

level_admits_a_broken_symmetry_reference() reads a DOUBLE_HYBRID_METHODS deny-list first, then gates on BROKEN_SYMMETRY_METHOD_TYPES = ['dft'] and BROKEN_SYMMETRY_METHODS = ['hf', 'hf3c', 'rhf', 'uhf', 'rohf'], rather than by exempting 'wavefunction' — hf carries that method type, and an exempt level makes job_scf_reference_is_restricted return None, which would stop the single point being recorded at all. Double hybrids are excluded because a double hybrid's energy is not its determinant's energy: it carries an MP2 term expanded about the KS determinant, which is the same construction that keeps a correlated wavefunction level on its restricted reference, and BS-MP2 about a spin-broken reference is the pathology this gate exists to prevent.

E0 therefore takes its geometry and ZPE from one reference and its electronic energy from another. check_scf_reference_consistency reports that: output.yml records reference_mismatch: true with the mixed-reference message, and T1 stays loud.

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

12 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, rebased onto current main: 3214 passed, 43 skipped, 5 failed. The 5 are arc/job/adapters/torch_ani_test.py, failing on an unconfigured TANI_PYTHON, which is None locally; 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.28%. Comparing base (90f8dff) to head (800cc8c).
⚠️ Report is 1 commits behind head on main.

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #1014      +/-   ##
==========================================
+ Coverage   64.80%   65.28%   +0.48%     
==========================================
  Files         119      120       +1     
  Lines       39997    40731     +734     
  Branches    10338    10519     +181     
==========================================
+ Hits        25920    26592     +672     
- Misses      11094    11123      +29     
- Partials     2983     3016      +33     

Flag	Coverage Δ
functionaltests	65.28% <ø> (+0.48%)	⬆️
unittests	65.28% <ø> (+0.48%)	⬆️

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.

[alongd]

Second review, at afc47f92cd0c06357fba0225e065d6b75e7ba1af against merge-base 79540b621010ed715237a8bb75b89aa897b574a6 — 42 files, +17205/−77, 19 source files across 9 packages. Same method as before: eight parallel passes over a pinned diff, none seeing another's conclusions.

The round-one blocker is fixed, and I verified it by execution rather than by reading. A GaussianAdapter built for a TS carrying an adopted verdict with checkfile=None — exactly the state adopt_stability_orbitals leaves — writes:

#P opt=(calcfc,maxcycle=100,noeigentest,ts)  guess=mix ub3lyp/def2tzvp   IOp(2/9=2000)

Both halves present: the u prefix and the symmetry-breaking directive. And the directive is what matters — PySCF B3LYP/def2-TZVP on this PR's own fixture geometry:

E(RKS)                            = -196.490056815    (fixture G16 RB3LYP: -196.490056259)
E(UKS | seeded from RKS density)  = -196.490056815    <S^2> = -0.000000   <- the round-one defect
E(UKS | guess=mix analogue)       = -196.510082249    <S^2> =  0.863368
                                                        difference: 12.5661 kcal/mol

The path the code now takes reaches the lower solution. Mutating adopt_stability_orbitals to preserve the checkfile — i.e. reintroducing the old bug — kills test_the_re_optimization_reads_the_orbitals_the_analysis_relaxed_into, so the behaviour is pinned rather than accidental.

Also fixed, each verified: the levelname→levelno test-isolation flake (162 pass with reaction_test.py first, same process); parse_wavefunction_stability now declared on the base ESSAdapter; specific_job_type: stability is no longer a no-op; both unguarded int() sites now gate on is_str_int; the three docstrings now state E_projected < E_BS < E_restricted; s2_ls now derives from the target multiplicity via target_low_spin_s_squared. Two round-one findings are refuted outright by the new code: record_scf_reference now fires for a single-level project (scheduler.py:1593), and check_spin_contamination (:2316) does compare ⟨S²⟩ against expectation.

And the ZPE splice is genuinely closed. spawn_post_opt_jobs returns at scheduler.py:2424-2429 before enqueueing freq/sp/IRC/rotors and before the if 'optfreq' in job_name freq check at :2440, so even the combined case discards its restricted Hessian; the re-optimization is job_type='opt' even when the original was optfreq, so a fresh freq is enqueued on re-entry. E_elect and ZPE both come from the re-optimized unrestricted surface.

arc/settings/submit.py is clean — all +10 lines are string literals across the six ORCA templates, no format interpolation, no globs, no widened transfer surface. The new fixtures are clean too: the ORCA ones contain no paths at all, the Gaussian ones only a generic sanitized scratch placeholder. Round one's fixture-identifier note does not recur.

That is a serious response to a hard finding. What follows is what the new surface brought with it.

---

Blockers

1. Molpro cannot honour an adopted verdict, and Gaussian-opt + Molpro-sp is the standard arrangement

arc/job/adapters/molpro.py:51-53 and :233 · confidence 9 (template read, input generation executed)

The template hardcodes the reference block:

{hf;${shift}
 maxit,999;
 wf,spin=${spin},charge=${charge};
}

${restricted}${method}

input_dict['spin'] = self.multiplicity - 1 (:229), so a singlet gives wf,spin=0 — a closed-shell RHF. :233 then prefixes u onto the method, not the reference. A real write_input_file() for an adopted singlet TS at CCSD(T)-F12/cc-pVTZ-F12 emits wf,spin=0 together with uccsd(t)-f12;. Nothing anywhere changes the reference.

Scenario: singlet H-abstraction TS, opt at ωB97X-D/def2-TZVP in Gaussian, RHF→UHF instability adopted, re-opt in Gaussian reaches the broken-symmetry geometry, sp in Molpro. Molpro either errors or returns an RHF-based CCSD(T)-F12 energy at the broken-symmetry geometry — neither the restricted answer (wrong geometry) nor the unrestricted one (wrong reference) — while output.yml records source: derived, verdict_restricted: False.

warn_on_collapsible_unrestricted_reference does fire here, which makes this a provenance defect rather than a silent one. But the run still publishes the number, and this is the common configuration, not a corner case.

The design question underneath it: is a log warning the right response when the ESS cannot honour the adoption? I'd argue no. The check is cheap and available before the sp is spawned. Refusing to adopt when the sp adapter is not in SYMMETRY_BREAKING_CAPABLE_ESS is strictly better than adopting and then reporting that the adoption didn't take — especially on a 200-species run where nobody reads the log.

2. BrokenSym 1,1 is emitted with the low-spin multiplicity

arc/job/adapters/orca.py:576-578 · confidence 8 (input generation executed; no ORCA binary available here)

The emitted input carries * xyz 0 1 together with BrokenSym 1,1. ORCA builds a broken-symmetry state by first converging the high-spin determinant of Ms1+Ms2 unpaired electrons, then localizing and flipping — this PR's own docstring at orca.py:341-343 says exactly that. At Mult 1 there are no singly-occupied orbitals to localize. ORCA's manual examples put the high-spin multiplicity on the coordinate line (BrokenSym 5,5 with xyz 0 11).

So either ORCA errors on every adopted-verdict job, or it ignores the directive and the UKS singlet collapses to the RKS solution — which is round one's defect reappearing on the fallback path. This is ARC's only symmetry-breaking mechanism for the Gaussian-opt/ORCA-sp arrangement and for a failed .gbw download.

Needs a two-minute ORCA smoke test on a cluster, which I did not run. I searched for a proven local BrokenSym/STABPerform recipe and found none, so there is no precedent to weigh against this.

3. A follow that ends still-unstable throws away the solution it found

arc/parser/adapters/orca.py:329,341-347 · confidence 9 (executed on a log assembled from the PR's own fixture blocks)

followed_to_stable=False forces s_squared_after_follow=None, which collapses a restricted instability to 'unattributed_instability', which is not adopted — so no re-optimization runs and ARC publishes the restricted energy.

ORCA allows five follow attempts (orca_rhf_uhf_instability_singlet_ts.out line 727: The wavefunction in try 1/ 5 is unstable), so "instability found, followed, still slightly unstable" is routine for a biradicaloid singlet TS. The fixture measures the cost itself: -196.344571580 (RHF, line 689) → -196.364789362 (followed, line 776) = 12.7 kcal/mol discarded.

The sector question and the "is the restricted reference the ground state?" question are not the same question, and only the second should gate adoption.

4. No TS guess after the first is ever measured

arc/scheduler.py:1866 · confidence 9 (executed grep: set at :1866, initialized species.py:411, restored :970, read :1822 — no reset anywhere in arc/)

carry_stability_verdict_across_ts_switch is otherwise careful — it clears the summary, the info fragment and four warning messages, which refutes round one's "nothing clears it" note. But stability_analysis_ran is never reset, so:

• Guess A unstable and adopted, B carried: B runs unrestricted from its first job on A's evidence. If B is a genuinely closed-shell saddle from a different search engine, its published energy is needlessly spin-contaminated and nothing can discover that. This is a regression relative to the merge-base, where B would have been computed restricted and correctly.
• Guess A stable, verdict dropped, B never measured: if B's restricted reference is unstable, ARC publishes the restricted energy with wavefunction_stability: null and no warning — indistinguishable in output.yml from "measured and stable."

The docstring at :3703 states this as intentional ("A carried verdict is therefore never contradicted by a later one"). I think that is the wrong trade: one extra stable=(rext,noopt) per switch, against an unbounded and unsignalled energy error.

5. A zero-byte .gbw is adopted as an SCF guess, defeating the scheduler's own guard

arc/job/adapters/orca.py:267-269 · confidence 9 (executed)

ssh.download_file swallows IOError and leaves a zero-byte local file. The scheduler knows this and explicitly refuses to adopt it — scheduler.py:1299-1301: "is empty, which is what a failed download leaves behind. Not adopting it…". The ORCA adapter's local-path fallback reads the same directory with os.path.isfile and no size check, re-adopting exactly the file the scheduler just rejected. reads_orbital_guess() (:290-292) also tests only isfile.

Planting a 0-byte input.gbw and rebuilding the adapter yields:

checkfile: .../opt_a1/input.gbw   size: 0
reads_orbital_guess: True
!MORead
%moinp "guess.gbw"

The 2>/dev/null in the submit templates makes the server-side cp input.gbw silent when ORCA died before writing it, so this is the quiet path. Fix: mirror the scheduler's guard with os.path.getsize(...) > 0 in both places.

(Found independently by two passes — one by execution, one by reading.)

6. Pipe-mode jobs may bypass the reference decision entirely

arc/job/pipe/pipe_planner.py, pipe_coordinator.py, arc/scheduler.py:1160 · confidence 6 — read-only, and explicitly not measured

Pipe-mode SP jobs bypass run_job(), so warn_on_collapsible_unrestricted_reference() and SCF-reference recording are skipped. With pipe_settings.enabled, a batch of adopted singlet TSs routed to an sp in a third ESS would get an unrestricted-but-spin-symmetric start, collapse to the restricted energy, have post_sp_actions() record e_elect, and emit no collapsed-reference warning and no scf_reference evidence.

warn_on_collapsible_unrestricted_reference's own docstring concedes a piped single point is not reported. arc/job/pipe/ contains no occurrence of checkfile, restricted or is_species_restricted.

I am flagging this as unmeasured: two passes reached it by reading and neither traced how a piped job composes its reference. If piped jobs do bypass the reference decision, it is a hole of the same size as round one's finding. That trace is the one I'd most want before merge.

---

Then

7. The "do not read a checkpoint" sentinel is not sticky

arc/job/adapters/gaussian.py:220-222 · confidence 9 (executed)

adopt_stability_orbitals expresses "this species must not read a checkpoint" as checkfile = None. But __init__ has a fallback that re-supplies one from whatever sits in the job's directory:

job1 local_path: .../calcs/TSs/TS/opt_a1
job2 local_path: .../calcs/TSs/TS/opt_a1
job2 .checkfile after __init__: .../calcs/TSs/TS/opt_a1/check.chk
job2 route: #P opt=(calcfc,maxcycle=100,noeigentest,ts)  guess=read ub3lyp/def2tzvp

set_file_paths (arc/job/adapter.py:394-395) reuses an existing directory and never clears it, and job_num comes from a global CSV line count that _rotate_csv_if_needed(max_lines=10000) resets to zero — so opt_a1, opt_a2 … can repeat within one project past 10,000 jobs, which a 200-species campaign reaches. It requires the same species to hold a job at both the pre- and post-rotation number, so it is not common; when it fires, the re-optimization reads a stale checkpoint under guess=read and round one's defect is back.

To be clear about what was checked and refuted: ESS troubleshooting does not restore the checkfile — scheduler.py:4697-4698 assigns species.checkfile = job.checkfile, and the re-optimization job's checkfile is None at that point. I initially believed otherwise; execution refuted it. Directory reuse is the live route, not troubleshooting.

Cheapest fix regardless: pass an explicit "do not read orbitals" flag rather than relying on None as a sentinel that three separate code paths are willing to overwrite.

8. The new readable_checkfile protocol is implemented for exactly one adapter

arc/job/adapter.py:165-166 · confidence 9 (executed)

Base defaults are check_file_name = 'check.chk', guess_file_name = 'check.chk' — Gaussian's names. Only OrcaAdapter/OrcaNebAdapter override them. Every other registered adapter inherits the Gaussian names:

JobEnum.ase/cfour/molpro/psi4/qchem/terachem/torchani/xtb ... check_file_name='check.chk'

TeraChem actually writes teracheck.chk, so readable_checkfile on a TeraChem adapter accepts a Gaussian checkfile and refuses TeraChem's own — inverted. Psi4 (commit c47d18ccc, "Refuse a checkfile written by another ESS") inherits 'check.chk' too, so it refuses ORCA's input.gbw and still accepts Gaussian's, which it also cannot read.

Not a regression — at the merge-base arc/job/adapters/common.py:166 was a bare obj.checkfile = checkfile with no guard at all, so these adapters are exactly where they were. But this PR is where the abstraction is introduced and is the natural place to finish it.

Relatedly, readable_checkfile decides by base name only, so a traversal path is accepted verbatim:

foreign /etc/passwd            -> None
traversal basename input.gbw   -> <home>/../etc/../etc/input.gbw

The "attacker" is the user's own input file and species_dict['checkfile'] was already trusted, so the real severity is low — but the docstring presents the method as the boundary that decides whether a checkfile may be read. If it is a trust boundary, add a project-directory containment check; if it is only an ESS-format sniff, say so, so no later caller leans on it.

9. arc/checks/spin.py — the arithmetic is now right, and it still has no consumers

confidence 10

The Yamaguchi implementation is correct, including the non-singlet generalization: the numerator is (s2_bs - s2_ls), verified by a derivation-independent two-state round trip (a bare ⟨S²⟩_BS numerator would be off by 0.375 Ha on a doublet). target_low_spin_s_squared returns S(S+1) for multiplicities 1-6 and refuses 0, negative and non-numeric input. 36 tests pass.

But this is the second round in which 238 lines plus a 354-line test file ship with zero callers outside arc/checks/__init__.py:3's package-listing import — and, more tellingly, nothing in ARC computes a high-spin energy at the broken-symmetry geometry, so three of the five arguments have no producer anywhere in the PR. check_spin_contamination (scheduler.py:2316) is the natural caller and doesn't use it. The interface is fixed by a docstring today and will have to change when a caller appears.

Recommendation: land spin.py + spin_test.py in the follow-up that consumes them. The rest of this PR stands without it.

Two small real defects while it is here: :158 accepts s2_bs up to 0.1 below the spin-pure target and returns an anti-projection above E_BS (executed: s2_bs=0.70 → -99.9836 against E_BS = -100.0), violating the inequality asserted at :127; and :55,80 accept half-integer multiplicities (1.5 → 0.3125) that name no spin state, contrary to :62.

10. Quadratic backtracking in one new regex

arc/parser/adapters/orca.py:315 · confidence 9 (executed)

re.match(r'\s*E\(\s*\d+\)\s*=\s*([-+]?\d*\.?\d+(?:[EeDd][-+]?\d+)?)\s*Eh', line)

\d*\.?\d+ over a common character class followed by a mandatory Eh that can fail is the classic ambiguous-quantifier shape. Measured clean O(N²):

2000  0.055
4000  0.221
8000  0.87
16000 3.671

~370 s at 160k characters on one line, and _get_lines_from_file has no line-length bound. check_stability_job catches exceptions but cannot catch a hang, and the scheduler loop is single-threaded.

Reachability is low — ORCA never writes that shape, so it needs someone with write access to the log ARC parses. Worth noting it is the only one of the new regexes with the defect: the two new Gaussian ones measured 0.000–0.001 s at 32k chars.

11. The reference flip is species-global, never scoped to the level that measured it

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

REFERENCE_AGNOSTIC_METHOD_TYPES exempts force field, composite and semiempirical; wavefunction is not exempt, so a DFT verdict rewrites a CCSD(T)-F12 sp. For a mildly RHF-unstable TS, RHF-CCSD(T)-F12 is frequently the better number — CCSD(T) recovers the static correlation the DFT reference could not, while a broken-symmetry UHF-CCSD(T) inherits the contamination. There is no knob to decline per level.

This is a policy question rather than a defect, and it is the maintainers' to decide. I would scope adoption to method_type == 'dft' and make the correlated-level behaviour explicit.

12. Smaller

• arc/scheduler.py:2171 — invalidates_analytic_freq is consumed only into a warnings string. The two verdicts that set it True (internal_instability, and external_instability on an unrestricted reference) are precisely the two ARC never adopts, so the freq is released and the analytic Hessian runs on a wavefunction the code just declared outside its domain of validity.
• arc/scheduler.py:2123 — an ORCA stability job that dies after printing the wavefunction is unstable is gated out by job_status != 'done' and never re-run (:1232), so a readable verdict is discarded and the species proceeds restricted with nothing in output.yml.
• arc/parser/adapters/orca.py:329 — followed_to_stable = len(blocks) > 1 and blocks[-1]['verdict'] == 'stable' never checks that blocks[0] was unstable, so two concatenated stable analyses yield True with an ⟨S²⟩ from a wavefunction the log never relaxed into. Unreachable through ARC's own composition, so a contract defect on public parser API rather than a live path. The log carries a direct anchor: Restart requested: Orbitals will be transformed, printed once per follow.
• arc/scheduler.py:1974 — the followed_to_stable guard in adopt_stability_orbitals is circular for ORCA: it is only reached when the verdict is external_instability, which orca.py:341-342 can only produce when followed_to_stable was already True. Worth saying out loud so nobody reads the two as independent confirmations.
• arc/species/species.py:214-218 vs :967-968 — the docstring says derived_stability_verdict "is never inputted by the user", but from_dict reads it unconditionally with no shape validation. A hand-written verdict in a restart file silently flips every job to unrestricted with no measurement behind it. Same channel as any other restart field, so low severity, but the docstring's claim is not enforced.
• arc/scheduler.py:1593 — round one's sp_level == opt_level blindness is refuted, but the restart sub-branch that falls back to output[label]['paths']['geo'] passes no job, so reference_mismatch stays null for a restarted single-level project.
• check_spin_contamination compares and warns, but the comparison changes no decision, and wavefunction_stability: stable with a 1.0 deviation land in the same output.yml block with no cross-reference.
• Essay docstrings. Several new docstrings run 40+ lines with ad-hoc ALL-CAPS section headers (TWO MECHANISMS PREVENT THAT, WHAT IS REPORTED WHERE). The content is genuinely load-bearing physics that isn't recoverable from the code — but it names SYMMETRY_BREAKING_CAPABLE_ESS and STABILITY_CAPABLE_ESS in prose, with nothing keeping that in sync when a third ESS is added. Consider moving the "why" into advanced.rst, which this PR already extends by +167 lines.

---

Cost

Not a defect — the re-optimization gate is sound. stability_reoptimized is set before run_job, persisted, and reset by nothing: not delete_all_species_jobs, not switch_ts, not the restart load. Three passes independently failed to construct a fail-open path, and the torn-write window (a hard kill between the verdict being persisted at :2178 and the flag being set) is closed by release_held_stability_work() re-entering exactly once on the next tick.

But the delta versus round one is large and worth stating plainly: one adoption is not one extra job. The re-optimization re-enters spawn_post_opt_jobs, which re-queues freq, sp, every rotor scan, and 2×IRC for a TS on the new geometry — necessary, since the old stationary point is gone. On top of that, every Gaussian-adopted species has checkfile = None for the duration, so the re-opt and its whole cascade start SCF cold rather than from guess=read. The two compound.

The qualifying population is not narrow either: every TS, plus every restricted-reference DFT/HF opt on a stability-capable ESS.

---

Checked and sound — please don't re-verify

• Round-one P0 fixed, by execution (route line + PySCF energies above), and pinned by a test.
• ZPE splice closed — E_elect and ZPE both come from the re-optimized surface.
• Re-optimization gate does not fail open — three passes, three routes, same conclusion.
• All five new ARCSpecies attributes (derived_stability_verdict, scf_references, stability_analysis_ran, stability_pending_opt_job, stability_reoptimized) round-trip correctly through as_dict() → YAML → species_dict=. No written-but-unread or read-but-unwritten field.
• checkfile = None survives a restart — as_dict() omits the key and the constructor defaults it, so no stale resurrection.
• Genuine pre-PR restart fixtures load clean — all 4 species of 1_restart_thermo/restart.yml, every new attribute at its documented default, no respawn.
• arc/settings/submit.py clean — no interpolation, no globs, no widened transfer. All six ORCA templates covered, none skipped.
• New fixtures carry no account, hostname, node name, scratch layout, licence number or credential.
• orca_rhf_uhf_instability_no_restart_crash.out is substantively exercised, not shape-only — tests assert verdict, restricted, followed_to_stable, exact lowest_eigenvalue and specific Nones.
• parse_wavefunction_stability now declared on the base; all 9 parser subclasses enumerated, the 7 without an override fall through to None correctly.
• Both round-one unguarded int() sites fixed; the new code does not repeat the pattern.
• A foreign ORCA .gbw handed to the Gaussian adapter is refused and falls through to guess=mix — the cross-ESS guess=read crash I hypothesised does not exist.
• !MORead and BrokenSym cannot both fire — if/else, and the predicate independently returns None when a guess is read.
• Parser exceptions cannot kill the scheduler — both consumers catch bare Exception.
• ORCA download scope is correct — freq/sp/scan only upload a guess and never pull back a .gbw they can't produce.
• Mutants killed across the ORCA followed_to_stable boundaries, the BrokenSym electron-count gate, _is_finite_s2's zero boundary, s_squared_expected_from_multiplicity's singlet boundary, and the stability_reoptimized guard.

Surviving mutants (coverage gaps, not live bugs): spin.py:148 <→<= and :158 <→<= both leave 36/36 passing; orca.py multiplicity != 1→> 1 leaves 26/26 (behaviourally equivalent on reachable input).

---

What I could not reach

• Whether guess=mix and a followed instability converge to the same broken-symmetry solution. The PR's own orca.py:345-346 says "the two can converge to different broken-symmetry solutions." I tried to settle it numerically and failed — my PySCF reconstruction of the followed orbitals collapsed back to the restricted solution, so that arm is unfaithful. What I can say is that guess=mix reaches a solution 12.57 kcal/mol below restricted. Whether Gaussian and ORCA publish the same energy for the same TS at the same level is open, and for a change to a published physical quantity that feels like a measurement the PR should carry rather than a docstring caveat.
• Whether ORCA accepts BrokenSym 1,1 at Mult 1 (blocker 2) — no ORCA binary available; needs a cluster smoke test.
• Whether stable=(rext,noopt) can emit RHF -> CRHF — if it can, adopted_reference_is_unrestricted has no relaxation-class gate, while the ORCA adapter has exactly that guard (derived_instability_breaks_spin_symmetry, orca.py:390) wired into the ESS that does not decide adoption. One or clause if it turns out to be reachable.
• The pipe (blocker 6) — the trace I'd most want before merge.
• No ESS job was submitted at any point in this review.

[calvinp0]

Second review, at afc47f92cd0c06357fba0225e065d6b75e7ba1af against merge-base 79540b621010ed715237a8bb75b89aa897b574a6 — 42 files, +17205/−77, 19 source files across 9 packages. Same method as before: eight parallel passes over a pinned diff, none seeing another's conclusions.

The round-one blocker is fixed, and I verified it by execution rather than by reading. A GaussianAdapter built for a TS carrying an adopted verdict with checkfile=None — exactly the state adopt_stability_orbitals leaves — writes:

#P opt=(calcfc,maxcycle=100,noeigentest,ts)  guess=mix ub3lyp/def2tzvp   IOp(2/9=2000)

Both halves present: the u prefix and the symmetry-breaking directive. And the directive is what matters — PySCF B3LYP/def2-TZVP on this PR's own fixture geometry:

E(RKS)                            = -196.490056815    (fixture G16 RB3LYP: -196.490056259)
E(UKS | seeded from RKS density)  = -196.490056815    <S^2> = -0.000000   <- the round-one defect
E(UKS | guess=mix analogue)       = -196.510082249    <S^2> =  0.863368
                                                        difference: 12.5661 kcal/mol

The path the code now takes reaches the lower solution. Mutating adopt_stability_orbitals to preserve the checkfile — i.e. reintroducing the old bug — kills test_the_re_optimization_reads_the_orbitals_the_analysis_relaxed_into, so the behaviour is pinned rather than accidental.

Also fixed, each verified: the levelname→levelno test-isolation flake (162 pass with reaction_test.py first, same process); parse_wavefunction_stability now declared on the base ESSAdapter; specific_job_type: stability is no longer a no-op; both unguarded int() sites now gate on is_str_int; the three docstrings now state E_projected < E_BS < E_restricted; s2_ls now derives from the target multiplicity via target_low_spin_s_squared. Two round-one findings are refuted outright by the new code: record_scf_reference now fires for a single-level project (scheduler.py:1593), and check_spin_contamination (:2316) does compare ⟨S²⟩ against expectation.

And the ZPE splice is genuinely closed. spawn_post_opt_jobs returns at scheduler.py:2424-2429 before enqueueing freq/sp/IRC/rotors and before the if 'optfreq' in job_name freq check at :2440, so even the combined case discards its restricted Hessian; the re-optimization is job_type='opt' even when the original was optfreq, so a fresh freq is enqueued on re-entry. E_elect and ZPE both come from the re-optimized unrestricted surface.

arc/settings/submit.py is clean — all +10 lines are string literals across the six ORCA templates, no format interpolation, no globs, no widened transfer surface. The new fixtures are clean too: the ORCA ones contain no paths at all, the Gaussian ones only a generic sanitized scratch placeholder. Round one's fixture-identifier note does not recur.

That is a serious response to a hard finding. What follows is what the new surface brought with it.

Blockers

1. Molpro cannot honour an adopted verdict, and Gaussian-opt + Molpro-sp is the standard arrangement

arc/job/adapters/molpro.py:51-53 and :233 · confidence 9 (template read, input generation executed)

The template hardcodes the reference block:

{hf;${shift}
 maxit,999;
 wf,spin=${spin},charge=${charge};
}

${restricted}${method}

input_dict['spin'] = self.multiplicity - 1 (:229), so a singlet gives wf,spin=0 — a closed-shell RHF. :233 then prefixes u onto the method, not the reference. A real write_input_file() for an adopted singlet TS at CCSD(T)-F12/cc-pVTZ-F12 emits wf,spin=0 together with uccsd(t)-f12;. Nothing anywhere changes the reference.

Scenario: singlet H-abstraction TS, opt at ωB97X-D/def2-TZVP in Gaussian, RHF→UHF instability adopted, re-opt in Gaussian reaches the broken-symmetry geometry, sp in Molpro. Molpro either errors or returns an RHF-based CCSD(T)-F12 energy at the broken-symmetry geometry — neither the restricted answer (wrong geometry) nor the unrestricted one (wrong reference) — while output.yml records source: derived, verdict_restricted: False.

warn_on_collapsible_unrestricted_reference does fire here, which makes this a provenance defect rather than a silent one. But the run still publishes the number, and this is the common configuration, not a corner case.

The design question underneath it: is a log warning the right response when the ESS cannot honour the adoption? I'd argue no. The check is cheap and available before the sp is spawned. Refusing to adopt when the sp adapter is not in SYMMETRY_BREAKING_CAPABLE_ESS is strictly better than adopting and then reporting that the adoption didn't take — especially on a 200-species run where nobody reads the log.

2. BrokenSym 1,1 is emitted with the low-spin multiplicity

arc/job/adapters/orca.py:576-578 · confidence 8 (input generation executed; no ORCA binary available here)

The emitted input carries * xyz 0 1 together with BrokenSym 1,1. ORCA builds a broken-symmetry state by first converging the high-spin determinant of Ms1+Ms2 unpaired electrons, then localizing and flipping — this PR's own docstring at orca.py:341-343 says exactly that. At Mult 1 there are no singly-occupied orbitals to localize. ORCA's manual examples put the high-spin multiplicity on the coordinate line (BrokenSym 5,5 with xyz 0 11).

So either ORCA errors on every adopted-verdict job, or it ignores the directive and the UKS singlet collapses to the RKS solution — which is round one's defect reappearing on the fallback path. This is ARC's only symmetry-breaking mechanism for the Gaussian-opt/ORCA-sp arrangement and for a failed .gbw download.

Needs a two-minute ORCA smoke test on a cluster, which I did not run. I searched for a proven local BrokenSym/STABPerform recipe and found none, so there is no precedent to weigh against this.

3. A follow that ends still-unstable throws away the solution it found

arc/parser/adapters/orca.py:329,341-347 · confidence 9 (executed on a log assembled from the PR's own fixture blocks)

followed_to_stable=False forces s_squared_after_follow=None, which collapses a restricted instability to 'unattributed_instability', which is not adopted — so no re-optimization runs and ARC publishes the restricted energy.

ORCA allows five follow attempts (orca_rhf_uhf_instability_singlet_ts.out line 727: The wavefunction in try 1/ 5 is unstable), so "instability found, followed, still slightly unstable" is routine for a biradicaloid singlet TS. The fixture measures the cost itself: -196.344571580 (RHF, line 689) → -196.364789362 (followed, line 776) = 12.7 kcal/mol discarded.

The sector question and the "is the restricted reference the ground state?" question are not the same question, and only the second should gate adoption.

4. No TS guess after the first is ever measured

arc/scheduler.py:1866 · confidence 9 (executed grep: set at :1866, initialized species.py:411, restored :970, read :1822 — no reset anywhere in arc/)

carry_stability_verdict_across_ts_switch is otherwise careful — it clears the summary, the info fragment and four warning messages, which refutes round one's "nothing clears it" note. But stability_analysis_ran is never reset, so:

* **Guess A unstable and adopted, B carried:** B runs unrestricted from its first job on A's evidence. If B is a genuinely closed-shell saddle from a different search engine, its published energy is needlessly spin-contaminated and nothing can discover that. This is a regression relative to the merge-base, where B would have been computed restricted and correctly.

* **Guess A stable, verdict dropped, B never measured:** if B's restricted reference _is_ unstable, ARC publishes the restricted energy with `wavefunction_stability: null` and no warning — indistinguishable in `output.yml` from "measured and stable."

The docstring at :3703 states this as intentional ("A carried verdict is therefore never contradicted by a later one"). I think that is the wrong trade: one extra stable=(rext,noopt) per switch, against an unbounded and unsignalled energy error.

5. A zero-byte .gbw is adopted as an SCF guess, defeating the scheduler's own guard

arc/job/adapters/orca.py:267-269 · confidence 9 (executed)

ssh.download_file swallows IOError and leaves a zero-byte local file. The scheduler knows this and explicitly refuses to adopt it — scheduler.py:1299-1301: "is empty, which is what a failed download leaves behind. Not adopting it…". The ORCA adapter's local-path fallback reads the same directory with os.path.isfile and no size check, re-adopting exactly the file the scheduler just rejected. reads_orbital_guess() (:290-292) also tests only isfile.

Planting a 0-byte input.gbw and rebuilding the adapter yields:

checkfile: .../opt_a1/input.gbw   size: 0
reads_orbital_guess: True
!MORead
%moinp "guess.gbw"

The 2>/dev/null in the submit templates makes the server-side cp input.gbw silent when ORCA died before writing it, so this is the quiet path. Fix: mirror the scheduler's guard with os.path.getsize(...) > 0 in both places.

(Found independently by two passes — one by execution, one by reading.)

6. Pipe-mode jobs may bypass the reference decision entirely

arc/job/pipe/pipe_planner.py, pipe_coordinator.py, arc/scheduler.py:1160 · confidence 6 — read-only, and explicitly not measured

Pipe-mode SP jobs bypass run_job(), so warn_on_collapsible_unrestricted_reference() and SCF-reference recording are skipped. With pipe_settings.enabled, a batch of adopted singlet TSs routed to an sp in a third ESS would get an unrestricted-but-spin-symmetric start, collapse to the restricted energy, have post_sp_actions() record e_elect, and emit no collapsed-reference warning and no scf_reference evidence.

warn_on_collapsible_unrestricted_reference's own docstring concedes a piped single point is not reported. arc/job/pipe/ contains no occurrence of checkfile, restricted or is_species_restricted.

I am flagging this as unmeasured: two passes reached it by reading and neither traced how a piped job composes its reference. If piped jobs do bypass the reference decision, it is a hole of the same size as round one's finding. That trace is the one I'd most want before merge.

Then

7. The "do not read a checkpoint" sentinel is not sticky

arc/job/adapters/gaussian.py:220-222 · confidence 9 (executed)

adopt_stability_orbitals expresses "this species must not read a checkpoint" as checkfile = None. But __init__ has a fallback that re-supplies one from whatever sits in the job's directory:

job1 local_path: .../calcs/TSs/TS/opt_a1
job2 local_path: .../calcs/TSs/TS/opt_a1
job2 .checkfile after __init__: .../calcs/TSs/TS/opt_a1/check.chk
job2 route: #P opt=(calcfc,maxcycle=100,noeigentest,ts)  guess=read ub3lyp/def2tzvp

set_file_paths (arc/job/adapter.py:394-395) reuses an existing directory and never clears it, and job_num comes from a global CSV line count that _rotate_csv_if_needed(max_lines=10000) resets to zero — so opt_a1, opt_a2 … can repeat within one project past 10,000 jobs, which a 200-species campaign reaches. It requires the same species to hold a job at both the pre- and post-rotation number, so it is not common; when it fires, the re-optimization reads a stale checkpoint under guess=read and round one's defect is back.

To be clear about what was checked and refuted: ESS troubleshooting does not restore the checkfile — scheduler.py:4697-4698 assigns species.checkfile = job.checkfile, and the re-optimization job's checkfile is None at that point. I initially believed otherwise; execution refuted it. Directory reuse is the live route, not troubleshooting.

Cheapest fix regardless: pass an explicit "do not read orbitals" flag rather than relying on None as a sentinel that three separate code paths are willing to overwrite.

8. The new readable_checkfile protocol is implemented for exactly one adapter

arc/job/adapter.py:165-166 · confidence 9 (executed)

Base defaults are check_file_name = 'check.chk', guess_file_name = 'check.chk' — Gaussian's names. Only OrcaAdapter/OrcaNebAdapter override them. Every other registered adapter inherits the Gaussian names:

JobEnum.ase/cfour/molpro/psi4/qchem/terachem/torchani/xtb ... check_file_name='check.chk'

TeraChem actually writes teracheck.chk, so readable_checkfile on a TeraChem adapter accepts a Gaussian checkfile and refuses TeraChem's own — inverted. Psi4 (commit c47d18ccc, "Refuse a checkfile written by another ESS") inherits 'check.chk' too, so it refuses ORCA's input.gbw and still accepts Gaussian's, which it also cannot read.

Not a regression — at the merge-base arc/job/adapters/common.py:166 was a bare obj.checkfile = checkfile with no guard at all, so these adapters are exactly where they were. But this PR is where the abstraction is introduced and is the natural place to finish it.

Relatedly, readable_checkfile decides by base name only, so a traversal path is accepted verbatim:

foreign /etc/passwd            -> None
traversal basename input.gbw   -> <home>/../etc/../etc/input.gbw

The "attacker" is the user's own input file and species_dict['checkfile'] was already trusted, so the real severity is low — but the docstring presents the method as the boundary that decides whether a checkfile may be read. If it is a trust boundary, add a project-directory containment check; if it is only an ESS-format sniff, say so, so no later caller leans on it.

9. arc/checks/spin.py — the arithmetic is now right, and it still has no consumers

confidence 10

The Yamaguchi implementation is correct, including the non-singlet generalization: the numerator is (s2_bs - s2_ls), verified by a derivation-independent two-state round trip (a bare ⟨S²⟩_BS numerator would be off by 0.375 Ha on a doublet). target_low_spin_s_squared returns S(S+1) for multiplicities 1-6 and refuses 0, negative and non-numeric input. 36 tests pass.

But this is the second round in which 238 lines plus a 354-line test file ship with zero callers outside arc/checks/__init__.py:3's package-listing import — and, more tellingly, nothing in ARC computes a high-spin energy at the broken-symmetry geometry, so three of the five arguments have no producer anywhere in the PR. check_spin_contamination (scheduler.py:2316) is the natural caller and doesn't use it. The interface is fixed by a docstring today and will have to change when a caller appears.

Recommendation: land spin.py + spin_test.py in the follow-up that consumes them. The rest of this PR stands without it.

Two small real defects while it is here: :158 accepts s2_bs up to 0.1 below the spin-pure target and returns an anti-projection above E_BS (executed: s2_bs=0.70 → -99.9836 against E_BS = -100.0), violating the inequality asserted at :127; and :55,80 accept half-integer multiplicities (1.5 → 0.3125) that name no spin state, contrary to :62.

10. Quadratic backtracking in one new regex

arc/parser/adapters/orca.py:315 · confidence 9 (executed)

re.match(r'\s*E\(\s*\d+\)\s*=\s*([-+]?\d*\.?\d+(?:[EeDd][-+]?\d+)?)\s*Eh', line)

\d*\.?\d+ over a common character class followed by a mandatory Eh that can fail is the classic ambiguous-quantifier shape. Measured clean O(N²):

2000  0.055
4000  0.221
8000  0.87
16000 3.671

~370 s at 160k characters on one line, and _get_lines_from_file has no line-length bound. check_stability_job catches exceptions but cannot catch a hang, and the scheduler loop is single-threaded.

Reachability is low — ORCA never writes that shape, so it needs someone with write access to the log ARC parses. Worth noting it is the only one of the new regexes with the defect: the two new Gaussian ones measured 0.000–0.001 s at 32k chars.

11. The reference flip is species-global, never scoped to the level that measured it

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

REFERENCE_AGNOSTIC_METHOD_TYPES exempts force field, composite and semiempirical; wavefunction is not exempt, so a DFT verdict rewrites a CCSD(T)-F12 sp. For a mildly RHF-unstable TS, RHF-CCSD(T)-F12 is frequently the better number — CCSD(T) recovers the static correlation the DFT reference could not, while a broken-symmetry UHF-CCSD(T) inherits the contamination. There is no knob to decline per level.

This is a policy question rather than a defect, and it is the maintainers' to decide. I would scope adoption to method_type == 'dft' and make the correlated-level behaviour explicit.

12. Smaller

* **`arc/scheduler.py:2171`** — `invalidates_analytic_freq` is consumed only into a warnings string. The two verdicts that set it True (`internal_instability`, and `external_instability` on an unrestricted reference) are precisely the two ARC never adopts, so the freq is released and the analytic Hessian runs on a wavefunction the code just declared outside its domain of validity.

* **`arc/scheduler.py:2123`** — an ORCA stability job that dies after printing `the wavefunction is unstable` is gated out by `job_status != 'done'` and never re-run (`:1232`), so a readable verdict is discarded and the species proceeds restricted with nothing in `output.yml`.

* **`arc/parser/adapters/orca.py:329`** — `followed_to_stable = len(blocks) > 1 and blocks[-1]['verdict'] == 'stable'` never checks that `blocks[0]` was unstable, so two concatenated `stable` analyses yield `True` with an ⟨S²⟩ from a wavefunction the log never relaxed into. Unreachable through ARC's own composition, so a contract defect on public parser API rather than a live path. The log carries a direct anchor: `Restart requested: Orbitals will be transformed`, printed once per follow.

* **`arc/scheduler.py:1974`** — the `followed_to_stable` guard in `adopt_stability_orbitals` is circular for ORCA: it is only reached when the verdict is `external_instability`, which `orca.py:341-342` can only produce when `followed_to_stable` was already True. Worth saying out loud so nobody reads the two as independent confirmations.

* **`arc/species/species.py:214-218` vs `:967-968`** — the docstring says `derived_stability_verdict` "is never inputted by the user", but `from_dict` reads it unconditionally with no shape validation. A hand-written verdict in a restart file silently flips every job to unrestricted with no measurement behind it. Same channel as any other restart field, so low severity, but the docstring's claim is not enforced.

* **`arc/scheduler.py:1593`** — round one's `sp_level == opt_level` blindness is refuted, but the restart sub-branch that falls back to `output[label]['paths']['geo']` passes no `job`, so `reference_mismatch` stays `null` for a _restarted_ single-level project.

* **`check_spin_contamination`** compares and warns, but the comparison changes no decision, and `wavefunction_stability: stable` with a 1.0 deviation land in the same `output.yml` block with no cross-reference.

* **Essay docstrings.** Several new docstrings run 40+ lines with ad-hoc ALL-CAPS section headers (`TWO MECHANISMS PREVENT THAT`, `WHAT IS REPORTED WHERE`). The content is genuinely load-bearing physics that isn't recoverable from the code — but it names `SYMMETRY_BREAKING_CAPABLE_ESS` and `STABILITY_CAPABLE_ESS` in prose, with nothing keeping that in sync when a third ESS is added. Consider moving the "why" into `advanced.rst`, which this PR already extends by +167 lines.

Cost

Not a defect — the re-optimization gate is sound. stability_reoptimized is set before run_job, persisted, and reset by nothing: not delete_all_species_jobs, not switch_ts, not the restart load. Three passes independently failed to construct a fail-open path, and the torn-write window (a hard kill between the verdict being persisted at :2178 and the flag being set) is closed by release_held_stability_work() re-entering exactly once on the next tick.

But the delta versus round one is large and worth stating plainly: one adoption is not one extra job. The re-optimization re-enters spawn_post_opt_jobs, which re-queues freq, sp, every rotor scan, and 2×IRC for a TS on the new geometry — necessary, since the old stationary point is gone. On top of that, every Gaussian-adopted species has checkfile = None for the duration, so the re-opt and its whole cascade start SCF cold rather than from guess=read. The two compound.

The qualifying population is not narrow either: every TS, plus every restricted-reference DFT/HF opt on a stability-capable ESS.

Checked and sound — please don't re-verify

* **Round-one P0 fixed**, by execution (route line + PySCF energies above), and pinned by a test.

* **ZPE splice closed** — E_elect and ZPE both come from the re-optimized surface.

* **Re-optimization gate does not fail open** — three passes, three routes, same conclusion.

* **All five new `ARCSpecies` attributes** (`derived_stability_verdict`, `scf_references`, `stability_analysis_ran`, `stability_pending_opt_job`, `stability_reoptimized`) round-trip correctly through `as_dict()` → YAML → `species_dict=`. No written-but-unread or read-but-unwritten field.

* **`checkfile = None` survives a restart** — `as_dict()` omits the key and the constructor defaults it, so no stale resurrection.

* **Genuine pre-PR restart fixtures load clean** — all 4 species of `1_restart_thermo/restart.yml`, every new attribute at its documented default, no respawn.

* **`arc/settings/submit.py` clean** — no interpolation, no globs, no widened transfer. All six ORCA templates covered, none skipped.

* **New fixtures carry no account, hostname, node name, scratch layout, licence number or credential.**

* **`orca_rhf_uhf_instability_no_restart_crash.out` is substantively exercised**, not shape-only — tests assert `verdict`, `restricted`, `followed_to_stable`, exact `lowest_eigenvalue` and specific `None`s.

* **`parse_wavefunction_stability` now declared on the base**; all 9 parser subclasses enumerated, the 7 without an override fall through to `None` correctly.

* **Both round-one unguarded `int()` sites fixed**; the new code does not repeat the pattern.

* **A foreign ORCA `.gbw` handed to the Gaussian adapter is refused** and falls through to `guess=mix` — the cross-ESS `guess=read` crash I hypothesised does not exist.

* **`!MORead` and `BrokenSym` cannot both fire** — if/else, and the predicate independently returns `None` when a guess is read.

* **Parser exceptions cannot kill the scheduler** — both consumers catch bare `Exception`.

* **ORCA download scope is correct** — freq/sp/scan only upload a guess and never pull back a `.gbw` they can't produce.

* **Mutants killed** across the ORCA `followed_to_stable` boundaries, the `BrokenSym` electron-count gate, `_is_finite_s2`'s zero boundary, `s_squared_expected_from_multiplicity`'s singlet boundary, and the `stability_reoptimized` guard.

Surviving mutants (coverage gaps, not live bugs): spin.py:148 <→<= and :158 <→<= both leave 36/36 passing; orca.py multiplicity != 1→> 1 leaves 26/26 (behaviourally equivalent on reachable input).

What I could not reach

* **Whether `guess=mix` and a followed instability converge to the same broken-symmetry solution.** The PR's own `orca.py:345-346` says "the two can converge to different broken-symmetry solutions." I tried to settle it numerically and failed — my PySCF reconstruction of the followed orbitals collapsed back to the restricted solution, so that arm is unfaithful. What I can say is that `guess=mix` reaches _a_ solution 12.57 kcal/mol below restricted. Whether Gaussian and ORCA publish the same energy for the same TS at the same level is open, and for a change to a published physical quantity that feels like a measurement the PR should carry rather than a docstring caveat.

* **Whether ORCA accepts `BrokenSym 1,1` at `Mult 1`** (blocker 2) — no ORCA binary available; needs a cluster smoke test.

* **Whether `stable=(rext,noopt)` can emit `RHF -> CRHF`** — if it can, `adopted_reference_is_unrestricted` has no relaxation-class gate, while the ORCA adapter has exactly that guard (`derived_instability_breaks_spin_symmetry`, `orca.py:390`) wired into the ESS that does _not_ decide adoption. One `or` clause if it turns out to be reachable.

* **The pipe** (blocker 6) — the trace I'd most want before merge.

* No ESS job was submitted at any point in this review.

2. BrokenSym 1,1 with * xyz 0 1 - So this one I disagree with. I ran the job with ** ORCA 6.0.0, !UKS B3LYP def2-TZVP TightSCF defgrid3, * xyz 0 1, %scf ... BrokenSym 1,1 end. ORCA builds the high-spin determinant internally from the Na,Nb operans - the log prints High-Spin SCF calculation and S(High-Spin) = 1.0 - then flips. <S2> runs 2.009324 → 0.864748 → 0.864748, energy -196.364789363862, matching to 1e-9 Eh both the !MORead-from-broken-symmetry-orbitals run and the solution ORCA reaches by following the instability. The coordinate line keeps the target multiplicity; the operancs supply the high-spin state. No change made.

Also, emission is hard gated on multiplicity != 1, but charge is not - closed shell cation with an even electron count would get * xyz +1 1 with BrokenSym 1,1. Untested.

**1. Molpro - True, it is a problem. So, Molpro is not the limitation. {uhf; wf, spin=0} is valid Molpro and ROTATE, i.sym, j.sym, 45 is how one requests a broken-symmetry singled. In this instance, ARC's adapter is the limitation. It hardcodes (molpro.py:50) {hf in the template, one occurrence, never varying — and byte-identical on main, so this PR didn't introduce it. Note also {hf} isn't RHF: arc/testing/trsh/molpro/insufficient_memory.out:199 shows {hf; wf,spin=1} → PROGRAM * RHF-SCF (OPEN SHELL), i.e. ROHF, and the u at :233 selects Molpro's UCCSD(T), which is the ROHF-orbital CC. A blanket {hf}→{uhf} swap would move every radical in every ARC run off ROHF.

So I measured the alternativ propelry. On the C5H10 singlet TS at cc-pVDZ

{uhf; wf,spin=0} -194.974539 <S2> 0.000000 collapses to RHF
{uhf; wf,spin=2} then {uhf; spin=0; start,} -194.974539 <S2> 0.000000 collapses to RHF
{uhf; rotate,20.1,21.1,45} -195.077240 <S**2> 1.006489 breaks symmetry

and on 90°-twisted ethylene (D2d, deliberately symmetric) ROTATE with the naive n/2 index fails until symmetry,nosym is added, after which it works. So the recipe exists: nosym + {uhf} + rotate,{n//2}.1,{n//2+1}.1,45, scoped to adopted species, using count_electrons.

Then, I ran it through the singlet point (proving why the gate stays - will come up with a future stack PR for it):

reference SCF CCSD(T)-F12 T1 D1
RHF -194.974539 -196.064366 0.0410 0.2172
broken-symmetry UHF -195.077240 -195.840652 0.0124 0.0594
triplet ROHF -195.066148 -196.053302 0.0118 0.0369

The broken-symmetry reference is 64 kcal/mol lower at SCF and 140 kcal/mol higher after coupled cluster, recovering 0.35 Eh less correlation. You cannot substitute one for the other in a barrier.

And the second row is the one that decided it: T1 drops from 0.041 to 0.0124 - under ARC's 0.015 warning threshold. The orbitals absorb the static correlation, so the diagnostic that exists to detect a multireference species stops detecting it. ARC already flags this TS correctly today; adopting would have replaced a flagged wrong number with a quietly-unflagged one.

So, the gate is Scheduler.stability_verdict_can_be_honoured, evaluated when the verdict is recorded (before adoption takes effect), testing only the three levels the E0 is built from and skipping reference-agnostic ones. The verdict is still measured and recorded for every species - only the adoption is refused, and only for a verdict ARC would act on. All-Gaussian and all-ORCA are provably unchanged.

True regarding STABILITY_CAPABLE_ESS/SYMMETRY_BREAKING_CAPABLE_ESS said "ESS that cannot" when they meant "adapter ARC has not taught". They're now STABILITY_ANALYSIS_ADAPTERS/SYMMETRY_BREAKING_ADAPTERS, documented at the definition, and the warning now names the physics: an open-shell singlet is a two-determinant state, a broken-symmetry reference approximates rather than describes it, and a multireference treatment (CASSCF → MRCI/CASPT2) is what such a species calls for.

3. Regarding ending unstable - yes, agreed and fixed. The sector is read off any solition the log relacxed into, not only a stable one; ORCA re-converges before each analysis, so the <S2> it stoppped on belongs to a converged determinant whichever try that was. A single-analysis log still measures nothing and stays unattributed_inability. THe item-12 followed_to_unstable contract defect is fixed in the same expression - it now requires the first block to have been unstable, so concatenated stable analyses aren't read as a follow. Declined the Restart requested: anchor as version-fragile.

**4. stability_analysis_ran - agreed in part. Cleared when the verdict is dropped, kept when it's carried. Your second failure mode is real and the reset closes it. The first one though, `guess=mix and BrokenSym are initial guesses so on a genuinely closed shell saffle, the unrestricted SCF relaces back and E(UKS) = E(RKS). Measuring a carried verdict would also let a stable verdict on the unrestricted reference replace the adopted one mid-species, flipping the reference under queued jobs.

**5. Zero byte .gbw - agreed and fixed. The guard went into JobAdapter.readable_checkfile, which now refuses a path naming no file and a file of size zero - covering _initialize_adapter, both adapters' species-checkfile fallbacks, Psi4 and TeraChem in one place. Both adapters job-directory fallbacks were bypassing it entirely and now go through it.

**6. The pipe - traced and it does not bypass the reference deciison. _pending_pipe_sp/_freq/_irc are populated only inside spawn_post_opt_jobs, after the stability gate returns, so no piped job is composed before the verdict is in. build_species_leaf_task puts species.as_dict() in the payload, as_dict() emits derived_stability_verdict, and pipe_worker._run_adapter rebuilds the species and calls job_factory. Measured end to end on a real payload: the verdict survives the round trip, restricted_used=False, and both guess=mix and BrokenSym 1,1 are present. What the pipe genuinely skips is provenance - warn_on_collapsible_unrestricted_reference is called only from run_job, and _ingest_species_sp writes e_elect without post_sp_actions, so no record_scf_reference and no spin-contamination check. Both documented in advanced.rst

7-12 (excl. 11). Sentinel made sticky via species_may_read_previous_orbitals rather than a new kwarg through every adapter signature. TeraChem's check_file_name declared and routed through readable_checkfile; Psi4 left, it reads no guess. spin.py's amplification clamps at zero so E_projected ≤ E_BS, and non-integer multiplicities are refused. The quadratic regex is now [-+]?(?:\d+(?:\.\d*)?|\.\d+)(?:[EeDd][-+]?\d+)? — 0.007 s at 160k chars, was ~370 s; my first attempt reproduced the defect and the new test caught it. The dying-ORCA-job case is fixed: check_stability_job gates on the log existing rather than job_status == 'done'. species.py's docstring corrected.

**11. Reference flip scope - agreed/fixed. You were right about the 1.99kcal/mol and that CCSD(T) is not variational. With Molpro 2026, it showed it was larger than 1.99 for the C5H10 singlet TS:

reference             SCF            CCSD(T)-F12     T1        D1
RHF                   -194.974539    -196.064366     0.0410    0.2172
broken-symmetry UHF   -195.077240    -195.840652     0.0124    0.0594
triplet ROHF          -195.066148    -196.053302     0.0118    0.0369

As we see, the broken symmetry reference is 64kcal/mol lower at SCF and 140kcal/mol higher after coupled cluster, recovering 0.35Eh less correlation. The fourth column shows that T1 falls from 0.0410 to 0.0124, below ARC's own 0.015 threshold.

Adoption is now gated positively on the level BROKEN_SYMMETRY_METHOD_TYPES = ['dft'] and BROKEN_SYMMETRY_METHODS = ['hf', 'rhf', 'uhf', 'rohf'], behind level_admits_a_broken_symmetry_reference(). Not by adding 'wavefunction' to REFERENCE_AGNOSTIC_METHOD_TYPES, which fails twice: hf resolves to method_type == 'wavefunction' (deduce_method_type puts it in wave_function_methods), so it would silence Hartree-Fock too; and a reference-agnostic level makes job_scf_reference_is_restricted return None, so record_scf_reference would never record a single point and the mixed-reference detection would be disabled entirely.

opt and freq still adopt, so the re-optimisation onto the broken-symmetry geometry and the ZPE taken from it are unchanged. The single point stays restricted and emits no BrokenSym. E0 therefore mixes references, which is the lesser inconsistency and one the branch already had machinery for - check_scf_reference_consistency now fires, output.yml reports reference_mismatch: true with MIXED_SCF_REFERENCE_MESSAGE, and T1 stays loud. check_scf_reference_consistency's docstring said "an adopted stability verdict does not reach this check", which this makes false; it is rewritten.

One consequence to be aware of: the documented default - Gaussian DFT opt/freq with a Molpro ccsd(t)-f12 single point - previously refused adoption outright and reported it. It now adopts for the geometry and the ZPE, keeps the energy restricted, and reports the mismatch instead. advanced.rst is updated, including deleting the measured BrokenSym-moves-the-energy-by-1.99-kcal/mol claim, which described a path that no longer exists.

Also, another thign was found and fixed. stable=(rext,noopt)` can report RHF -> CRHF, and adopted_reference_is_unrestricted had no relaxation-class gate - such a verdict would have run as a real unrestricted determinant chasing a complex solution. Adoption now requires derived_instability_breaks_spin_symmetry(species) is not False.

[alongd]

Round 3 — scoped to the broken-symmetry adoption gate

Reviewed only what changed since afc47f9 in response to findings 1 and 11:
level_admits_a_broken_symmetry_reference(), BROKEN_SYMMETRY_METHOD_TYPES,
BROKEN_SYMMETRY_METHODS, and Scheduler.stability_verdict_can_be_honoured.
Worktree at f95cd83; arc/job/adapters/common_test.py arc/scheduler_test.py arc/output_test.py → 308 passed in 8.11s.

Verdict

Keep the gate. Scoping adoption to levels whose energy is their determinant's energy is
the right call, and refusing to expand a correlated wavefunction about a symmetry-broken
reference is correct. The mixed-reference E0 it produces for the documented default
(Gaussian DFT opt/freq + Molpro ccsd(t)-f12 sp) is the lesser of the two inconsistencies,
and — checked, not assumed — it is reliably detected: molpro.py:233 calls
is_restricted(self), so the sp job carries a real restricted_used bool,
post_sp_actions records it with reference_key='sp' (scheduler.py:3993),
check_scf_reference_consistency fires, and output.py:354 guards reference_mismatch
against a None on either side. Restart and TS-switch survival are clean (details in 5).

Three things need to change before it merges, and the first two are both about the
justification, not the mechanism.

---

1. advanced.rst promises a reporting channel that does not exist

docs/source/advanced.rst:86-87, new in this PR:

ARC reports that in the log, in the species' output.yml warnings and in the run
summary rather than re-running the species.

It does not reach the run summary. arc/main.py:753-763:

for label, output in self.output.items():
    if output['convergence']:
        status_dict[label] = True
        logger.info(f'Species {label} converged successfully')
        if output.get('wavefunction_stability'):
            logger.info(f'  Wavefunction stability: {output["wavefunction_stability"]}')
        logger.info('\n')
    elif not label.startswith('IRC_'):
        ...
        if 'warnings' in output and output['warnings']:
            logger.info(f'  and warnings: {output["warnings"]}')

output['warnings'] — which is where MIXED_SCF_REFERENCE_MESSAGE lives — is printed
only in the failure branch. A species that mixes references converges, so it prints
Species X converged successfully and one stability line, and the mismatch is not in the
summary at all.

This is the whole configuration the gate was designed for. And it is a three-line fix in
a block this PR already edits — the diff against 735de21 shows main.py changed exactly
these lines to add the stability line. Please add the warnings alongside it, or delete the
claim from the docs.

Related, and this is the honest answer to "does anything downstream consume E0 without
seeing reference_mismatch": nothing consumes it. grep -rn 'scf_references| reference_mismatch|MIXED_SCF_REFERENCE' arc/ --include=*.py outside scheduler.py returns
only arc/output.py and arc/species/species.py. arc/processor.py never reads
warnings or scf_references; the thermo library entry and the rate coefficient are
written regardless. That is a defensible choice — ARC doesn't invalidate on this class of
warning anywhere — but combined with the summary gap it means the only in-band trace of a
mixed-reference barrier is one logger.warning at the moment the sp lands, thousands of
lines up a 200-species log. The summary line is what makes it real.

2. The 140 kcal/mol measurement should not ship as documentation

advanced.rst:76-80 states it as a fact. I can't run Molpro, but the three rows refute each
other by arithmetic:

reference	SCF	CCSD(T)-F12	E_corr
RHF	−194.974539	−196.064366	−1.089827
BS-UHF	−195.077240	−195.840652	−0.763412
ROHF triplet	−195.066148	−196.053302	−0.987154

The BS-UHF reference absorbs 0.1027 Eh (64.4 kcal/mol) of static correlation at the SCF
level, yet the correlated treatment loses 0.3264 Eh (204.8 kcal/mol). 0.2237 Eh goes
missing. CCSD(T) is not variational, but a 140 kcal/mol spread between two single-reference
CC treatments of the same singlet state at the same geometry and basis means at least one of
them is not a CCSD(T)-F12 energy of that state — not when T1(RHF) = 0.041 puts the
multireference character in "watch it" territory, not "the method has failed by 0.2 Eh"
territory.

The sharpest form: BS-UHF recovers 0.22 Eh less correlation than the ROHF triplet row,
which is also an open-shell F12 treatment of the same molecule in the same basis. Same code
family, same electron count. A spin-relaxed reference legitimately recovers a bit less; it
does not lose a fifth of a Hartree.

The likely mechanism is ARC-relevant, so worth naming: Molpro's ccsd(t)-f12 is the
closed-shell F12 code and uccsd(t)-f12 is the ROHF-orbital open-shell one. Neither takes a
genuinely spin-broken UHF determinant as its reference. {uhf; rotate,...} orbitals fed to
either gives a CC expansion about something that is not the BS determinant.

The gate does not need this number. The T1/D1 row does all the work and is independently
credible: 0.0410 → 0.0124 and D1 0.2172 → 0.0594 is exactly the documented behaviour of a
reference that has absorbed static correlation into its orbitals, and the diagnostic that
exists to catch a multireference species stops catching it. Note that the ROHF triplet row
demonstrates your point more cleanly than the BS row does — T1 = 0.0118 with a perfectly
ordinary total energy, i.e. spin relaxation suppresses T1 without any of the 140 kcal/mol
weirdness. Cite that, drop or heavily qualify the 140 kcal/mol and 0.35 Eh claims, and the
argument gets stronger, not weaker.

3. Double hybrids admit a broken-symmetry reference, and by your own argument they must not

advanced.rst:71-76 states the principle:

The energy of those levels IS the energy of their SCF determinant... A correlated
wavefunction level keeps its restricted reference, because its energy is a correlation
expansion built about a spin-adapted reference rather than the energy of that reference.

Executed against the shipped gate:

b2plyp         type=dft            admits=True
b2plypd3       type=dft            admits=True
dsd-pbep86     type=dft            admits=True
wb97x-2        type=dft            admits=True
hf             type=wavefunction   admits=True
hf-3c          type=wavefunction   admits=False
ccsd(t)-f12    type=wavefunction   admits=False

A double hybrid's energy is not its determinant's energy — it carries an MP2 term expanded
about the KS determinant. That is the same construction the docstring excludes, and your own
objection applies with full force: BS-MP2 about a spin-broken KS reference is the textbook
pathology. BROKEN_SYMMETRY_METHOD_TYPES = ['dft'] lets B2PLYP, DSD-PBEP86 and ωB97X-2 do
precisely what CCSD(T)-F12 is forbidden to do.

Not hypothetical for ARC: double hybrids are a normal sp_level choice, and for a
double-hybrid sp the mismatch check would also report reference_mismatch: false — one
reference throughout — while the number is built the way this gate exists to prevent.

Suggest excluding double hybrids explicitly (a DOUBLE_HYBRID_METHODS deny-list checked
before the method_type test is the smallest fix). hf-3c landing on admits=False is the
mirror image and much less consequential — worth a line either way so the two lists are
symmetric.

4. Confirmed as you described it

deduce_method_type really does put hf in 'wavefunction' — arc/level.py:344-371,
'hf' matches by substring — so Level(method='hf').method_type == 'wavefunction' and
admits=True only because of the BROKEN_SYMMETRY_METHODS fallback. Your reason for gating
positively rather than adding 'wavefunction' to REFERENCE_AGNOSTIC_METHOD_TYPES holds on
both counts: it would silence HF, and a reference-agnostic level makes
job_scf_reference_is_restricted return None, so record_scf_reference would never record
the sp and the mismatch detection would go dark entirely. Verified by reading both paths.

5. Restart and TS switch — checked, and clean

Both survive:

• scf_references round-trips through the restart file (species.py:829-830 writes it when
  non-empty, :972 restores it) and job.restricted_used is restored at
  scheduler.py:5067-5070. A species that adopts for opt/freq, gets restarted, and then runs
  its sp still reports the mismatch.
• carry_stability_verdict_across_ts_switch resets species.scf_references = dict()
  (:3879) and strips MIXED_SCF_REFERENCE_MESSAGE (:3881-3884), so a freq reference
  from an abandoned guess can never be compared against a new guess's sp. And a verdict
  carrying reference_change_available: False can't leak into the rebuilt dict at
  :3901-3905, because adopted_reference_is_unrestricted returns False for it and the
  drop branch fires first. That's tighter than it needed to be; good.

6. Smaller

• UNREACHABLE_REFERENCE_MESSAGE is not stripped on a TS switch. The list at :3881 is
  [MIXED_SCF_REFERENCE_MESSAGE, INVALID_ANALYTIC_FREQ_MESSAGE, SPIN_CONTAMINATION_MESSAGE, COLLAPSED_REFERENCE_MESSAGE]. UNREACHABLE_REFERENCE_MESSAGE, appended at :2252, isn't
  in it, so it outlives the guess it described while output[label]['wavefunction_stability']
  is set to None on the same pass. The method's own docstring says "THE TWO RECORDS ARE
  REDUCED TOGETHER".
• post_sp_actions at scheduler.py:1626-1629 still passes no job. Round-2 item 12.
  The blast radius is bounded — that branch only runs when sp_level == opt_level, where a
  reference mismatch is impossible by construction — so the effect is reference_mismatch: null in output.yml for a restarted single-level project rather than a wrong number.
  Worth a comment saying so, since the guard reads like an oversight.
• 'rhf' in BROKEN_SYMMETRY_METHODS is unreachable: Level(method='rhf', basis='cc-pvdz')
  raises IndexError: list index out of range at level.py:450 because no adapter registers
  rhf. Pre-existing (level.py is untouched here), not yours to fix, but the entry is dead.
  'uhf'/'rohf' are constructible but moot, since adoption only fires for a species whose
  reference was restricted.

---

Ranked for merge: (1) fix the run-summary claim — code or docs, either is fine, but they
have to agree; (2) strike or qualify the 140 kcal/mol and 0.35 Eh claims and lean on T1/D1;
(3) exclude double hybrids. (6) is cleanup. The gate itself, and the decision to accept a
mixed-reference E0 over a correlated energy expanded about a symmetry-broken determinant, I'd
merge as designed.

[calvinp0]

1. Run summary - Fixed in code. output['warnings'] now prints in the ocnverged branched alongside the stability line, so a mixed reference species that converges normally shows the mismatch. Test added. Your point that nothing downstream consumes reference_mismatch is confirmed and I've left that as-is: ARC doesn't invalidate on this class of warning anywhere, and inventing a new invalidation path here would be inconsistent.
2. The 140 kcal/mol measurement -withdrawn, and you were right about the mechanism. I checked the arithmetic before deleting anything. The rows are internally consistent (E_corr = E_CC − E_SCF to ~1e-14, so not a transcription error) and still wrong: BS-UHF recovers 0.2237 Eh less correlation than the ROHF triplet and 0.3264 less than RHF, while RHF and ROHF-triplet totals agree to 0.011 Eh - two credible rows and a 20× outlier. And the direction is wrong: BS-UHF is variationally 64 kcal/mol below RHF at SCF, so CC on it cannot land 140 kcal/mol above. The prose also miscopied its own table - 0.35 Eh against 0.3264.
   Your mechanism holds: ccsd(t)-f12 is the closed-shell program and uccsd(t)-f12 the ROHF-orbital one - the "U" names spin-unrestricted amplitudes, not orbitals - so neither takes a spin-broken determinant and the difference subtracted two different references' energies.
   One correction to your constructive half: the T1/D1 pair you suggested leaning on (0.0410 → 0.0124, D1 0.2172 → 0.0594) is the BS-UHF row. A diagnostic read off a calculation that isn't an energy of the state is worth no more than the energy, so I withdrew those with it. The docs now argue from the ROHF triplet row alone - exactly your point, and it's the stronger version: T1 0.0410 → 0.0118, below threshold, with the total energy moving under 7 kcal/mol. The docs also say plainly that no BS-UHF correlated energy is quoted because none exists, and why.
3. Double hybrids - excluded. DOUBLE_HYBRID_METHODS deny-list read before the method_type test, as you suggested. b2plyp, b2plypd3, dsd-pbep86, wb97x-2, xyg3, pbe0dh and the rest now return False; ordinary DFT is unchanged. Your reasoning is the same one the docstring already used for correlated wavefunction levels, and it applied here with full force — a double hybrid's energy is not its determinant's energy. ARC had no double-hybrid list anywhere; the only acknowledgement was a comment in data/ess_methods.yml whose grouping has since drifted (PW6B95, M08HX, MN15 sit inside the marked block and aren't double hybrids). hf-3c added to BROKEN_SYMMETRY_METHODS so the two lists are symmetric.

6a. UNREACHABLE_REFERENCE_MESSAGE added to the strip list. It's always safe to strip: that warning is raised only on a verdict stamped reference_change_available: False, which the drop branch already reads, so such a verdict never survives the switch anyway.
6b. Docstring note added recording that the job-less caller is the restart path, reachable only where sp_level == opt_level, and that the cost is a null rather than false reference_mismatch.
6c. 'rhf' isn't dead — kept. The IndexError you hit comes from deduce_software running when no software is given. rhf is registered for TeraChem (data/ess_methods.yml:136), and Level(method='rhf', basis='cc-pvdz', software='terachem') builds fine, types as wavefunction and reaches the gate with admits=True. Removing the entry would break a reachable case.

Separately, your item 3 turned up a main bug worth its own fix. Checking why cam-b3lyp reported admits=False, it turns out deduce_method_type matches semiempirical_methods = ['am', 'pm', …] as substrings, so 'am' inside 'cam' fires. Every Coulomb-attenuated functional — cam-b3lyp, camh-b3lyp, cam-qtp00, lc-camb3lyp — types as semiempirical. That matters beyond this gate: is_species_restricted early-returns True for semiempirical before any multiplicity check, so a doublet radical at CAM-B3LYP currently runs restricted (b3lyp → restricted=False, cam-b3lyp → restricted=True, same species). level.py is untouched here so I haven't fixed it in this PR.

[Copilot]

Pull request overview

Copilot reviewed 42 out of 43 changed files in this pull request and generated 1 comment.

Suppressed comments (1)

docs/source/input_reference.rst:268

• This paragraph still frames the stability job as consuming the checkfile a frequency job used and a Hessian that was already built. Since stability runs before the frequency calculation, this should be phrased in terms of what the frequency job will use/build so the timeline is accurate.