Displace the TS along the normal mode in Cartesian, not mass-weighted, coordinates

[calvinp0]

⚠️ This PR changes a physical quantity. It has had a quantum-chemistry review in addition to a code review. The calibration argument in Part 3 is the part to read closely.

Base: main. This branch absorbs #967 (frequency-log geometry frame) and #968 (ESSs that report no normal modes), both flattened into this PR and now closed. Three commits, no file touched by more than one of them.

All three fix defects at the head of analyze_ts_normal_mode_displacement(), and each of them misplaces the same displaced geometry x + A·d that the entire analysis is measured on. They are split by cause, not by file.

#	cause	symptom
1	modes reported as (None, None) reach a subscript	uncaught TypeError ends the whole ARC run
2	geometry and modes come from different coordinate frames	verdict is not rotation-invariant
3	the displacement is mass-weighted where it must be Cartesian	measures a different quantity, then reorders the bonds

---

Part 1 — an ESS that reports no modes must not kill the run

• parse_normal_mode_displacement() returns the sentinel (None, None) when it cannot read displacements, and every consumer in ARC treated that tuple as data. Only 3 of ARC's 9 concrete adapters return real data (Gaussian, xtb, YAML); CFour, Molpro, Orca, Psi4, QChem and TeraChem return the sentinel for every file, and Gaussian, xtb and YAML return it for a file with no modes.
• analyze_ts_normal_mode_displacement() guarded the parse with except NotImplementedError, which no adapter ever raises, so the guard was dead code and (None, None) fell through to normal_mode_disp[0] → TypeError: 'NoneType' object is not subscriptable. Nothing up the chain catches it: Scheduler.check_freq_job → post_freq_actions → check_ts → check_normal_mode_displacement → the analysis holds no try/except, nor does Scheduler.schedule_jobs, Scheduler.__init__, ARC.execute or ARC.py::main(). One Orca frequency job was enough to end the run.
• The same defect sat at a second site. arc/job/trsh.py::trsh_negative_freq() unpacked the sentinel and evaluated len(normal_modes_disp) → TypeError: object of type 'NoneType' has no len(). Reproduced on the shipped fixtures freq/orca_neg_freq_ts.out, freq/orca_example_freq.log, freq/orca6_example.out, freq/CH2O_freq_molpro.out, freq/C2H6_freq_QChem.out, freq/CH2O_freq_terachem.dat. It is reachable for any non-TS species: Scheduler.check_freq_job() and parse_composite_geo() both call troubleshoot_negative_freq() when a frequency job converges with an imaginary frequency and trsh_ess_jobs is on. Frequencies parse for those files while displacements do not (orca_neg_freq_ts.out reports 15 frequencies, minimum −1271.62 cm⁻¹), so the troubleshooter is entered on exactly the files whose displacements are missing. It now reaches its existing "Could not troubleshoot negative frequency" branch and the species is left to the ordinary ESS troubleshooter.
• New parser.get_normal_mode_displacement() returns the frequencies and displacements, or None, warning and naming the ESS via the existing determine_ess(). It is keyed on the returned value, not on a list of which adapters implement the method, so an adapter that gains parsing is picked up with no change here, and a supported ESS whose log happens to hold no modes is skipped by the same path. It takes a log file path rather than a JobAdapter because trsh_negative_freq() receives a path, and returns the frequencies alongside the displacements because that caller needs both.
• Why not fix it in make_parser(). Making (None, None) trigger the documented NotImplementedError is the more correct contract but does not fix the bug: raise_error defaults to False and no production caller passes True, so callers would still receive the sentinel. Raising unconditionally would change the return contract of all fifteen parse_* entry points at once, including ts.py::get_rxn_zone_atom_indices(), which passes raise_error=False precisely to keep going.
• Why the helper lives in arc/parser/parser.py — beside parse_normal_mode_displacement() whose sentinel it interprets and beside determine_ess() which it calls. Both call sites already import that module on main (arc/checks/nmd.py as from arc.parser import parser, arc/job/trsh.py as from arc.parser.parser import ...), so one shared implementation is reached at module level with no new import edge. Defining it in arc/checks/nmd.py would have made arc/job/trsh.py import arc.checks.nmd, an edge main does not have and one that closes a cycle, since nmd.py's imports reach arc.job.trsh back through arc/__init__.py's eager package imports (CodeQL py/unsafe-cyclic-import). Measured, not asserted: every module's top-level imports were walked on origin/main and on this branch (following if TYPE_CHECKING and try/except, excluding function bodies) — the production import graph is identical; the only difference anywhere is the test module arc/checks/nmd_test.py. Import order was also checked both ways in fresh interpreters run from an isolated directory.

---

Part 2 — the geometry must come from the frequency log

• The geometry came from reaction.ts_species.get_xyz(), i.e. final_xyz, which only the opt handlers (Scheduler.parse_opt_geo, parse_composite_geo) write and which check_freq_job never refreshes. The modes came from the frequency job's log, and Gaussian reports them in that job's own standard orientation. When get_xyz() falls through to initial_xyz, a conformer or a TS guess, the frames are unrelated altogether.

• Measured over 996 real benchmark frequency jobs: 82% agree to within 1°, but 17.6% differ by ~180° — exact rigid rotations (fit RMSD ~1e-16), not numerical noise. Adding a mode vector expressed in one frame to a coordinate expressed in another corrupts the cross term of every displaced bond length, so the rotation-invariance of a bond length does not rescue the result. Any non-identity rotation breaks it, proper or improper; a pure translation is harmless.

• get_ts_xyz_in_normal_mode_frame() sources the geometry from the same file the modes are parsed from. It falls back to the species geometry, with a warning, when that file yields no geometry, when parsing it raises, or when the parsed element symbol sequence differs from the species' one, so mismatched atoms are never silently compared. The fallback is one-directional: when the species has no geometry to compare against, the parsed geometry is used rather than discarded — it is the geometry the analysis exists to obtain, and the caller immediately does len(ts_xyz['symbols']).

• Per-ESS behaviour, verified against the shipped fixtures:

  ESS / path	outcome
  Gaussian	parse_geometry() returns the last Standard orientation: block — the frame the Frequencies -- vectors are printed in. On freq/TS_CH4_OH.log it reproduces the species geometry to a 0.000° rotation (fit RMSD 1e-16), so this is a no-op on the existing tests
  nosymm	with no Standard orientation:, falls back to Input orientation:, which is the mode frame when Gaussian does not reorient and is also the geometry ARC submitted, i.e. final_xyz
  composite / optfreq	parse_composite_geo() and parse_opt_geo() set final_xyz from the same log check_freq_job() then reads, so those frames already agreed
  xtb	modes come from a sibling g98.out and parse_geometry() returns None for the freq output.out of normal_mode/HO2 and normal_mode/TS_0, so the fallback preserves previous behaviour exactly
  Orca, Q-Chem, Molpro, TeraChem, CFOUR, Psi4	return (None, None), so the analysis cannot run either way. Molpro's parse_geometry() raises TypeError, which the fallback absorbs rather than turning into a new failure mode

• Independently confirmed by the Eckart condition. Modes expressed in their own frame satisfy Σ mₐ(rₐ − R_com) × dₐ = 0, and that quantity is frame-dependent, so it directly tests whether the pairing is consistent: 73% of affected jobs exceed 3σ in the opt frame, versus 0 of 111 in the freq frame.

• Fixed in nmd.py, not by refreshing final_xyz in the scheduler, which keeps the blast radius to the one consumer that needs the mode frame. final_xyz is the geometry ARC reports, saves to the restart and output files, passes to Arkane and feeds to every subsequent job; rewriting it with a reoriented copy would rotate all of those for no benefit.

• Correcting the frame exposed a precondition that was never enforced. The forming, breaking and changed bond indices index into the concatenated reactant geometry, while a TS geometry may order its atoms differently — heavy atoms first, say, where the reactant concatenation interleaves them. The indices then address the wrong atoms entirely, and the resulting bond lengths can happen to agree with the expected pattern, so the check returned a confident verdict about atoms it was not looking at. Sharpening the geometry turns some of those accidental passes into equally unfounded rejections. is_ts_atom_order_consistent_with_reactants() compares the TS element symbol sequence against get_reactants_xyz(), and the analysis now returns None when they differ. Canonicalising the atom order belongs at ingestion in the TS adapters and is not attempted here.

None, not False, and why it is load-bearing

Scheduler.post_freq_actions() calls switch_ts() when ts_checks['NMD'] is False. Returning False for either new skip path would discard converged transition states on the strength of indices known to be meaningless, or search for a replacement that ARC cannot read either way. None is what distinguishes "could not check" from "checked and failed", and is also the value populate_ts_checks() initialises the entry to, so an unrun check and an uncheckable one now agree. An ESS that cannot report modes is telling us nothing about the TS; it must not count as evidence against it.

---

Part 3 — the displacement was mass-weighted where it must be Cartesian

get_displaced_xyzs() scaled each atom's displacement by √mₐ before adding it to a Cartesian geometry. That is a unit error — it adds vectors from two different spaces. Removing it alone would shrink the probe step 3.46× for carbon and start rejecting genuine transition states, so the amplitude moves with it: 0.25 → DEFAULT_AMPLITUDE = 0.9, since 0.25 × √12 = 0.87 ≈ 0.9. The recalibration restores the original author's effective scale rather than overriding it — keeping 0.25 is what would break it. The two changes are one change; shipping either alone is wrong.

Why √m has to go — five independent lines

1. Gaussian prints Cartesian displacements. Per Vibrational Analysis in Gaussian (Ochterski, https://gaussian.com/vib/), and provable from the file itself: mass-weighted eigenvectors L satisfy Σ|L|² = 1; the Cartesian form is lₐ = Lₐ/√mₐ; the reduced mass is defined μ ≡ 1/Σₐ|lₐ|²; Gaussian prints dₐ = √μ·lₐ, so Σ|d|² = 1. Therefore Σₐ mₐ|dₐ|² = μ exactly — an identity that holds only for Cartesian unit-normalised vectors. Checked against the reduced mass the same files print, over 136 real modes: median error 0.77%, 99.2% within 2%.
2. So the old line added vectors from two different spaces.
3. It moved the centre of mass. A vibration is orthogonal to translation (Σₐ mₐdₐ = 0). COM drift per unit amplitude: 0.001 Å median as printed vs 0.037 Å median / 0.299 Å max with √m — a 35× discrepancy.
4. It reordered the bonds. Correlation with the first-order rate û·(dᵢ−dⱼ) over 927 bonds: 0.99999 as printed vs 0.789 with √m. For a test whose entire job is to rank reactive against spectator bonds, reordering is fatal.
5. The algebra shows it was never a rescaling — √mᵢ·dᵢ − √mⱼ·dⱼ = √mᵢ(dᵢ − dⱼ) + (√mᵢ − √mⱼ)·dⱼ. The second term is proportional to the absolute displacement of one atom, not the relative one, so a rigidly translating fragment (true rate exactly zero) reports non-zero for any heteronuclear bond. Measured |spurious|/|true| over 456 bonds: median 2.46, greater than 1 for 63.4% of bonds; the best single scale factor per structure still leaves a median residual of 20.6%.

√m did not scale the measurement — it measured a different quantity. No choice of amplitude could compensate. That is why this is a correctness fix, not a tuning disagreement.

The weights argument is removed from get_displaced_xyzs() rather than defaulted off, because no per-atom weight is meaningful when a mode is added to a Cartesian coordinate. The weights analyze_ts_normal_mode_displacement() computes still scale the compared bond lengths, which is a separate use and is left unchanged.

Why the amplitude moves with it

• A bond counts as reactive only if its length changes by more than 5% of the bond length, and that change is proportional to the amplitude — so the amplitude sets the effective threshold on the underlying rate.

• Corroborated on IRC-validated structures. Acceptance floors, bisected, verdict monotone in amplitude:

  structure	ν (cm⁻¹)	acceptance floor
  reaction_08/freq_a5382	−1743.2	0.676
  reaction_08/freq_a9411	−1743.2	0.676
  r3_07/freq_a1758	−1466.8	0.680

  All three are rejected at 0.5 — and these are the transition states reaction_08 and r3_07 published kinetics for. 0.9 clears the 0.68 floor with margin.

• Upper edge. On the labelled C3H7 fixture, separation holds from 0.1 through 1.5 and breaks at 2.0. Independently, agreement between the finite-step bond change and the first-order rate stays above r = 0.99 for all 109 structures up to amplitude 1.0, degrading for 1.8% at 1.25, 11% at 1.5 and 24% at 2.0 as second-order contamination starts reordering bonds. The usable band is roughly 0.68 – 1.25; 0.9 sits inside it.

• The amplitude is the module constant DEFAULT_AMPLITUDE, used by both the analysis default and check_normal_mode_displacement() in arc/checks/ts.py, so the two cannot drift apart again.

What 0.9 does not do

• It does not separate genuine from doubtful. The last genuine acceptance is at 0.680 and the first doubtful one at 0.682 — 0.002 apart. That is a coincidence, not a threshold, and the constant should not be read as one.
• Two structures are admitted at 0.75 and 0.9 alike and need IRC to adjudicate: r3_05/freq_a779 (ν = −662.8) and r3_15/freq_a666 (ν = −333.6). r3_15 looks wrong on mode character — labelled an H-migration, yet μ = 3.22 and the migrating H is not among the top three movers.
• The required amplitude tracks how much of the imaginary mode sits in the secondary reactive coordinate, which is family-dependent: a clean single-H-transfer control has a floor of 0.467, while the ring-forming cases need 0.68. A single global constant may be the wrong shape long-term — a follow-up, not this PR.

---

How it was verified

• Both new physical tests fail on the unfixed code, with the predicted numbers:
  • test_get_displaced_xyzs_conserves_the_center_of_mass → 0.03234 not less than 0.005
  • test_get_displaced_xyzs_moves_heavy_atoms_less_than_hydrogens → 0.008660254 != 0.0025, i.e. 0.008660 / 0.0025 = 3.4641 = √12, the carbon inflation measured directly.
• A third test pins the separation itself — the property the amplitude governs, which the other two do not cover because they assert displacements rather than verdicts. On the labelled C3H7 fixture (TS3 is the correct TS for iC3H7 <=> nC3H7; TS4 and TS7 are wrong saddle points on the same PES) it is two-sided: DEFAULT_AMPLITUDE = 0.02 → AssertionError: False != True : TS3 (correct TS rejected); = 2.0 → AssertionError: True != False : TS4 (wrong saddle accepted).
• Part 2: frames agree, a deliberately rotated geometry, an element symbol mismatch, a log with no geometry, a log that fails to parse, a species with no geometry at all, the atom-order guard (unit and end to end), and a scheduler-level test that an unknown (None) verdict does not trigger a TS switch.
• Part 1: the helper's own contract is tested in arc/parser/parser_test.py over a Gaussian, a YAML and an xtb log for the parsed case and over the six unreadable fixtures plus freq/yml_no_freqs.yml for the None case; nmd_test.py covers the analysis returning None rather than raising for the six unreadable fixtures and still reaching a bool verdict for an xtb log; trsh_test.py asserts trsh_negative_freq() returns its four empty lists over those six fixtures rather than raising, proved failing before the fix with the TypeError quoted above.

What changed in the existing tests

• 5 frozen golden-coordinate comparisons were removed. Every one had been computed from the buggy formula, so "correcting" them would have re-frozen numbers derived from the thing under repair. They are replaced by the definitional assertion x ± A·d, which pins the full coordinate array (atol=1e-10) and cannot drift with the implementation, plus the physical tests above.
• The σ golden in the integration test became a bound, assertGreater(float(sigma), 10 * nmd.SIGMA_THRESHOLD), i.e. σ > 30. The removed byte-exact golden was σ = 14.7891 — the value the √m code produces, so it had to be rewritten rather than checked by this PR; the corrected code gives σ ≈ 122 on this fixture (105.7 – 128.0 across the platforms measured — BLAS differences move the trailing digits, which is why it is a bound and not a golden). The bound sits far above the regression value and far below the corrected one, so it fails loudly if the √m weighting ever comes back, without enshrining a float.
• Downstream scalar goldens (baseline, std, the reactive bond diffs) were recomputed rather than removed, and checked for a preserved verdict.

---

Deliberately not addressed here

• A second, independent √m displacer. arc/species/converter.py::displace_xyz() does the same operation with the same weighting (use_weights=True by default); its production caller is trsh_negative_freq(), which perturbs a geometry to seed a fresh optimization. The Part 3 argument does not transfer unchanged: that code does not measure relative bond rates, it only needs a step that leaves the saddle, and its amplitude ladder (0.25 → 2.5) is itself calibrated against the weighted step. Correcting it is a separate change with its own calibration.
• The surviving bond-length weighting. get_weights_from_xyz() remains live; its result now feeds only get_bond_length_changes() and get_bond_length_changes_baseline_and_std(). Whether that weighting is correct is an open question this PR neither settles nor endorses. (mᵢmⱼ)^¼ is not √μ and is not a reduced-mass proxy — it is √(wᵢwⱼ) with w = √m, the geometric mean of √mass; normalised to C–C it gives C–H 0.538 where √μ gives 0.394. Measured over 99 structures it inflates the reactive bond in 47.5% of cases and suppresses in only 5% (median ratio 1.000), so it biases toward false accepts, not against X–H.
• Two absolute-scale constants interact with the amplitude and are not recalibrated. STD_FLOOR = 1e-4 clamps the spread: 53/99 structures have a MAD below it (19/99 have a single background bond, MAD ≡ 0), so for those σ is (min − baseline)/floor and therefore scales with the amplitude rather than being invariant to it. DIRECTIONALITY_MIN_DELTA = 0.005 Å is likewise absolute. A principled pass should address all three together.
• arc/checks/ts.py::get_rxn_zone_atom_indices() passes the same sentinel into get_rms_from_normal_mode_disp(). Left alone because it needs a decision about what an empty reaction zone should mean for the caller, not a guard — unlike the two sites fixed here, which both already had a "give up cleanly" branch that the crash was jumping over. It can reuse get_normal_mode_displacement() when it is fixed.
• An adjacent pre-existing bug: the analysis takes normal_mode_disp[0] unconditionally and never checks the sign of the corresponding frequency, so on a structure with no imaginary frequency it silently analyses the lowest real mode and returns a confident verdict. The guard would be to return None unless freqs[0] < 0.
• The atom-count gate just above the new guards, if n_ts != n_expected: return False, is the same class of missing precondition but still returns False on this branch. It is addressed separately in Support reactions where a species participates more than once (A + A) #974.
• A scoring redesign. Scoring each bond by the first-order rate dr_ij/ds = û_ij·(dᵢ−dⱼ) is amplitude-free and mass-weighting-free, which would remove the calibration coupling in Part 3. It does not subsume Part 2: both criteria are already exactly invariant to a common rigid rotation, whereas Part 2 concerns relative misorientation between a geometry taken from one file and a mode taken from another.

What I searched for

• Every routine in ARC that displaces a geometry along a normal mode, by behaviour rather than name. There are two: get_displaced_xyzs() (fixed here) and converter.py::displace_xyz() (scoped out above). They are not consolidated because their callers want different things — one measures, one perturbs.
• Whether weights=True was a deliberate, documented choice: the full history of arc/checks/nmd.py, every docstring, and the tests. The commit that introduced it has an empty body and nothing justifies it.
• Whether a helper already resolved a geometry into the frame of a job's parsed modes — trsh_negative_freq does it inline for a different purpose; there was no reusable helper to import.
• Whether any "is this ESS supported for this parse method" predicate or registry of adapter capabilities exists. None does — settings['supported_ess'] lists the programs ARC can run, not what each parser adapter can read. determine_ess() is reused for the ESS name in the warning.
• Reused get_element_mass from arc/common.py rather than adding a mass lookup, and routed vector-length maths in the new tests through arc/species/vectors.py::get_vector_length rather than inlining np.linalg.norm.

🤖 Generated with Claude Code

---

Also folds in feature_nmd_family_recipe

The reaction atom map is only determined up to permutations of chemically equivalent atoms, and a valid-but-permuted map inflates the map-derived formed/broken sets with spurious pairs no real mode can satisfy. On CCO[O] <=> C=C + [O]O the map crosses the four ethylene spectator hydrogens, yielding 5 formed + 6 broken instead of 1 formed + 2 broken + 1 changed — and NMD rejected the genuine −1074 cm⁻¹ concerted-elimination saddle.

ARCReaction.get_reactive_bonds_from_family() now derives the reactive bonds straight from the RMG family recipe and r_label_map, in reactant index space, using only forward-discovered product dicts — a reverse-discovered r_label_map is in product space and produced three non-existent bonds, two of them H–H, on a Retroene case. get_bond_change_candidates() tries that first and keeps the map-derived set as a fallback, which also lets NMD run at all when atom mapping fails. Two Gaussian frequency fixtures ship with it.

One interaction with this PR's recalibration, worth noting. The folded branch's HO₂-elimination test pinned amplitude=0.25. With the √mass factor removed and the default moved to DEFAULT_AMPLITUDE = 0.9, the genuine saddle is rejected at 0.25 — measured good/wrong = False/False at 0.25, then True/False at 0.5, 0.676, 0.9, 1.0 and 1.25. The test now probes at nmd.DEFAULT_AMPLITUDE and still discriminates: with get_bond_change_candidates reduced to the map-derived candidate, the genuine saddle is rejected at 0.9.

[codecov]

Codecov Report

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

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #970      +/-   ##
==========================================
+ Coverage   64.55%   64.66%   +0.11%     
==========================================
  Files         119      119              
  Lines       39788    39872      +84     
  Branches    10307    10326      +19     
==========================================
+ Hits        25684    25785     +101     
+ Misses      11123    11105      -18     
- Partials     2981     2982       +1     

Flag	Coverage Δ
functionaltests	64.66% <ø> (+0.11%)	⬆️
unittests	64.66% <ø> (+0.11%)	⬆️

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.

[Copilot]

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

[Copilot]

Pull request overview

Copilot reviewed 9 out of 11 changed files in this pull request and generated 2 comments.