Return every atom map per reaction and cluster them into distinct reaction channels - #1024
Return every atom map per reaction and cluster them into distinct reaction channels#1024kfir4444 wants to merge 8 commits into
Conversation
A worked example with two real channels
The reactant is These are not symmetry-equivalent. All six paths are predicted by the family recipe, so the absolute reference of §"Validity filtering" did the
Both transition states have been located and optimized independently, at b2plyp/def2tzvp, by supplying Feeding those two hand-written maps to |
Codecov Report✅ All modified and coverable lines are covered by tests. 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
Flags with carried forward coverage won't be shown. Click here to find out more. ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
444f651 to
3f2346f
Compare
There was a problem hiding this comment.
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.pyto 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 viaprepare_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.
| if inc_vals is not None: | ||
| atom_map = [value + inc_vals for value in atom_map] | ||
| return atom_map |
| 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 |
| 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))) |
Stacked on #1023 — please review that first; this branch's base will collapse once it lands.
Problem
map_reactionreturns one atom map per reaction. Several places in the pipeline compute more than oneand discard the rest:
map_two_speciesscores every superimposable backbone candidate, then keeps the lowest-RMSD one andbreaks ties by taking the last tied candidate — an arbitrary choice.
map_rxnreturns on the first templateproduct_dictthat succeeds, and only advances to the next onewhen 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. Eachclass is one chemically distinct reaction channel (one TS to search for); its
degeneracyis that channel'sreaction path degeneracy.
Exposed as
ARCReaction.atom_map_clusters— lazily computed, cached, never triggered byatom_map, and notpersisted by
as_dict.atom_mapitself is untouched, so every existing consumer is unaffected.The equivalence relation
Two maps are equivalent iff
sigma_2 = beta . sigma_1 . alphafor somealphainAut(R),betainAut(P)— the double cosets. That reduces to a single group. WritingC(sigma)for the set of bonds ofRthat break, form, or change order, and
R'forRwith those changes applied, any validsigmais anisomorphism
R' -> P. So equal changed-bond sets differ by an element ofAut(P), and:Aut(P)never appears. Automorphisms are enumerated in-module by colour-pruned backtracking over thecore 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:
Verified against known values:
Cyclohexane exercises both halves: 6 equivalent carbons times 2 hydrogens each.
Validity filtering
changed_bondscomputes a center for any element-preserving bijection, so scrambled maps would otherwiseopen 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_bondsover eachproduct_dict'sr_label_map), falling back toa 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 anend-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.
Corpus effect
testing/errs_no_li.yml, 452 reactions, 14 workers, 200 s each: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_reactionresolves 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_clustersis opt-in andnever triggered by
atom_map. Caps exist at three levels (MAX_AUTOMORPHISMS,MAX_ENUMERATED_MAPS,map_rxn_all'smax_maps) and all log when hit; a truncated automorphism group is reported throughMapCluster.truncated, since it would over-split rather than merely slow things down.