MECP/MECI paired-state penalty scanner

[maxscheurer]

Closes #229

Add a paired-state MECP/MECI penalty-function scanner on top of the single-target NuclearGradientScanner from #221. It evaluates two electronic surfaces at one geometry from a single SCF + single ADC and exposes both energies and gradients to a Levine–Coe–Martínez penalty objective for geomeTRIC — no derivative couplings required.

Implementation plan posted as a comment below.

[maxscheurer]

Implementation Plan — Paired-state MECP/MECI penalty scanner

Analysis

Issue 229 adds a paired-state scanner on top of the single-target NuclearGradientScanner (PR 230). It evaluates two electronic surfaces at one geometry from one SCF + one ADC and returns both (energy, gradient) pairs, which a Levine–Coe–Martínez penalty objective combines into one (energy, gradient) for geomeTRIC — no derivative couplings.

Key exploration findings that shape the plan:

• nuclear_gradient is a pure function and is reusable as-is by calling it twice on two Excitations from the same run_adc result.
• The hard novel part is paired root tracking with a distinctness guard: the present _select_excitation computes one score vector and returns one root; nothing prevents two independent trackers collapsing onto the same root near the seam.
• geomeTRIC ships a built-in ConicalIntersection(Engine) implementing the exact LCM penalty — useful only as a validation oracle, not as the integration path: adcc's calc_new is not a geometric.engine.Engine subclass and production goes through PySCF's as_pyscf_method + geometric_solver.optimize callback bridge. Routing the penalty objective through that same existing bridge keeps adcc decoupled and testable without geomeTRIC (matching the "unit tests with controlled energies/gradients" acceptance criterion).
• Tests gate on adcc.backends.available() (PySCF) and importlib.util.find_spec("geometric"); end-to-end runs are separate (functionality_*_test.py), unit runs are PySCF-only (geomopt_scanner_test.py).

Resolved open question

Integration path → penalty-objective-only via the existing as_pyscf_method bridge. adcc owns a pure, geomeTRIC-free LCM penalty function and an MECPObjective wrapper returning a combined (energy, gradient), driving geomeTRIC exactly like today's single-surface scanner. This makes the penalty math unit-testable with controlled data and avoids coupling to geomeTRIC engine internals (Engine, dirname, ConicalIntersection). geomeTRIC's ConicalIntersection is used only as a numerical oracle in tests. The default penalty is the Levine–Coe–Martínez smoothed penalty (sigma, alpha parameters) with a documented raw energy-difference fallback mode.

Deliverables

1. Paired scanner evaluating two surfaces per geometry with paired root tracking + distinctness guard.
2. Pure penalty objective + an MECPObjective wrapper for the geomeTRIC bridge.
3. End-to-end MECI smoke test (slow/skipped) + fast pure-unit coverage of penalty math and paired-engine mechanics.
4. Docs subsection + runnable example.

Files to create / modify

• adcc/gradients/scanner.py — light, behaviour-preserving refactor: lift the geometry/SCF plumbing (_coords_array, _mol_at, _run_scf-style SCF lifecycle) and guards (_import_pyscf, _is_pyscf_scf, _PyscfModules, _check_ground_state_level, _mp_level, _safe_ao_ndarray) into reusable private helpers; the existing NuclearGradientScanner keeps identical behaviour (guarded by its current tests). No new user-facing API here.
• adcc/gradients/paired_scanner.py (NEW) — PairedStateGradientScanner: one SCF + one ADC + two nuclear_gradient calls; plural tracking state (previous_descriptors, previous_indices, last_trackings, last_excitations, last_gradients); joint distinct root selection (_select_excitations reusing the already-computed score vectors, picking two distinct indices with an energy-gap tie-break near degeneracy and a tracking_min_score/tracking_min_gap per channel); a distinctness guard that raises/diagnoses when the two tracked roots cannot be kept apart; __call__ returns ((e_lower, g_lower), (e_upper, g_upper)); calc_new returning both pairs; supports MECI (two excited indices from one ADC solve) and MECP (one GroundStateTarget MP2 + one excited root — gradient handled via nuclear_gradient(LazyMp) and LazyMp.energy(2), MP2 enforced by _check_ground_state_level).
• adcc/gradients/mecp.py (NEW) — mecp_penalty(e_lower, g_lower, e_upper, g_upper, *, sigma, alpha) pure function returning (E_pen, g_pen) (LCM smoothed penalty, matching geomeTRIC's formula for oracle cross-checks); raw edelta mode fallback. MECPObjective: callable wrapping a PairedStateGradientScanner + penalty into a single (energy, gradient) and a calc_new honoring the geomeTRIC flattened dict contract. No geomeTRIC import.
• adcc/gradients/__init__.py — export PairedStateGradientScanner, mecp_penalty, MECPObjective, and paired target dataclasses.
• adcc/__init__.py — top-level PairedStateGradientScanner (and a paired_nuclear_gradient_scanner convenience alias) mirroring NuclearGradientScanner placement.
• adcc/tests/geomopt_mecp_test.py (NEW, PySCF-only) — fast unit coverage.
• adcc/tests/functionality_geomopt_mecp_test.py (NEW, PySCF + geomeTRIC, slow/skipped) — end-to-end smoke.
• docs/gradients.rst — expand the one-line MECI note into a "Minimum-energy crossing points (MECP/MECI)" subsection with a worked example.
• examples/optimization/pyscf_adcc_mecp.py (NEW) — runnable twisted-ethylene MECI optimization via the as_pyscf_method bridge over MECPObjective.

Testing approach (mandatory)

Fast unit (geomopt_mecp_test.py, PySCF-only, mimics geomopt_scanner_test.py gates):

• Penalty math with controlled (e, g) pairs: sign/value match hand-computed LCM formula; g_penalty matches a finite-difference derivative of the scalar objective; gradient → 0 at exact degeneracy (e_lower == e_upper, any gradients equal); behaviour of sigma/alpha; MECPObjective.__call__/calc_new shape & units contract; raw edelta fallback.
• Paired-engine mechanics: a paired excited/excited call returns distinct finite energies and (natoms,3) gradients for both channels from one SCF+ADC; calc_new returns both; both channels match adcc.nuclear_gradient(excitation_i) independently.
• Distinctness guard: feed a fake reordered/near-degenerate candidate set and assert the two selected roots stay distinct; assert last_trackings carries per-channel diagnostics; assert unresolvable collapse raises.
• Target normalisation: MECI states=(i,j) and MECP ground MP2 + excited paths; MP2-only guard on the ground side; state_index out of range raises; run_adc_kwargs forwarded with native names.
• Reuse scanner singleton regressions (coordinate shape, invalid follow, unconverged SCF) where the paired scanner shares the plumbing.

End-to-end (functionality_geomopt_mecp_test.py, gated on find_spec("geometric") + PySCF, slow):

• Twisted ethylene MECI smoke: short restricted SCF near the 90°-twisted geometry, n_singlets>=2, paired (0,1), maxsteps capped 3–5, convergence_grms/gmax mirrors the existing geomopt tests; assert the energy gap between the two tracked surfaces decreases across steps and both gradients stay finite. Fixture authored fresh (no ethylene fixture exists today).

Acceptance criteria (mapped to issue checklist)

• Paired scanner returns both (energy, gradient) pairs — item 1.
• Excited/excited (MECI) from one SCF + one ADC ✓
• Ground/excited (MECP) from one SCF + LazyMp + ADC ✓ (MP2 enforced)
• Stateful tracking of both states with a distinctness guard keeping them off the same root ✓
• Penalty objective helper combining both into one (energy, gradient) for geomeTRIC, no derivative couplings ✓
• Unit tests for penalty gradient with controlled energies/gradients ✓
• Unit tests that paired scanner returns distinct energies/gradients on both channels ✓
• At least one end-to-end MECI smoke (slow/skipped) + fast unit coverage ✓
• Docs subsection + example ✓

Risks

• Tracker collapse at the seam (highest): two overlap trackers can both lock the same root near a true CI. Mitigation: joint distinct selection + energy-gap tie-break + per-channel TrackingResult diagnostics + tracking_min_score/tracking_min_gap; test the guard with synthetic near-degenerate candidate sets.
• Backward compat of NuclearGradientScanner: the plumbing refactor must leave its passing tests green — guarded by geomopt_scanner_test.py.
• Slow CI: cap maxsteps, reuse the skipif find_spec("geometric") + --mode full conventions; pair the real optimization with fast pure-unit coverage.
• geomeTRIC contract: the MECPObjective must emit exactly {energy: float, gradient: (3N,) Bohr} and (energy, (N,3)) for the as_pyscf_method bridge or the optimizer diverges silently — covered by the shape/units unit tests.
• API fit: states=(i,j) MECI vs lower="mp2", upper=0 MECP — pinned below (two normalised forms); resolved as a single PairedStateGradientScanner accepting states=(i,j) for MECI and lower=/upper= target specs for MECP.

---

Plan created by mach6

[maxscheurer]

Progress Update

Implemented the paired-state MECP/MECI penalty scanner from the plan.

• adcc/gradients/paired_scanner.py — PairedStateGradientScanner subclasses NuclearGradientScanner (no parent refactor); one SCF + one ADC + two nuclear_gradient calls; MECI (states=(i,j)) and MECP (lower="mp2", upper=k, MP2 enforced) paths; joint distinct root selection with a distinctness guard that keeps the two tracked roots apart near the seam and raises on collapse; plural diagnostics (last_trackings, last_excitations, ...) in slot order.
• adcc/gradients/mecp.py — pure numpy mecp_penalty (Levine–Coe–Martínez smoothed, n_states2=1 for the two-state pair) + MECPObjective driving geomeTRIC via the existing as_pyscf_method bridge with no derivative couplings. geomeTRIC stays optional.
• adcc/gradients/__init__.py, adcc/__init__.py — export PairedStateGradientScanner, PairedExcitedStateTarget, PairedGroundExcitedStateTarget, mecp_penalty, MECPObjective.
• adcc/tests/geomopt_mecp_test.py — 34 fast PySCF-only unit tests: finite-difference gradient consistency, hand-computed LCM value, degeneracy collapse to the average surface, raw alpha=0 mode, parameter validation, paired-engine mechanics (distinct energies/gradients, MECI + MECP), distinctness guard with synthetic collapse/too-few-states, normalisation/validation, and a direct cross-check against geomeTRIC's live ConicalIntersection engine (atol 1e-13) — which caught an early n_states2=4-vs-1 bug.
• adcc/tests/functionality_geomopt_mecp_test.py — 3 end-to-end MECI/MECP smoke tests (PySCF+geomeTRIC, maxsteps capped) on a fresh twisted-ethylene fixture.
• docs/gradients.rst — "Minimum-energy crossing points" subsection; examples/optimization/pyscf_adcc_mecp.py — runnable twisted-ethylene MECI example.

Verified: 66 tests pass across the 4 geomopt modules; flake8/ruff clean.

Open question from the discussion: whether to also offer an optional MECPObjectiveConicalIntersection adapter delegating to geomeTRIC's upstream optimizer path for users who want it — deferred to a follow-up.

Commit: b08c5ed

---

Progress tracked by mach6

[maxscheurer]

Code Review

Multi-agent review of PR #231 (MECP/MECI paired-state penalty scanner, closes #229). Five agents ran: code-reviewer, error-auditor, test-reviewer, completeness-checker, simplifier. All findings below have confidence ≥ 80.

Critical

(none)

Important

1. [error-auditor] mecp_penalty raises ZeroDivisionError at exact degeneracy when alpha == 0 — adcc/gradients/mecp.py:96-103. The "raw energy-difference fallback mode" (alpha=0) divides by denom = e_dif + alpha. At exact degeneracy (e_dif == 0), e_pen = sigma*e_dif*e_dif/(denom*1) is 0.0/0.0; because e_dif is a Python float (coerced via float(...)), this raises an unguarded ZeroDivisionError mid-optimization rather than reducing to the average surface. The docstring claims unconditionally "At exact degeneracy the penalty and its gradient vanish, leaving E_obj = E_avg..." — this is only true for alpha > 0. Fix: short-circuit if e_dif == 0.0 or (alpha == 0.0 and abs(e_dif) < ~1e-300): return float(e_avg), g_avg, or raise an explicit ValueError with context, and scope the docstring claim to alpha > 0 (or document the alpha=0 singular point as user-error).

2. [test-reviewer] MECP excited-root overlap-tracking branch is untested — adcc/gradients/paired_scanner.py (_select_excited_root). The overlap branch (previous descriptor present, argsort of scores, TrackingResult, threshold enforcement) never executes: every MECP test uses follow="index" and only the first geometry call. Add a two-call overlap test (seeded scanner, displaced geometry, follow="overlap") asserting last_trackings[1] is populated and index/switched are correct.

3. [test-reviewer] tracking_min_score / tracking_min_gap enforcement untested for the paired scanner — adcc/gradients/paired_scanner.py (_select_excitations, _select_excited_root). The paired scanner duplicates the single-surface threshold guards, but the paired test file contains no tracking_min_score/tracking_min_gap references. Add paired analogues of test_overlap_tracking_below_min_score_raises / test_overlap_tracking_ambiguous_gap_raises for both the MECI joint path and the MECP single-root path.

4. [test-reviewer + error-auditor] tracking_min_gap is silently inert at n_states == 2 — adcc/gradients/paired_scanner.py:405-419 (also flagged by code-reviewer). In _select_excitations, when n == 2 (minimum accepted by the MECI form), remaining is empty for both slots and gap is set to float("inf"), so tracking_min_gap can never trigger regardless of the user setting. The inline comment "the distinctness guard means there is always at least one other candidate here" is wrong at n == 2. Fix: correct the comment and/or document that tracking_min_gap is only effective when n_singlets >= 3; optionally set gap = 0.0 for the empty-remaining case so the guard fires by default. The minimal-2-state config (MECI example uses n_singlets=4, tests use 3) means this is a latent rather than active failure.

Suggestions

5. [test-reviewer] Singular point (alpha=0, e_dif=0) is untested — adcc/gradients/mecp.py (mecp_penalty). test_penalty_alpha_zero_is_raw_squared_difference_mode uses e_dif=0.2; degeneracy tests use the default alpha=0.025. Add a test asserting the actual behavior at (alpha=0, e_dif=0) regardless of which fix-in-1 is chosen (document NaN/raise, or the guarded average-surface return).

6. [test-reviewer] End-to-end smoke does not assert the gap/objective decreases — adcc/tests/functionality_geomopt_mecp_test.py. The MECI smoke collects seen_gaps but only asserts len >= 2 and finiteness; the MECP smoke only checks the final pair is finite/ordered. Add np.isfinite(gradient) per step inside energy_and_gradient and assert the final gap ≤ initial gap (or non-increasing objective) — current assertions amount to "optimizer didn't crash", not "driving toward the seam".

7. [error-auditor] Non-finite energies/gradients flow silently through the pipeline — adcc/gradients/mecp.py:88-103, paired_scanner.py:173-182. No np.isfinite check at any boundary; np.argsort([nan, 1.0]) == [1, 0] silently reorders lower/upper slot assignment, and NaN propagates to MECPObjective.calc_new's dict. A borderline ADC solve near a true seam is the most likely origin of inf/NaN residuals. Fix: assert finite (e_obj, g_obj) at the boundary with a contextual RuntimeError.

8. [test-reviewer] last_pair_order energy-sorting reorder branch untested — adcc/gradients/paired_scanner.py (__call__). In every existing paired test slot order already equals energy order (root 0 < root 1, MP2 ground < excited), so the last_pair_order == (1, 0) branch and the slot-vs-sorted distinction on last_energies/last_gradients are never exercised. Add a test starting where the excited root is below the ground, or inject a synthetic energy-inverted pair.

9. [test-reviewer] last_trackings per-channel gap/switched/previous_index not asserted — adcc/tests/geomopt_mecp_test.py::test_paired_excited_overlap_tracking_seeds_and_tracks_both_slots. Asserts slot0.index != slot1.index and best_score ≈ 1.0 only; the diagnostics that would catch a silent paired root flip (the distinctness guard's raison d'être) are unchecked.

10. [test-reviewer] MECP lower int / GroundStateTarget / bad-type / out-of-range upper branches untested — adcc/gradients/paired_scanner.py (_normalise_paired_target, _select_excited_root). Only the str lower branch and the upper bad-type branch are tested. Add tests for lower=GroundStateTarget(level=2), lower=2 (int), bad-type lower=1.5, and MECP upper out-of-range.

11. [test-reviewer] Unresolvable-collapse raise paths in the paired guard untested — adcc/gradients/paired_scanner.py (_select_excitations). The one-free-slot pick is None raise and the states=(0,0) pure-index collapse raise are untested; test_distinctness_guard_keeps_slots_apart_under_collapse exercises the both-track path and test_distinctness_guard_raises_when_too_few_states covers the <2 excitations raise.

12. [test-reviewer] Misleading test name — adcc/tests/geomopt_mecp_test.py::test_penalty_alpha_zero_is_raw_squared_difference_mode. With alpha=0 the penalty is linear in E_dif (sigma*e_dif), not squared. Rename to e.g. test_penalty_alpha_zero_is_raw_energy_difference_mode.

13. [test-reviewer] Joint tie-break toward smaller excitation-energy gap untested — adcc/gradients/paired_scanner.py (_select_excitations, joint double loop). The tie = abs(total-best_score) <= eps and tie_gap < best_tie_gap branch (bias toward the seam) never triggers under the existing tests. Add a synthetic candidate set with two distinct pairs tied on combined overlap but differing in excitation-energy gap.

14. [simplifier] Collapse the 0- and 1-free-slot branches in _select_excitations — adcc/gradients/paired_scanner.py:438-457. A single for s in free_slots loop (empty for 0 free slots, one iteration for 1) is behavior-identical to the two-branch if len(free_slots) == 0 / elif == 1 structure and removes the separate other_idx/chosen[other] plumbing. (Confidence 80.)

15. [simplifier] Remove dead used/available filtering in the joint branch — adcc/gradients/paired_scanner.py:463-475. The joint-branch is only reached when both slots have score vectors, so both fixed entries are None, used is always empty, and len(available) < 2 is unreachable (top-of-method n >= 2 guard already holds). Iterate range(n) directly; behavior unchanged. (Confidence 92.)

16. [simplifier] Replace energies append-loop with a comprehension — adcc/gradients/paired_scanner.py:232-237. Matches the comprehension style already used directly above for grad_results/grads; behavior identical. (Confidence 85.)

17. [simplifier] Drop duplicate np.asarray(g_obj) in MECPObjective.__call__ — adcc/gradients/mecp.py:191-192. Compute g_obj = np.asarray(g_obj) once and reuse for last_gradient and the return value. (Confidence 80.)

Strengths

• Correct and idiomatic throughout. The MECP energy/gradient level-consistency concern is a non-issue: _check_ground_state_level pins level==2, nuclear_gradient always dispatches to MP2 for any LazyMp, and the LazyMp is constructed without a level arg (defaults to MP2) — energy and gradient use the same MP2 amplitudes.
• Distinctness guard (MECI joint maximisation over distinct ordered pairs + energy-gap tie-break toward the seam + explicit raise on collapse; MECP ground surface uniquely defined by the SCF reference) is correct. states=(0,0) in index mode raises correctly.
• LCM penalty math is correct: e_pen and g_pen_scale match the analytic derivative; _N_STATES2=1 matches geomeTRIC's pair-count normalization; the geomeTRIC ConicalIntersection oracle cross-check passes at atol=1e-13.
• All 34 fast unit tests pass; docs and example are non-stub and accurate; exports mirror the single-surface scanner placement.
• Gating conventions (skipif PySCF missing; slow/gated find_spec("geometric") end-to-end smoke) match geomopt_scanner_test.py.
• completeness-checker independently verified all 9 issue-MECP/MECI optimization: penalty-function paired-state scanner (no derivative couplings) #229 acceptance criteria are met, and the 3 open questions resolved (LCM default; joint distinctness strategy; MECPObjective exposure rather than the geomeTRIC built-in driver).

Agents run: code-reviewer, error-auditor, test-reviewer, completeness-checker, simplifier

---

Reviewed by mach6

[maxscheurer]

Review Assessment

Review comment: #231 (comment)

Classifications

Finding	Classification	Reasoning
1. mecp_penalty ZeroDivisionError at alpha=0, exact degeneracy	genuine	Confirmed at mecp.py:96-103. denom = e_dif + alpha; with e_dif == 0.0 and alpha == 0.0 it's a Python-float 0.0/0.0 → unguarded ZeroDivisionError mid-optimization. The docstring's unconditional "vanishes by construction" claim is only true for alpha > 0.
2. MECP excited-root overlap-tracking branch untested	genuine	Confirmed. Every MECP test uses follow="index" and the first call only; _select_excited_root's overlap branch (scores argsort, TrackingResult, thresholds) never executes. This is new testable code added by the PR.
3. tracking_min_score/tracking_min_gap enforcement untested for paired scanner	genuine	Confirmed. No references to either threshold in the paired test file, while the paired scanner duplicates both raise sites.
4. tracking_min_gap silently inert at n_states == 2	genuine	Confirmed at paired_scanner.py:405-419. With n == 2, remaining == [] → gap = inf → guard never fires; the inline comment ("always at least one other candidate") is wrong. Latent (example uses n_singlets=4, tests 3) but a real contract break.
5. Singular point (alpha=0, e_dif=0) untested	genuine	Tied to finding 1. Tests cover degeneracy only with alpha=0.025 and alpha=0 only with nonzero e_dif. The crashing point is uncovered; the test must assert whichever behaviour finding 1's fix lands on.
6. End-to-end smoke doesn't assert gap/objective decreases	nitpick	The assertion loosened is real but the smoke tests intentionally do tiny 2-3 step runs; a strict monotonic-decrease claim may be brittle over so few geomeTRIC steps. Worth tightening to "final gap ≤ initial gap" + per-step np.isfinite(gradient), but not a correctness blocker.
7. Non-finite energies/gradients flow silently through pipeline	genuine	Confirmed. No np.isfinite checks in __call__ or the objective; np.argsort([nan, 1.0]) == [1, 0] silently reorders lower/upper, then NaN propagates into MECPObjective.calc_new. A borderline-but-valid ADC solve near the seam is the most likely origin.
8. last_pair_order reorder branch untested	genuine	Confirmed. All paired tests have slot order == energy order (ground < excited, root 0 < root 1); the last_pair_order == (1, 0) branch and the slot-vs-sorted distinction are unexercised.
9. last_trackings per-channel gap/switched/previous_index not asserted	genuine	Confirmed. The overlap test asserts distinct indices and best_score only; the diagnostics that would catch a silent paired root flip are unchecked.
10. MECP lower int/GroundStateTarget/bad-type/out-of-range upper branches untested	genuine	Confirmed. Only lower="mp3" (str) and non-int upper type are tested; the int, GroundStateTarget, bad-lower-type, and out-of-range MECP upper paths are uncovered.
11. Unresolvable-collapse raise paths (one-free-slot pick is None; states=(0,0) untested	genuine (partial)	The states=(0,0) positional distinctness-guard raise is genuinely untested. The one-free-slot pick is None path appears unreachable (with len(excitations) >= 2 and one fixed index, at least one different candidate always exists) — the assessor flags it as dead-second-defence code; worth either adding a test for the reachable states=(0,0) case or documenting/removing the unreachable branch.
12. Misleading test name ..._raw_squared_difference_mode	nitpick	Confirmed but harmless. Name says "squared"; alpha=0 gives a linear sigma*e_dif. Body/assertions are correct. Rename at convenience.
13. Joint tie-break toward smaller excitation-energy gap untested	genuine	Confirmed. No test constructs equal combined overlap scores to trigger the tie_gap comparison in the both-slots-track joint loop.
14. Collapse 0/1-free-slot branches into one loop	nitpick	Behavior-preserving but pure style; current branches are readable and correct.
15. Remove dead used/available filtering in joint branch	nitpick	Confirmed redundant (both fixed are None in the else branch, so used is empty and len(available) < 2 is unreachable after the top-of-method n >= 2 guard). Cleanup only.
16. Replace energies append-loop with a comprehension	nitpick	Pure style; behavior identical.
17. Drop duplicate np.asarray(g_obj) in MECPObjective.__call__	nitpick	Confirmed redundant; harmless.

Counts: 12 genuine, 6 nitpicks, 0 false positives, 0 deferred.

Action Plan

1. [finding 1, highest priority] Fix mecp_penalty at alpha=0, e_dif=0: short-circuit to return float(e_avg), g_avg (the documented limit) OR raise an explicit ValueError with context; scope the docstring's "vanishes by construction" claim to alpha > 0.
2. [finding 5] Add a regression test asserting the actual behaviour at (alpha=0, e_dif=0) — whichever path finding 1's fix takes.
3. [finding 4] Decide paired tracking_min_gap semantics for exactly two states: either set gap = 0.0 when remaining is empty (so the guard fires by default), raise an explicit error, or document that tracking_min_gap requires n_singlets >= 3. Correct the wrong inline comment.
4. [finding 7] Add np.isfinite validation for (energies, gradients) in PairedStateGradientScanner.__call__ before the argsort, raising a contextual RuntimeError so NaN doesn't silently reorder slots or propagate into the optimizer.
5. [finding 3] Add paired-scanner threshold tests for tracking_min_score and tracking_min_gap on both the MECI joint path and the MECP single-root path.
6. [finding 2] Add a two-call MECP overlap test (seeded scanner, displaced geometry, follow="overlap") asserting last_trackings[1] is populated and index/switched are correct.
7. [finding 8] Add a test starting where the excited root is below the ground (or an injected synthetic energy-inverted pair) asserting last_pair_order == (1, 0) and energy-sorted return while last_energies stays in slot order.
8. [finding 9] Assert gap, previous_index, and switched on both paired TrackingResults in the existing overlap test.
9. [finding 10] Add MECP target-normalisation tests for lower=GroundStateTarget(level=2), lower=2 (int), bad lower type, and out-of-range MECP upper.
10. [finding 11] Add a states=(0,0) second-call test asserting the distinctness-guard raise. Separately, document or remove the unreachable one-free-slot pick is None branch.
11. [finding 13] Add a synthetic candidate set with two distinct pairs tied on combined overlap but differing in excitation-energy gap, asserting the smaller-gap pair is chosen.
12. [optional] Findings 12, 14, 15, 16, 17 (rename + simplifier cleanups) — address opportunistically; none are merge blockers.

No false positives; no items deferred (all genuine findings concern testable code added by this PR and should ship with it). No tracking issues needed.

---

Assessment by mach6

[maxscheurer]

Progress Update

Addressed the review/assessment findings and converted the example to the S0/S1 MECP.

Code fixes

• mecp_penalty: short-circuits the alpha=0 / exact-degeneracy singular point to the average surface (was an unguarded ZeroDivisionError from 0.0/0.0); covers the sub-DBL_MIN underflow too. Raises a contextual RuntimeError on non-finite objective output. Docstring "vanishes by construction" claim now scoped to alpha > 0 with the alpha=0 limit explained. (assessment items 1 & 7, penalty side)
• PairedStateGradientScanner.__call__: validates np.isfinite of both surface energies and gradients before the energy argsort, which was silently reordering lower/upper slots on NaN. (item 7, scanner side)
• tracking_min_gap: no longer silently inert when only two excited states are computed. On both the MECI joint path (_select_excitations) and the MECP single-root path (_select_excited_root) it now raises a clear RuntimeError when tracking_min_gap > 0 and a per-channel gap is unmeasurable; the default (== 0) still passes and reports gap = inf. The wrong inline comment ("always at least one other candidate") is corrected. (item 4)

Tests added (+22 fast units in geomopt_mecp_test.py, now 56 total)

• Exact-degeneracy alpha=0 collapse to average; near-degeneracy finiteness; non-finite energy/gradient rejection (items 5 & 7).
• Paired tracking_min_score / tracking_min_gap raises for both MECI and MECP (item 3).
• Two-state tracking_min_gap raises-with-context vs inert-at-default with gap == inf (item 4).
• MECP excited-root overlap tracking across two calls (item 2) and per-channel gap/previous_index/switched diagnostics (item 9).
• Energy-sort reorder via synthetic surface injection: last_pair_order == (1, 0) with last_energies kept in slot order, plus the already-sorted default (item 8).
• MECP lower as GroundStateTarget/int, bad-type lower, out-of-range MECP upper (item 10).
• states=(0, 0) index-mode collapse → distinctness-guard raise (item 11).
• Joint seam-bias tie-break picks the smaller excitation-energy-gap pair from a controlled-score synthetic candidate set (item 13).

Example (examples/optimization/pyscf_adcc_mecp.py): now locates the S0/S1 MECP (MP2 ground S0 + first tracked excited singlet S1) instead of the S1/S2 MECI — lower="mp2", upper=0, n_singlets=1, docstring and printed labels updated; verified the construction runs and returns finite S0/S1 surfaces.

Verification: 120 geomopt tests pass (56 mecp fast + 27 scanner + 34 contraction + 3 end-to-end geomeTRIC MECI/MECP smoke); flake8 + ruff clean; example and modules compile.

Not addressed this batch (nitpicks, follow-up): assessment items 14–17 (simplifier cleanups), item 12 (test rename), and the likely-unreachable one-free-slot pick is None defensive raise in _select_excitations.

Commit: b417787

---

Progress tracked by mach6

[maxscheurer]

Progress Update

Switched the worked example to a spin-flip ADC(2) triplet MECI, which avoids the convergence trouble of the ground/excited MECP near the diradical seam.

• examples/optimization/pyscf_adcc_mecp.py: triplet unrestricted HF reference (scf.UHF, spin=2) on 90°-twisted ethylene + spin-flip ADC(2) (n_spin_flip=4), optimising the MECI between the first two spin-flip states (states=(0, 1)). The first spin-flip state recovers the singlet ground state S0 and the next one the first excited singlet S1, so the S0/S1 crossing is reached as a single-solve excited/excited MECI -- no LazyMp ground surface, no derivative couplings. Verified the geomeTRIC penalty optimisation runs every step cleanly (the prior ground/excited MECP crashed); S0/S1 are nearly degenerate at the triplet geometry (gap ≈ 2.6 mEh).
• docs/gradients.rst: added a "Spin-flip ADC for S0/S1 conical intersections" note in the MECP/MECI subsection pointing users at this setup and the example.

The end-to-end functionality smoke tests are intentionally left on their existing MECI/MECP forms this batch (independent of the example); a follow-up can switch them to SF-ADC if we want full coverage of the new setup.

Fast 56-test geomopt_mecp_test suite still green; flake8/ruff clean.

Commit: 73e5f6c

---

Progress tracked by mach6

[maxscheurer]

Progress Update

Root-caused the SF-ADC example divergence (final 9e-2 Eh gap) and fixed it. The gap is now driven to ~5.9e-4 Eh (~0.016 eV) — a proper converged MECI.

Diagnosis. The example used follow="overlap" for the MECI pair. Density-overlap tracking is designed for single-surface optimization (stay on the same physical state as geometry changes). For a MECI pair it does the wrong thing: the seam is defined by a degeneracy of the energy-ordered two lowest roots, not by "the same two state characters you started with." Once S0/S1 densities blur at the seam, the overlap tracker could not tell them apart and flipped slot 1 onto root 2 (a higher valence state), so the optimizer started chasing S1/S2 and the gap exploded.

Cross-check against other codes.

• geomeTRIC's ConicalIntersection engine (engine.py:2129) does no identity tracking itself: it calls one sub-engine per state, re-sorts all returned energies by value (np.argsort), and penalises the energy-ordered gap with the same LCM penalty MECPObjective already uses. Its sub-engines are expected to return "roots 0 and 1" — i.e. the adiabatically lowest roots.
• Levine–Coe–Martínez (the paper + Q-Chem's PENALTY_FUNCTION) is identity-agnostic by the same construction: minimise the gap of the energy-ordered pair, driven harder with larger σ.
• Q-Chem's STATE_FOLLOW (overlap of attachment/detachment densities, with an energy-window fallback) is an optional separate feature for single-state excited-state geometry optimisation — explicitly not the MECP path. It is the analogue of our scanner's follow="overlap", and the manual itself only nudges toward it for hard spin-contamination cases.

The fix (and the right mental model).

• follow="index" for a MECI pair — positional roots (0, 1), the scanner energy-sorts them and MECPObjective penalises the gap. This matches the contract geomeTRIC's penalty driver assumes from its sub-engines.
• follow="overlap" stays the right tool for single-surface NuclearGradientScanner (stay on "the bright ππ*" across a twist); it just should not be used to label a crossing pair.
• Also raised the penal strength σ to 200 (tighter seam enforcement) and tightened geomeTRIC grms=3e-5 / gmax=1e-4; kept follow="index".

Verified end-to-end on twisted-ethylene SF-ADC(2): converges cleanly to gap ≈ 5.9e-4 Eh with no root flip across all ~40 geomeTRIC steps. A matching note ("use follow="index" for a MECI pair, not follow="overlap"") is added to the MECP/MECI docs subsection.

This is a usage fix, not a scanner bug: the single-surface overlap tracker is correctly implemented for its purpose; the error was applying it where the adiabatic ordering should be preserved instead. No scanner code changed.

Commit: 6b93dd9

---

Progress tracked by mach6

[maxscheurer]

Progress Update

Implemented option B: the scanner now warns loudly when follow="overlap" is used for a MECI pair, so the footgun that left the SF-ADC example at a 90 mEh gap is flagged at construction time rather than discovered as a diverged optimisation.

• paired_scanner.py: in __init__, if follow == "overlap" and the paired target is a PairedExcitedStateTarget (excited/excited MECI), emit a UserWarning explaining that density-overlap tracking fights the adiabatic energy-ordering that defines the seam and can flip a slot onto a higher root, and pointing at follow="index" (geomeTRIC's ConicalIntersection sub-engine contract). The warning is specific to the MECI pair form — follow="overlap" stays a legitimate choice for ground/excited MECP (only one root is tracked) and for single-surface NuclearGradientScanner, so it is not raised there.
• geomopt_mecp_test.py: three new tests — assert the warning fires for a MECI pair with follow="overlap", asserts it does not fire for MECP, and asserts follow="index" (the recommended MECI setting) is silent. The eight existing tests that deliberately exercise the overlap-tracking mechanism on a MECI pair are marked @pytest.mark.filterwarnings("ignore::UserWarning") since they are testing the (still-valid) machinery, not the recommendation.

Why warning, not raising: a hard ValueError would forbid users who genuinely know what they are doing (e.g., well-separated roots, or wanting to inspect overlap diagnostics near the seam). The warning surfaces the recommendation without removing the capability; the docs note (previous commit) explains the "why".

All 59 fast geomopt_mecp_test tests pass; geomopt_scanner_test still green (27); flake8/ruff clean on all changed .py files. The option C energy-continuity tie-break (a faithful reimplementation of Q-Chem's STATE_FOLLOW fallback for when a user does want overlap tracking on a pair) remains as a sensible follow-up.

Commit: 70f19e0

---

Progress tracked by mach6

[maxscheurer]

Code Review (round 2)

Second multi-agent review of PR #231 (MECP/MECI paired-state penalty scanner, closes #229), run after the round-1 fixes (commits b417787, 73e5f6c, 6b93dd9, 70f19e0). Five agents ran: code-reviewer, error-auditor, test-reviewer, completeness-checker, simplifier. All findings below have confidence ≥ 80.

Round-1 findings — confirmed resolved

• Prior finding 1 (mecp_penalty ZeroDivisionError at alpha=0, exact degeneracy): resolved — short-circuits to the average surface; docstring claim scoped.
• Prior finding 4 (tracking_min_gap silently inert at n==2): resolved on both MECI and MECP paths — now raises with context when tracking_min_gap > 0 and a per-channel gap is unmeasurable.
• Prior finding 7 (np.argsort silently reorders slots on NaN): resolved — __call__ validates finite energy + gradient per slot before the sort; mecp_penalty raises on non-finite objective output.
• MECI follow="overlap" footgun: resolved — __init__ emits a specific UserWarning for an excited/excited MECI pair; 8 machinery tests filter it, 3 dedicated tests pin the warning scope.

Critical

(none)

Important

1. [code-reviewer] Docs MECI snippet still shows follow="overlap", contradicting the PR's own recommendation — docs/gradients.rst, the states=(0, 1) ... follow="overlap" MECI code block in the "Minimum-energy crossing points" section. The .. note:: printed immediately below now says "For a MECI pair use follow="index" ... not follow="overlap"", and commits 6b93dd9/70f19e0 established follow="overlap" is a MECI footgun that diverged the SF-ADC example to a 90 mEh gap. The most copy-pasteable item on the page instructs the user to do exactly what the same page warns against. Fix: change this snippet's follow="overlap" to follow="index" (leave the MECP snippet's follow="overlap" — only one root is tracked there).

2. [error-auditor] alpha=0 exact-degeneracy guard silently produces a wrong gradient — discontinuous objective at the seam — adcc/gradients/mecp.py (guard ~L141-147; docstring ~L88-96). For alpha == 0, e_dif != 0 the formula gives g_pen_scale = sigma * (E_dif**2 + 0) / (E_dif**2) = sigma (a constant, independent of e_dif), so the correct continuous extension of the gradient at the seam is G_avg + sigma * G_dif — the penalty gradient does not vanish for alpha = 0; only the energy does. The guard, however, sets g_pen_scale = 0.0 at e_dif == 0.0 for any alpha, so exactly at the seam the gradient jumps from G_avg + sigma*G_dif (just off-seam) to G_avg (on-seam), dropping precisely the term that pins the optimiser to the crossing. This is pinned by the tests as intended behaviour (test_penalty_alpha_zero_is_raw_squared_difference_mode asserts +2*sigma*G_dif off-seam; test_penalty_alpha_zero_exact_degeneracy_collapses_to_average asserts G_avg at-seam), and the docstring's claim "At exact degeneracy the penalty and its gradient vanish … includes the alpha == 0 raw energy-difference mode" is false for the gradient at alpha = 0. Reachable at symmetry-enforced crossings (degenerate triplets, Jahn–Teller / high-symmetry MECIs — the PR's own example history was a SF-ADC triplet MECI). The default alpha=0.025 avoids this (the guard is mathematically correct for alpha > 0); impact is confined to users who opt into the documented alpha=0 "raw energy-difference mode". Fix: in the alpha == 0 near-seam branch keep e_pen = 0.0 but set g_pen_scale = sigma (the genuine alpha→0 limit) so the extension is continuous; scope the docstring to alpha > 0 and document the alpha=0 gradient limit as sigma * G_dif. Alternatively, raise on alpha=0 exact-degeneracy rather than silently returning the wrong Jacobian.

Suggestions

3. [test-reviewer] End-to-end smoke asserts "didn't crash", not "driving toward the seam" (prior finding 6, still open) — adcc/tests/functionality_geomopt_mecp_test.py. The MECI smoke collects seen_gaps but only asserts len >= 2 and finiteness; the MECP smoke doesn't collect gaps at all. Prior finding 6 was downgraded to a nitpick and not addressed — but the PR's own example did diverge to a 90 mEh gap under follow="overlap", and a non-directional smoke would not have caught that regression. Fix: assert seen_gaps[-1] <= seen_gaps[0] + 1e-9 (or final objective ≤ initial) plus per-step np.isfinite(gradient).

4. [code-reviewer + test-reviewer] End-to-end MECI smoke runs the warned-against follow="overlap" configuration unsilenced — adcc/tests/functionality_geomopt_mecp_test.py::test_paired_scanner_objective_drives_meci_optimization (~L90-99). Constructs a MECI pair with follow="overlap", which since 70f19e0 emits the UserWarning that says overlap tracking "can flip a slot onto a higher root, diverging the penalty optimisation." The 8 fast machinery tests were marked filterwarnings("ignore::UserWarning"); this slow smoke was not, so it spews the warning every run and exercises the exact footgun — only maxsteps=3 keeps it from blowing up. Fix: switch this smoke to follow="index" (consistent with the example and the docs note) or silence the warning and assert the two tracked roots stay distinct at every recorded step, not just the final one.

5. [test-reviewer] MECP single-state tracking_min_gap > 0 raise path untested — adcc/gradients/paired_scanner.py (_select_excited_root, elif self.tracking_min_gap > 0: raise). The MECP threshold tests all use n_singlets=3 and hit the len(order) > 1 path; the n_singlets=1 raise (the MECP example's config) is never executed. Fix: a MECP test with n_singlets=1, follow="overlap", tracking_min_gap > 0, two calls, asserting the "at least two computed excited states" raise.

6. [test-reviewer] last_trackings.gap assertion allows the silent-collapse sentinel to pass — adcc/tests/geomopt_mecp_test.py::test_paired_overlap_tracking_records_per_channel_diagnostics. Asserts np.isfinite(slot.gap) and slot.gap >= 0.0; but gap == 0.0 is the dangerous value (best ≈ second-best), and with n_singlets=3 at one geometry the realised gap is ≈ 1.0. Fix: tighten to assert slot.gap > 0.5 (or pytest.approx) so the diagnostic is shown to actually discriminate.

7. [error-auditor, low] The underflow guard comment over-states coverage — adcc/gradients/mecp.py. The secondary clause abs(e_dif) < 1e-300 does not cover the actual underflow regime (e_dif**2 flushes to 0 for |e_dif| < ~7e-162), so e_dif in [1e-300, 7e-162) with alpha == 0 raises an uncaught ZeroDivisionError. Physically unreachable (~1e-200 Eh gaps) and a loud crash rather than a silent one, so low severity — but the comment is off by ~162 orders of magnitude and the existing 1e-310 test sits inside the guard window, not in the gap.

8. [test-reviewer, trivial rename] Misleading test name ..._raw_squared_difference_mode — alpha=0 gives a linear sigma*E_dif, not squared. Body/assertions are correct; only the name lags. Rename at convenience.

9. [simplifier] Collapse the 0- and 1-free-slot branches into one greedy loop — adcc/gradients/paired_scanner.py (_select_excitations). The 0-free case is a no-op over the loop body; int(o) not in chosen is equivalent to int(o) != other_idx. Behavior-preserving, ~3 fewer lines. (Still valid from round 1.)

10. [simplifier] Dead used/available filtering in the joint branch — adcc/gradients/paired_scanner.py (both-track else). fixed is [None, None] there, so used is always empty and the len(available) < 2 guard is unreachable after the top-of-method n >= 2 check. Iterate range(n) directly. (Still valid from round 1.)

11. [simplifier] Replace the energies append-loop with a comprehension — adcc/gradients/paired_scanner.py (__call__). Direct conditional-expression comprehension over targets; behavior identical. (Still valid from round 1.)

12. [simplifier] Drop the double np.asarray(g_obj) in MECPObjective.__call__ — adcc/gradients/mecp.py. Materialise once, store + return the same reference; the two calls already alias the same object today. (Still valid from round 1.)

13. [simplifier, new] Use zip in the finiteness-guard loop in __call__ — adcc/gradients/paired_scanner.py. enumerate(zip(energies, grads)) removes the parallel grads[slot] lookup; behavior identical. (Introduced by the fix commit.)

Strengths

• The penalty math, joint distinct-root selection (ordered-pair maximisation with a seam-biased tie_gap tie-break), the MECP MP2-level enforcement, and the last_pair_order / last_energies slot-vs-sorted discipline are all correct and idiomatic.
• The geomeTRIC oracle cross-check at atol 1e-13 is a strong correctness anchor and already caught an early n_states2 bug.
• The follow="overlap" MECI warning is well-scoped (MECP and single-surface use are left legitimate) and backstopped by clear diagnostics.
• Round-1 fixes were landed precisely against the assessment; the unit suite (59 tests) is green and flake8/ruff clean.
• Documentation and a runnable SF-ADC(2) example are present and the example was driven to a genuinely-converged MECI (~0.016 eV gap).

Completeness

All 9 explicit acceptance criteria from issue #229 are fully implemented and verified against the changed files; all 4 plan deliverables are present and complete. The only completeness gap is low-severity and non-blocking: the plan promised a adcc.paired_nuclear_gradient_scanner(...) convenience alias (and the issue's own design sketch uses it), but it was never added (PairedStateGradientScanner is class-only, consistent with the existing single-surface NuclearGradientScanner which also has no lowercase function). Worth either adding the alias or amending the issue sketch in the docs.

Agents run: code-reviewer, error-auditor, test-reviewer, completeness-checker, simplifier

---

Reviewed by mach6

[maxscheurer]

Review Assessment

Review comment: #231 (comment)

Classifications

Finding	Classification	Reasoning
1. Docs MECI snippet still uses follow="overlap" contradicting its own note	genuine	Confirmed in docs/gradients.rst: the MECI example uses states=(0, 1) with follow="overlap", while the .. note:: immediately below says MECI pairs should use follow="index". Real user-facing contradiction.
2. alpha=0 exact-degeneracy guard drops the penalty gradient (discontinuous objective at seam)	genuine	Confirmed in adcc/gradients/mecp.py. Guard sets g_pen_scale = 0.0 for e_dif == 0.0 regardless of alpha; for alpha=0 the off-seam limit is g_pen_scale = sigma (constant), so the correct continuous gradient at the seam is G_avg + sigma*G_dif. Existing tests assert +sigma*G_dif off-seam and average-only at-seam, pinning the discontinuity as "intended" — but it is mathematically wrong for the documented alpha=0 mode and silently drops the seam-pinning force at the crossing.
3. E2E smoke asserts "didn't crash", not "driving toward seam"	nitpick	Accurate but these are intentionally tiny smoke tests; stronger gap-trend assertions would improve coverage, not a correctness bug.
4. E2E MECI smoke uses warned follow="overlap" unsilenced	nitpick	Confirmed; test hygiene/consistency with the docs note, not a runtime correctness issue.
5. MECP single-state tracking_min_gap > 0 raise untested	nitpick	Confirmed narrow coverage gap (tests use n_singlets=3).
6. last_trackings.gap >= 0.0 assertion too loose	nitpick	Confirmed; gap == 0.0 would pass, but the code is not wrong.
7. Underflow guard comment/coverage wrong for alpha=0 (abs(e_dif) < 1e-300 vs real ~7e-162)	genuine	Confirmed; ZeroDivisionError for e_dif in [1e-300, 7e-162) with alpha=0. Low severity (loud, physically unreachable ~1e-200 Eh) but real.
8. Misleading test name ..._raw_squared_difference_mode	nitpick	Confirmed; rename only.
9-12. Simplifier refactors (collapse free-slot branches; dead used/available filtering; energies append-loop → comprehension; drop double np.asarray(g_obj))	nitpick	All confirmed behavior-preserving cleanups; still valid from round 1.
13. Use zip in finiteness-guard loop	nitpick	Confirmed; cleaner but behavior is correct.
Completeness note: missing paired_nuclear_gradient_scanner lowercase alias	false-positive	The issue sketch used a lowercase helper name, but acceptance criteria require a paired scanner object. PairedStateGradientScanner matches the existing class-based NuclearGradientScanner API pattern (no lowercase alias exists for the original scanner either); consistent convention, non-blocking.

Counts: genuine = 3, nitpicks = 9, false positives = 1, deferred = 0.

Action Plan

1. [finding 2, highest priority] Fix mecp_penalty for alpha == 0 at exact degeneracy: in the alpha == 0 near-seam branch keep e_pen = 0.0 but set g_pen_scale = sigma (the genuine alpha → 0 limit) so the seam extension is continuous and the penalty gradient keeps pinning the crossing — instead of silently returning G_avg. Alternatively, reject the alpha == 0 exact-degeneracy case explicitly rather than dropping the Jacobian.
2. [finding 2] Update the alpha=0 docstring and the two alpha=0 tests so they no longer claim/right-assert that the gradient vanishes at exact degeneracy; scope the "vanishes by construction" claim to alpha > 0 and document the alpha=0 gradient limit as sigma * G_dif.
3. [finding 7] Fix the alpha=0 underflow path — ideally by special-casing the raw mode as e_pen = sigma * e_dif and g_pen_scale = sigma (which also resolves finding 2), or else by widening the guard threshold to ~1e-160. Correct the over-stated comment.
4. [finding 1] Change the MECI docs snippet in docs/gradients.rst from follow="overlap" to follow="index", matching the note below it and the fixed example.
5. [nitpicks 4, 3] Optionally: switch the e2e MECI smoke to follow="index" (or silence + assert distinct at every step), and add a directional seen_gaps[-1] <= seen_gaps[0] assertion.
6. [nitpicks 5, 6] Optionally: add the MECP single-state tracking_min_gap test; tighten the last_trackings.gap assertion to > 0.5.
7. [nitpicks 8-13] Optional cosmetic: rename the alpha=0 test; apply the simplifier cleanups (collapse free-slot branches; remove dead used/available; energies comprehension; drop double np.asarray(g_obj); zip in finiteness guard).

No deferred items need tracking issues.

---

Assessment by mach6

[maxscheurer]

Progress Update

Implemented the four round-2 review/assessment findings (genuine items 1, 2 + folded 7; nitpicks 3, 4).

Penalty math (assessment items 2 & 7 — genuine)

• mecp_penalty (adcc/gradients/mecp.py): the alpha == 0 "raw energy-difference mode" is now evaluated in closed form with no division — E_pen = sigma * E_dif / n_states2 and g_pen_scale = sigma / n_states2 — the genuine continuous extension of the smoothed penalty. The penalty energy still vanishes at exact degeneracy, but the penalty gradient tends to the constant sigma * G_dif (the force that pins the optimiser to the crossing) and is no longer silently dropped at E_dif == 0. The previous short-circuit set g_pen_scale = 0.0 at the seam, which produced a discontinuous objective that threw away the seam-pinning force exactly where it matters — a real bug at symmetry-enforced crossings (degenerate triplets, Jahn–Teller MECIs), reached via the documented alpha=0 mode.
  • Because the raw branch carries no division, the former ZeroDivisionError gap in [1e-300, ~7e-162) (where E_dif**2 flushes to zero in float64) is eliminated entirely, not just the abs(E_dif) < 1e-300 window — so item 7 is resolved by the same change. The over-stated underflow comment is corrected.
  • The alpha > 0 smoothed path is unchanged (denom = E_dif + alpha >= alpha > 0 everywhere); its "penalty and gradient vanish at the seam" guarantee still holds, and the docstring's "vanishes by construction" claim is now correctly scoped to alpha > 0, with the alpha == 0 gradient limit documented.

Tests

• test_penalty_alpha_zero_exact_degeneracy_collapses_to_average was pinning the discontinuity as "intended" (asserted the bare average surface). Renamed to ..._keeps_penalty_gradient and it now asserts the correct limit G_avg + sigma * G_dif at E_dif == 0.
• New test_penalty_alpha_zero_underflow_regime_is_finite (E_dif = 1e-200, squarely in the former crashing gap) — regression guard for item 7.
• The near_degeneracy and off-seam alpha=0 tests' comments were refreshed to reflect the no-division mechanism; the off-seam test (already asserting +sigma * G_dif) passes unchanged.

Docs (item 1 — genuine)

• docs/gradients.rst: the MECI snippet now uses follow="index" (two lowest adiabatic roots, energy-sorted), consistent with the .. note:: immediately below it, the fixed example, and the ConicalIntersection-sub-engine contract. The MECP snippet's follow="overlap" is left as-is (only one excited root is tracked there).

End-to-end smoke (items 3 & 4 — nitpicks)

• functionality_geomopt_mecp_test.py: the MECI smoke switched from follow="overlap" (the warned-against MECI footgun that diverged the SF-ADC example to a 90 mEh gap) to follow="index". The now-inapplicable last_trackings non-None assertions are replaced by a structural-distinctness check (last_pair energy-sorted).
• Both the MECI and MECP smokes now record seen_gaps and assert a per-step finite objective and a directional seen_gaps[-1] <= 2.0 * seen_gaps[0] + 1e-9 anti-divergence bound — so the tests catch a diverging optimisation rather than merely "didn't crash". Verified empirically the MECI optimizer genuinely drives the gap down (0.042 → 0.0064 Eh min over the 3-step run); the 2× factor catches a doubling yet stays robust over the tiny line search.

Verification: 60 fast geomopt_mecp_test tests pass (was 59; +1 underflow regression test), 27 geomopt_scanner_test pass, 3 end-to-end geomeTRIC smokes pass; flake8 + ruff clean on all changed files.

Not addressed this batch (nitpicks, sensible follow-up): assessment item 5 (MECP single-state tracking_min_gap test), item 6 (tighten last_trackings.gap assertion), item 8 (..._raw_squared_difference_mode test rename), items 9-13 (simplifier cleanups).

Commit: 1a63b59

---

Progress tracked by mach6

[maxscheurer]

Progress Update

Reduce example boilerplate: scanners and the MECP objective now accept a PySCF Mole directly, so the def energy_and_gradient(mol_at_step): return scanner(mol_at_step.atom_coords(unit="Bohr")) wrapper (and the docs' lambda m: objective(m.atom_coords(unit="Bohr"))) disappears. Per-step diagnostics are kept as an opt-in callback rather than baked into a hand-rolled wrapper.

Core (backward compatible)

• adcc/gradients/scanner.py: _coords_array now also accepts an object exposing atom_coords(unit="Bohr") (duck-typed, no new dependency); the array path is unchanged. NuclearGradientScanner gained an opt-in step_callback: Optional[Callable] invoked as cb(energy, gradient) at the end of every __call__. PairedStateGradientScanner inherits both via the parent's _coords_array.
• adcc/gradients/mecp.py: MECPObjective got the same step_callback hook, firing after last_pair is populated so a callback can read the seam gap.

The minimum reasonable setup is now literally:

method = as_pyscf_method(mol, scanner)        # ground/excited single-surface
mol_eq = geometric_solver.optimize(method, maxsteps=20)

or for a paired MECI/MECP:

objective = adcc.MECPObjective(paired, sigma=200.0)
method = as_pyscf_method(mol, objective)
mol_ci = geometric_solver.optimize(method, maxsteps=50, ...)

with an optional one-line monitor:

objective.step_callback = lambda e, g: print(f"E_pen={e:.8f}, |g|={np.linalg.norm(g):.4e}")

Examples

• examples/optimization/pyscf_adcc_geoopt.py: the def energy_and_gradient wrapper is gone; body is as_pyscf_method(mol, scanner).
• examples/optimization/pyscf_adcc_geoopt_excited.py: make_energy_gradient factory → report_step(scanner) callback factory assigned to scanner.step_callback; both stages pass as_pyscf_method(mol, scanner) directly. The per-step E / state / ω console output is preserved as the opt-in callback.
• examples/optimization/pyscf_adcc_mecp.py: def energy_and_gradient → objective.step_callback = _print_step; as_pyscf_method(mol, objective) directly. The per-step E_pen / |g| / gap print (which was load-bearing when diagnosing the follow="overlap" divergence) is preserved.

Docs

• docs/gradients.rst: the lambda m: objective(m.atom_coords(unit="Bohr")) becomes as_pyscf_method(mol, objective), and a short note documents the direct-Mole-acceptance plus the optional step_callback hook. The MECI snippet keeps follow="index" from the earlier fix.

Tests (5 new fast units)

• geomopt_scanner_test.py: single-surface test_scanner_accepts_pyscf_mole_directly (Mole⇄coords equivalence) + ..._step_callback_is_invoked_....
• geomopt_mecp_test.py: paired scanner + MECPObjective Mole-acceptance equivalence + test_mecp_objective_step_callback_is_invoked (asserts the callback fires after last_pair is set, with a readable seam gap).

The end-to-end smoke tests (functionality_geomopt_mecp_test.py) were intentionally left with their def energy_and_gradient, since those carry genuine per-step assertion logic (finite checks, the directional seen_gaps anti-divergence bound) — not pure plumbing.

Verification: 92 fast geomopt_* tests pass (was 87; +5), 3 e2e geomeTRIC smokes pass; flake8 + ruff clean on all changed files; examples byte-compile; docs rebuild clean (no Python lambda remains in the snippets — only the Greek λ in the theory math).

Commit: e151640

---

Progress tracked by mach6