Skip to content

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

Open
kfir4444 wants to merge 8 commits into
mainfrom
atom_map_clustering
Open

Return every atom map per reaction and cluster them into distinct reaction channels#1024
kfir4444 wants to merge 8 commits into
mainfrom
atom_map_clustering

Conversation

@kfir4444

Copy link
Copy Markdown
Collaborator

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

Copy link
Copy Markdown
Collaborator Author

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 Bot commented Aug 24, 2026

Copy link
Copy Markdown

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 64.67%. Comparing base (c4c2db9) to head (3f2346f).
⚠️ Report is 4 commits behind head on main.

Additional details and impacted files
@@            Coverage Diff             @@
##             main    #1024      +/-   ##
==========================================
+ Coverage   64.61%   64.67%   +0.05%     
==========================================
  Files         119      120       +1     
  Lines       39788    40294     +506     
  Branches    10307    10411     +104     
==========================================
+ Hits        25708    26059     +351     
- Misses      11101    11239     +138     
- Partials     2979     2996      +17     
Flag Coverage Δ
functionaltests 64.67% <ø> (+0.05%) ⬆️
unittests 64.67% <ø> (+0.05%) ⬆️

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 AI left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Comment thread arc/mapping/engine.py
Comment on lines +118 to +120
if inc_vals is not None:
atom_map = [value + inc_vals for value in atom_map]
return atom_map
Comment thread arc/reaction/reaction.py
Comment on lines +192 to +197
if self._atom_map_clusters is None \
and all(species.get_xyz(generate=False) is not None for species in self.r_species + self.p_species):
self._atom_map_clusters = map_reaction_clusters(rxn=self, backend='ARC')
if not self._atom_map_clusters:
logger.error(f"The requested ARC reaction {self} could not be atom mapped into channels.")
return self._atom_map_clusters
Comment thread arc/mapping/cluster.py
Comment on lines +658 to +664
for atom_map in atom_maps:
if len(atom_map) != r_graph.n_atoms:
logger.warning(f'Skipping an atom map of length {len(atom_map)} for {rxn}, '
f'expected {r_graph.n_atoms}.')
continue
scored.append((atom_map, changed_bonds(r_graph, p_graph, atom_map,
ignore_bond_orders=ignore_bond_orders)))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants