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.

Top of a 3-PR stack. Reading order: #967 → #968 → this. It sits on the other two only because all three touch the head of analyze_ts_normal_mode_displacement(); it is independent in substance and is the only one of the three carrying physics, so it is the one to revert alone if the calibration is ever disputed.

The 30-second version

• get_displaced_xyzs() multiplied each atom's normal mode displacement by √mₐ before adding it to a Cartesian geometry. That is a unit error — it adds vectors from two different spaces.
• Removing it, on its own, 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. 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.
• This is not a judgement that the original calibration was wrong. The two changes are one change; shipping either alone is wrong.

Part 1 — why √m has to go

Five independent lines. None of them is a preference.

1. Gaussian prints Cartesian displacements, not mass-weighted ones. 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. √m·d is a mass-weighted displacement; the coordinates it was added to are Cartesian.
3. It moved the centre of mass. A vibration is orthogonal to translation (Σₐ mₐdₐ = 0), so displacing along one must not translate the molecule. 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 and not a tuning disagreement.

Part 2 — 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.

• √m was inflating the step for heavy atoms, so 0.25 behaved like ~0.87 for a carbon-dominated reactive coordinate (a ring-closure C···C, where both reactive atoms are carbon, is the maximum-effect case). Remove the factor and keep 0.25, and that probe shrinks 3.46×.

• 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 now 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.

Part 3 — 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 and fails if the constant leaves the safe band in either direction:
  • DEFAULT_AMPLITUDE = 0.02 → AssertionError: False != True : TS3 (correct TS rejected)
  • DEFAULT_AMPLITUDE = 2.0 → AssertionError: True != False : TS4 (wrong saddle accepted)
• arc/checks/ passes in full — 68 tests.

What changed in the 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 arc/job/trsh.py::trsh_negative_freq, which perturbs a geometry to seed a fresh optimization. This PR does not touch it. The Part 1 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 one production caller is analyze_ts_normal_mode_displacement(), and its result now feeds only the bond-length machinery (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. Removing it is a separate follow-up.
• Three absolute-scale constants interact with the amplitude, and this PR recalibrates only one. 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.
• An adjacent pre-existing bug: analyze_ts_normal_mode_displacement() takes normal_mode_disp[0] unconditionally and get_normal_mode_displacement() 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.
• 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 2. It does not subsume Use the frequency log's geometry for the TS normal mode displacement … #967: both criteria are already exactly invariant to a common rigid rotation, whereas Use the frequency log's geometry for the TS normal mode displacement … #967 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, searched by behaviour rather than name. There are two: get_displaced_xyzs() (fixed here) and arc/species/converter.py::displace_xyz() (scoped out above, with reasons). They are not consolidated in this PR 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.
• 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.

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 64.23%. Comparing base (7d37f89) to head (a74b0af).

Additional details and impacted files

@@                      Coverage Diff                       @@
##           fix_nmd_unsupported_ess_guard     #970   +/-   ##
==============================================================
  Coverage                          64.22%   64.23%           
==============================================================
  Files                                119      119           
  Lines                              39564    39565    +1     
  Branches                           10264    10264           
==============================================================
+ Hits                               25410    25414    +4     
- Misses                             11179    11181    +2     
+ Partials                            2975     2970    -5     

Flag	Coverage Δ
functionaltests	64.23% <ø> (+<0.01%)	⬆️
unittests	64.23% <ø> (+<0.01%)	⬆️

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.