Return every atom map per reaction and cluster them into distinct reaction channels

[kfir4444]

Stacked on #1023 — please review that first; this branch's base will collapse once it lands.

Problem

map_reaction returns one atom map per reaction. Several places in the pipeline compute more than one
and discard the rest:

• map_two_species scores every superimposable backbone candidate, then keeps the lowest-RMSD one and
  breaks ties by taking the last tied candidate — an arbitrary choice.
• map_rxn returns on the first template product_dict that succeeds, and only advances to the next one
  when a stage fails.

So when a reaction genuinely has more than one distinct channel, ARC silently picks one of them by a
tie-break, and there is no way to ask how many equivalent paths a channel has.

What this adds

arc/mapping/cluster.py: enumerate every valid atom map, then group them into equivalence classes. Each
class is one chemically distinct reaction channel (one TS to search for); its degeneracy is that channel's
reaction path degeneracy.

Exposed as ARCReaction.atom_map_clusters — lazily computed, cached, never triggered by atom_map, and not
persisted by as_dict. atom_map itself is untouched, so every existing consumer is unaffected.

The equivalence relation

Two maps are equivalent iff sigma_2 = beta . sigma_1 . alpha for some alpha in Aut(R), beta in
Aut(P) — the double cosets. That reduces to a single group. Writing C(sigma) for the set of bonds of R
that break, form, or change order, and R' for R with those changes applied, any valid sigma is an
isomorphism R' -> P. So equal changed-bond sets differ by an element of Aut(P), and:

sigma_1 ~ sigma_2   <=>   C(sigma_1) and C(sigma_2) lie in the same Aut(R) orbit

Aut(P) never appears. Automorphisms are enumerated in-module by colour-pruned backtracking over the
core skeleton (hydrogens folded onto their parent), which avoids the factorial blow-up — isobutane has
1296 full-graph automorphisms and 6 core ones — and needs no new dependency.

Reaction path degeneracy

Degeneracy is not a count of what the enumeration found; that depends on how many template matches RMG
happens to generate, and under-reports. It is the orbit size of the reaction center under the full
automorphism group, which factorises:

degeneracy = |orbit of the collapsed center under Aut_core(R)|
             * product over core atoms p of  n_p! / (n_p - k_p)!

Verified against known values:

reaction	degeneracy
CH4 + OH -> CH3 + H2O	4
C2H6 + OH -> C2H5 + H2O	6
C3H8 + NH2 -> n-C3H7 + NH3	6
C3H8 + NH2 -> i-C3H7 + NH3	2
C2H5Cl -> C2H4 + HCl	3
cyclohexane + OH -> C6H11 + H2O	12

Cyclohexane exercises both halves: 6 equivalent carbons times 2 hydrogens each.

Validity filtering

changed_bonds computes a center for any element-preserving bijection, so scrambled maps would otherwise
open bogus clusters — for the Diels-Alder of butadiene with ethene the correct map changes 2 bonds while
scrambled ones change 10 to 14, and each opened its own cluster. Maps are validated against the family
recipe (ARCReaction.get_expected_changing_bonds over each product_dict's r_label_map), falling back to
a minimal-center heuristic where the recipe is unreadable. On the corpus the two agree.

Test

arc/mapping/cluster_test.py, 60 tests: graph construction, automorphism counts against known group orders
(benzene 12, naphthalene 4, neopentane 24, 2x benzene 288), Kekule bond-order handling, changed-bond sets
and their Aut(P) invariance, canonical keys, the degeneracy factorisation, both validity filters, and an
end-to-end CH4 + OH run.

Expected values are derived independently — group orders from graph theory, degeneracies from known
reaction path degeneracies — not read off the implementation. Each behaviour was mutation-checked: reverting
the code under test makes the corresponding test fail.

arc/mapping/ arc/reaction/ arc/checks/nmd_test.py  ->  217 passed

Corpus effect

testing/errs_no_li.yml, 452 reactions, 14 workers, 200 s each:

ok       330 (73.0%)      1 cluster:  329  (99.7%)
no_maps  120 (26.5%)      2 clusters:   1  (0.3%)

One cluster is overwhelmingly the norm. With reactant and product species both pinned, different
reacting sites generally give different product graphs, so naming the products pins the channel. More than
one cluster requires two Aut(R)-inequivalent centers yielding isomorphic products.

The single multi-cluster case is C10H9-3 <=> C10H9-4 (cyclopentadienyl-cyclopentadiene, |Aut(R)| = 1):
an H shift within one ring, or one that swaps the rings' roles, both give an isomorphic product. Both
centers are predicted by the family recipe, so neither is an enumeration artifact — though whether the
second corresponds to a physically sensible (large-ring) TS is a question for a chemist. map_reaction
resolves that ambiguity today by an arbitrary tie-break.

So with both wells specified the feature mostly confirms degeneracy; at 1-in-330 a multi-cluster result is
signal worth investigating rather than noise. Its sharper use is enumerating channels when products are
not fixed in advance.

Cost

Enumeration is substantially more expensive than return-on-first-success — it sweeps every product_dict,
both directions, and every backbone candidate per fragment. That is why atom_map_clusters is opt-in and
never triggered by atom_map. Caps exist at three levels (MAX_AUTOMORPHISMS, MAX_ENUMERATED_MAPS,
map_rxn_all's max_maps) and all log when hit; a truncated automorphism group is reported through
MapCluster.truncated, since it would over-split rather than merely slow things down.

[kfir4444]

A worked example with two real channels

CC[CH]C <=> CCC[CH2] (sec-butyl to n-butyl, intra_H_migration) is a clean case where the single-map
answer is a coin flip between two genuinely different transition states.

The reactant is C0-C1-C2(.)-C3. A hydrogen leaves a terminal methyl and lands on the radical carbon C2,
and either methyl works — both give n-butyl, since n-butyl is the same molecule numbered from either end:

C3 -> C2    1,2 shift, three-membered TS
C0 -> C2    1,3 shift, four-membered TS

These are not symmetry-equivalent. |Aut_core(R)| = 1 — sec-butyl has no automorphism exchanging its ends —
so no element of Aut(R) carries one center onto the other, and they must not merge into a single channel of
degeneracy 6. The clustering reports two channels of degeneracy 3 each, one per methyl:

maps enumerated : 8
clusters        : 2
   degeneracy 3   ('C0', 'H@C0', 1, 0), ('C2', 'H@C0', 0, 1)      1,3 shift
   degeneracy 3   ('C2', 'H@C3', 0, 1), ('C3', 'H@C3', 1, 0)      1,2 shift

All six paths are predicted by the family recipe, so the absolute reference of §"Validity filtering" did the
filtering here and neither channel is an enumeration artifact.

rxn.atom_map returns a representative of the 1,3-shift channel. The 1,2-shift channel is never reported.

Both transition states have been located and optimized independently, at b2plyp/def2tzvp, by supplying
the two atom maps to ARC by hand — declaring the product species twice and writing out each map, which is
precisely the manual workaround this PR removes the need for. The two barriers differ substantially, so the
channels are kinetically distinct rather than a bookkeeping artifact.

Feeding those two hand-written maps to cluster_atom_maps puts them in two different clusters, and the
automatic enumeration reproduces exactly the same two channel keys with no maps supplied. That equivalence is
pinned by test_sec_butyl_to_n_butyl_has_two_channels.

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 65.64%. Comparing base (f45edb0) to head (6e44935).
⚠️ Report is 20 commits behind head on main.

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #1024      +/-   ##
==========================================
+ Coverage   64.82%   65.64%   +0.81%     
==========================================
  Files         119      121       +2     
  Lines       40039    40924     +885     
  Branches    10350    10529     +179     
==========================================
+ Hits        25957    26863     +906     
+ Misses      11111    11046      -65     
- Partials     2971     3015      +44     

Flag	Coverage Δ
functionaltests	65.64% <ø> (+0.81%)	⬆️
unittests	65.64% <ø> (+0.81%)	⬆️

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]

Pull request overview

Adds an opt-in pathway to enumerate all valid atom maps for a reaction and cluster them into chemically distinct reaction channels (with per-channel reaction-path degeneracy), exposed via a new lazily computed ARCReaction.atom_map_clusters property while keeping existing atom_map behavior unchanged.

Changes:

• Introduces arc/mapping/cluster.py to enumerate atom maps across templates/directions and cluster them by reactant-symmetry–equivalent reaction centers, computing exact degeneracies.
• Extends the mapping pipeline to support multi-map enumeration (map_two_species_all, map_pairs_all, map_rxn_all) and improves reverse-discovered-template flipping via prepare_flipped_reaction.
• Adds comprehensive unit tests for clustering theory/implementation and targeted tests for reverse-discovered-template handling.

Reviewed changes

Copilot reviewed 6 out of 6 changed files in this pull request and generated 3 comments.

Show a summary per file

File	Description
arc/reaction/reaction.py	Exposes atom_map_clusters on ARCReaction and caches results.
arc/mapping/engine.py	Adds enumeration APIs (map_two_species_all, map_pairs_all) and refactors backbone scoring helpers.
arc/mapping/driver.py	Adds prepare_flipped_reaction and map_rxn_all to support flipped mapping and per-template multi-map enumeration.
arc/mapping/driver_test.py	Adds tests covering flipped reaction seeding and reverse-discovered-template mapping.
arc/mapping/cluster.py	New module implementing map enumeration, symmetry handling, clustering, and degeneracy computation.
arc/mapping/cluster_test.py	New test suite validating graph construction, automorphisms, center canonicalization, degeneracy, validity filters, and end-to-end clustering.

---

💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.

[calvinp0]

the n!/(n-k)! factor counts ordered assignments, but the centre is set - so when two hydrogens on the same parent are interchangeable it counts each pair twice.
CH2O <=> CO + H2 - 1,2 Insertion CO comes out as 2; RMG gives 1 and so does a direct orbit count under the docstring definition.

[kfir4444]

Confirmed, and fixed in 59631c4. You are exactly right about the cause: the falling factorial is the orbit of an ordered tuple of hydrogens, but the center is a set, so any permutation of hydrogens that reproduces the same set was being counted as a separate path.

Reproduced at unit level on CH2O, both C-H bonds breaking into H2:

reported=2   true orbit=1     <- both hydrogens in the same role
reported=2   true orbit=2     <- control, one hydrogen named

The fix divides the falling factorial by the stabilizer of the center under permutations of that parent's named hydrogens (hydrogen_stabilizer). For CH2O: 2!/0! = 2, one role class of size 2, 2/2! = 1. The control is untouched at 2, and a third test pins the case that must not collapse — two hydrogens on one carbon in different roles (one breaks, one changes order) stay two paths.

I cross-checked center_degeneracy against an exhaustive orbit count over the full group (core automorphisms x hydrogen bijections) rather than trusting the factorisation: CH2O 1 and 2, CH4+OH 4, C2H6+OH 6, cyclohexane+OH 12, propane 6 - all agree. End to end, CH2O <=> CO + H2 now reports degeneracy 1, matching RMG.

[calvinp0]

The degeneracy is computed from the one Lewis structure ARC holds, so a radical spread over equivalent sites looks less symmetric than it is. Measured through map_reaction_clusters, against RMG in the same forward direction:

reaction	RMG	here
allyl + CH4 → propene + CH3	8	4
NO2 + CH4 → HONO + CH3	8	4
c-C5H5 + CH4 → c-C5H6 + CH3	20	4
benzyl + CH4 → toluene + CH3	4	4

enumerated_degeneracy matches degeneracy in all four.

Benzel agrees becuase its radical carbon is already the carbon with two hydrogens and its ring sites give a different product.

[kfir4444]

Confirmed, and fixed in 97109c1. Your diagnosis was right and the mechanism is precisely the one you name: graph.invariant seeded the colour refinement with radical_electrons (and lone pairs and formal charge), all of which are properties of the one Lewis structure ARC happens to hold. Allyl's two terminal carbons were then distinguished only by which end the stored structure puts the radical on, which halved every degeneracy computed through it.

ignore_bond_orders already normalized bond orders for exactly this reason - Kekule structures - but stopped at the bonds. Resonance moves the radical, the lone pairs and the formal charge too, so effective_invariant now drops all of them when orders are ignored, leaving element and attached-hydrogen count, which resonance preserves.

Your table, re-measured:

reaction	RMG	before	now
allyl + CH4	8	4	8
NO2 + CH4	8	4	8
c-C5H5 + CH4	20	4	20
benzyl + CH4	4	4	4

Benzyl is the important control, and your explanation of why it already agreed is the reason it is in the test: its radical carbon is genuinely unique, so it must stay 4 and not inflate. It does.

The trade-off is documented on effective_invariant: two fragments differing only in charge or radical count, with identical element, connectivity and hydrogen count, are now seen as interchangeable. That is the correct reading for resonance forms of one species, and the price paid elsewhere.

[calvinp0]

for C2H5OO <=> C2H4 + HO2 and C2H5 + C2H5 <=> C2H6 + C2H4, the map that rxn.atom_map returns is not in any cluster. Its center has 11 and 10 changed bonds against the recipes 3 and 2. Would be good to log it?

[kfir4444]

Agreed, and it is now logged - warn_if_atom_map_is_unclustered, in 97109c1. It reports the map's center size against the channels', at warning level:

The atom map of <rxn> is in none of its N enumerated reaction channels: it changes 11 bonds,
against [3] for the channels. atom_map is most likely scrambled for this reaction.

Two deliberate constraints. It only inspects the atom map when it has already been computed - reading rxn.atom_map there would run the whole single-map pipeline as a side effect of clustering, which this module goes out of its way never to do, and there is a test pinning that it does not.

Second, and worth stating plainly: I think this is evidence about map_reaction, not about clustering. The channels are filtered against the family recipe, so a map outside all of them with 11 changed bonds against a recipe predicting 3 means the single-map pipeline returned a scrambled result. I hit the same thing independently on CH3 + CH3 <=> C2H6, where the sole enumerated map sends one methyl's hydrogens to the other methyl's carbon; handing the chemically correct map to the same code gives a clean one-bond center. So there is a real mapping defect upstream. I have not tried to fix map_rxn_all here - it is a much bigger change than this PR - but it is now visible rather than silent.

[calvinp0]

expected_reaction_centers returns None for every R_Recombination reaction, so only the remaining fileter is filter_minimal_centers.
CH3 + CH3 <=> C2H6 reports one channel with a 13-bond centre and degeneracy 36.

[kfir4444]

Confirmed and fixed in 97109c1, and the cause was recoverable. R_Recombination's recipe is FORM_BOND(*, *) - one label, two atoms - and get_expected_changing_bonds reads only the base label, so both endpoints resolved to the same index and the bond looked degenerate. But the label map does carry the second site under a suffixed key:

r_label_map: {'*': 0, '*_2': 4}      get_expected_changing_bonds -> ([], [(0, 0)])

resolve_duplicate_label recovers it, so the predicted center for CH3 + CH3 <=> C2H6 is now {(0, 4, 0, 1)} - the forming C-C bond, correct. Only when no suffixed sibling exists is the dictionary skipped, so the existing guard against a genuinely degenerate action still holds.

The degeneracy is fixed as a result: 36 -> 1.

The 13-bond center, though, is not a clustering artifact - it is a real defect in the enumeration, and the now-working recipe reference proves it by rejecting every enumerated map. The single map map_rxn_all produces sends one methyl's hydrogens to the other methyl's carbon; the correct map [0, 2, 3, 4, 1, 5, 6, 7] gives a one-bond center and degeneracy 1 through the same code path. Same root cause as your comment on rxn.atom_map.

Since the fallback to filter_minimal_centers cannot reject a uniformly scrambled enumeration, MapCluster now carries recipe_validated. A channel the recipe rejected is no longer presented as recipe-backed, so False next to an implausible center reads as "the enumeration is wrong here", not "this is a reaction channel".

[calvinp0]

1,2-Birad_to_alkene and Singlet_Val6_to_triplet only change bond orders, so with ignore_bond_orders=True the changed-bond set is empty and the cluster reports degeneracy=0. Was check on [CH2][CH2] <=> C=C. The empty centre also switches off both validity filters,

[kfir4444]

Confirmed on both halves, and fixed in 59631c4 and 97109c1.

The degeneracy first: an empty center returned 0, which is not a meaningful reaction path degeneracy for a reaction that exists. The orbit of the empty set under any group is the empty set, so the answer is 1, and it now returns 1.

That alone would have been papering over the real problem, which is the one you point at: with ignore_bond_orders=True the center of an order-only family is empty, and an empty center carries no information at all - nothing to cluster on, and, as you note, nothing for either validity filter to judge. So cluster_atom_maps now detects that every center came out empty and re-clusters honoring bond orders, which is the only setting that can describe these families.

On [CH2][CH2] <=> C=C:

center, orders ignored : []
center, orders honored : [(0, 1, 1.0, 2.0)]
result                 : 1 channel, degeneracy 1

The C-C bond going single to double is the reaction center, and one channel of degeneracy 1 is right - there is only one way to pair the two radical electrons. Note the recipe-based reference is unavailable with orders honored, since the recipe describes only breaking and formation, so validation falls back to the relative minimal-center filter there; that is called out in the code comment.

[calvinp0]

i think the clamp is unreachable?

[kfir4444]

You are right, it is unreachable, and it is gone in 59631c4.

named[parent] is built only from hydrogens whose parent is that atom, so it is a subset of the hydrogens on it: len(hydrogens) <= available. The loop runs offset over range(len(hydrogens)), so the smallest factor is available - (len(hydrogens) - 1) >= 1, and max(..., 1) can never bind. Measured across the graphs in the test suite the minimum factor is exactly 1, i.e. the clamp sits precisely at its boundary and never below it.

Replaced with a comment stating the subset relation, since that is the fact that makes the clamp unnecessary and is worth keeping in view.

[calvinp0]

Not sure if this is impactful on the work but RMG values for degeneracy differ by a factor of 2 whenever the reactants are identical - C2H5 + C2H5 <=> C4H10 is 1 here and 0.5 in RMG (RMG halves deliberately in reduce_same_reactant_degeneracy, Bishop & Laidler 1965)

[kfir4444]

Good catch, and it is a genuine convention difference rather than an error on either side, so I have documented it rather than changed the number.

What this reports is the orbit size of the reaction center under the full automorphism group of the reactant complex. RMG deliberately halves that for identical reactants in reduce_same_reactant_degeneracy (Bishop and Laidler 1965), so that a rate expression does not double count two indistinguishable colliding partners. Both are correct for their own purpose; they answer different questions.

Since ARC exposes this as a property of the channel rather than as a rate coefficient, halving it here would be wrong for anyone counting distinct transition states, and not halving it is wrong for anyone folding it straight into a rate. The center_degeneracy docstring now states which of the two it is, names reduce_same_reactant_degeneracy and the Bishop and Laidler reference, gives your C2H5 + C2H5 <=> C4H10 example as 1 here against 0.5 in RMG, and says explicitly that a caller building a rate coefficient for identical reactants has to apply the factor itself.

If ARC ever consumes this in a rate expression, that is the place it has to be applied, and I would rather it be an explicit step there than a silent factor buried in the orbit count.