From a24b3f4439732ff071be2c9b247687ffdd2c8294 Mon Sep 17 00:00:00 2001 From: kfir4444 Date: Sun, 23 Aug 2026 18:51:35 +0300 Subject: [PATCH 1/8] Seed the flipped reaction with a family so the reverse-discovery retry can map it --- arc/mapping/driver.py | 38 +++++++++++++++++++++++++++++++++++- arc/mapping/driver_test.py | 40 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 77 insertions(+), 1 deletion(-) diff --git a/arc/mapping/driver.py b/arc/mapping/driver.py index 395fdff36a..70fd719c5d 100644 --- a/arc/mapping/driver.py +++ b/arc/mapping/driver.py @@ -61,7 +61,7 @@ def try_mapping(r: ARCReaction) -> list[int] | None: except ValueError: return None if flip: - raw_map = try_mapping(rxn.flip_reaction()) + raw_map = try_mapping(prepare_flipped_reaction(rxn)) if raw_map is None: return None return check_atom_map_and_return(flip_map(raw_map)) @@ -80,6 +80,42 @@ def try_mapping(r: ARCReaction) -> list[int] | None: return check_atom_map_and_return(raw_map) +def prepare_flipped_reaction(rxn: ARCReaction) -> ARCReaction: + """ + Build the flipped reaction and give it a usable, forward-discovered family template. + + ``ARCReaction.flip_reaction`` resets the family, so the flipped copy re-derives its product dictionaries + lazily with the *default* family set. When the original reaction was matched only by a broader set - the + usual situation for a template discovered in reverse - that rederivation comes back empty and the flip + retry has nothing to map with. Widen the search in that case, and prefer a template that describes the + flipped reaction forward. + + Args: + rxn (ARCReaction): The reaction to flip. + + Returns: + ARCReaction: The flipped reaction, seeded with a family and product dictionaries where possible. + """ + flipped = rxn.flip_reaction() + try: + product_dicts = flipped.product_dicts or list() + except (ValueError, KeyError, AttributeError): + product_dicts = list() + if not any(not pd.get('discovered_in_reverse') for pd in product_dicts): + try: + widened = flipped.get_product_dicts(rmg_family_set='all') + except (ValueError, KeyError, AttributeError): + widened = list() + product_dicts = widened or product_dicts + forward_dicts = [pd for pd in product_dicts if not pd.get('discovered_in_reverse')] + product_dicts = forward_dicts or product_dicts + if product_dicts: + flipped.product_dicts = product_dicts + flipped.family = product_dicts[0]['family'] + flipped.family_own_reverse = product_dicts[0]['own_reverse'] + return flipped + + def check_atom_map_and_return(atom_map: list[int] | None) -> list[int] | None: """ Check if the atom map is valid and return it. diff --git a/arc/mapping/driver_test.py b/arc/mapping/driver_test.py index cd13bf6356..f131559ac4 100644 --- a/arc/mapping/driver_test.py +++ b/arc/mapping/driver_test.py @@ -860,6 +860,46 @@ def test_map_ho2_elimination_from_peroxy_radical(self): self.assertIn(tuple(rxn_4.atom_map[5: 8]), list(permutations([3, 4, 5]))) self.assertEqual(rxn_4.atom_map[8:], [6, 9]) + def test_prepare_flipped_reaction_seeds_a_forward_template(self): + """Test that the flipped reaction is given a family that matches it in the forward direction. + + ``flip_reaction`` resets the family, so the flipped copy would otherwise re-derive its product + dictionaries with the default family set, come back empty, and leave the flip retry with nothing. + """ + rxn = ARCReaction(r_species=[ARCSpecies(label='C2H5Cl', smiles='CCCl')], + p_species=[ARCSpecies(label='C2H4', smiles='C=C'), + ARCSpecies(label='HCl', smiles='Cl')]) + rxn.product_dicts = rxn.get_product_dicts(rmg_family_set='all') + flipped = prepare_flipped_reaction(rxn) + self.assertIsNotNone(flipped.family) + self.assertTrue(flipped.product_dicts) + self.assertFalse(flipped.product_dicts[0]['discovered_in_reverse']) + # The template products of a forward discovery are isomorphic to that reaction's own products. + template_products = flipped.product_dicts[0]['products'] + self.assertEqual(len(template_products), len(flipped.p_species)) + self.assertTrue(any(flipped.p_species[0].is_isomorphic(mol) for mol in template_products)) + + def test_map_reaction_with_a_reverse_discovered_template(self): + """Test mapping an elimination whose family only matches it in the addition direction. + + ``XY_Addition_MultipleBond`` matches HX elimination only in reverse, so the template's 'products' + are isomorphic to the reaction's reactants. Both ``get_template_product_order`` and + ``reorder_p_label_map`` compare against the reaction's products, so the forward attempt fails until + the reaction is flipped first. + """ + rxn = ARCReaction(r_species=[ARCSpecies(label='C2H5Cl', smiles='CCCl')], + p_species=[ARCSpecies(label='C2H4', smiles='C=C'), + ARCSpecies(label='HCl', smiles='Cl')]) + rxn.product_dicts = rxn.get_product_dicts(rmg_family_set='all') + rxn.family = rxn.product_dicts[0]['family'] + rxn.family_own_reverse = rxn.product_dicts[0]['own_reverse'] + self.assertEqual(rxn.family, 'XY_Addition_MultipleBond') + self.assertTrue(rxn.product_dicts[0]['discovered_in_reverse']) + atom_map = map_reaction(rxn=rxn, backend='ARC') + self.assertIsNotNone(atom_map) + self.assertEqual(sorted(atom_map), list(range(len(atom_map)))) + self.assertTrue(check_atom_map(rxn)) + def test_map_flipped_reaction(self): """Test the map_flipped_reaction() function.""" c2h5o3_xyz = {'coords': ((-1.3476727508427788, -0.49923624257482285, -0.3366372557370102), From 7f8bcc21fae8d8adadd274cdf02aeb6bcf7887f7 Mon Sep 17 00:00:00 2001 From: kfir4444 Date: Mon, 24 Aug 2026 14:38:38 +0300 Subject: [PATCH 2/8] Restrict the flipped reaction's product dictionaries to the family that was chosen --- arc/mapping/driver.py | 14 ++++++++++++-- arc/mapping/driver_test.py | 5 ++++- 2 files changed, 16 insertions(+), 3 deletions(-) diff --git a/arc/mapping/driver.py b/arc/mapping/driver.py index 70fd719c5d..63a5dbd034 100644 --- a/arc/mapping/driver.py +++ b/arc/mapping/driver.py @@ -101,17 +101,27 @@ def prepare_flipped_reaction(rxn: ARCReaction) -> ARCReaction: product_dicts = flipped.product_dicts or list() except (ValueError, KeyError, AttributeError): product_dicts = list() - if not any(not pd.get('discovered_in_reverse') for pd in product_dicts): + if all(pd.get('discovered_in_reverse') for pd in product_dicts): try: widened = flipped.get_product_dicts(rmg_family_set='all') except (ValueError, KeyError, AttributeError): widened = list() + logger.debug(f'Widened the family search for {flipped.label} to all families, ' + f'got {len(widened)} product dictionaries.') product_dicts = widened or product_dicts forward_dicts = [pd for pd in product_dicts if not pd.get('discovered_in_reverse')] product_dicts = forward_dicts or product_dicts if product_dicts: + # Keep only the chosen family's dictionaries. get_reaction_family_products() concatenates matches + # across every family without grouping them, while ``family`` can only name one, so leaving the + # whole list in place lets map_rxn's product_dict_index retry pair one family's recipe with + # another's label map. Every family labels its atoms *1/*2/*3, so that mispairing resolves against + # the wrong atoms and still yields a permutation, which is all check_atom_map_and_return verifies. + # TODO: replace with ARCReaction.restrict_product_dicts_to_family() once #978 lands. + family = product_dicts[0]['family'] + product_dicts = [pd for pd in product_dicts if pd['family'] == family] flipped.product_dicts = product_dicts - flipped.family = product_dicts[0]['family'] + flipped.family = family flipped.family_own_reverse = product_dicts[0]['own_reverse'] return flipped diff --git a/arc/mapping/driver_test.py b/arc/mapping/driver_test.py index f131559ac4..4e9e65b33a 100644 --- a/arc/mapping/driver_test.py +++ b/arc/mapping/driver_test.py @@ -871,9 +871,12 @@ def test_prepare_flipped_reaction_seeds_a_forward_template(self): ARCSpecies(label='HCl', smiles='Cl')]) rxn.product_dicts = rxn.get_product_dicts(rmg_family_set='all') flipped = prepare_flipped_reaction(rxn) - self.assertIsNotNone(flipped.family) + self.assertEqual(flipped.family, 'XY_Addition_MultipleBond') self.assertTrue(flipped.product_dicts) self.assertFalse(flipped.product_dicts[0]['discovered_in_reverse']) + # Only the chosen family's dictionaries survive, so a recipe can never be paired with another + # family's label map. + self.assertEqual({pd['family'] for pd in flipped.product_dicts}, {'XY_Addition_MultipleBond'}) # The template products of a forward discovery are isomorphic to that reaction's own products. template_products = flipped.product_dicts[0]['products'] self.assertEqual(len(template_products), len(flipped.p_species)) From cc43e8efd97dff30750c93bf3d6814d616b9968b Mon Sep 17 00:00:00 2001 From: kfir4444 Date: Sun, 23 Aug 2026 18:33:45 +0300 Subject: [PATCH 3/8] Return every superimposable atom map instead of only the best-scoring one --- arc/mapping/driver.py | 100 ++++++++++++++ arc/mapping/engine.py | 307 ++++++++++++++++++++++++++++++++++-------- 2 files changed, 352 insertions(+), 55 deletions(-) diff --git a/arc/mapping/driver.py b/arc/mapping/driver.py index 63a5dbd034..f9ed4de3b9 100644 --- a/arc/mapping/driver.py +++ b/arc/mapping/driver.py @@ -7,6 +7,7 @@ 3) If the reaction is supported by RMG, it is sent to the driver. Else, it is mapped with map_general_rxn. """ +from itertools import product from typing import TYPE_CHECKING from arc.common import logger @@ -23,6 +24,7 @@ label_species_atoms, make_bond_changes, map_pairs, + map_pairs_all, pairing_reactants_and_products_for_mapping, reorder_p_label_map, update_xyz, @@ -393,6 +395,104 @@ def map_rxn(rxn: ARCReaction, return atom_map +def map_rxn_all(rxn: ARCReaction, + backend: str = 'ARC', + product_dict_index_to_try: int = 0, + max_maps: int = 200, + max_per_pair: int | None = None, + ) -> list[list[int]]: + """ + The multiplicity-preserving counterpart of :func:`map_rxn`: run the family-guided pipeline for one + template product dictionary and return *every* atom map it can produce, rather than only one. + + The pipeline stages are identical to :func:`map_rxn` up to the point where fragment pairs are mapped. + There, :func:`map_pairs_all` returns all superimposable maps per fragment pair instead of the single + best-scoring one, and every combination across fragments is glued into a full-reaction map. + + Unlike :func:`map_rxn` this function does not recurse to the next product dictionary on failure - it + returns an empty list instead. Sweeping product dictionaries is the caller's job, see + ``arc.mapping.cluster.enumerate_atom_maps``. + + Args: + rxn (ARCReaction): An ARCReaction object instance belonging to an RMG reaction family. + backend (str, optional): Currently only supports ``'ARC'``. + product_dict_index_to_try (int, optional): The index of the reaction family product dictionary to use. + max_maps (int, optional): Stop after gluing this many combinations. + max_per_pair (int, optional): Keep at most this many candidate maps per fragment pair. + + Returns: + list[list[int]]: The distinct atom maps produced for this product dictionary, possibly empty. + """ + pdi = product_dict_index_to_try + reactants, products = rxn.get_reactants_and_products(return_copies=False) + reactants, products = copy_species_list_for_mapping(reactants), copy_species_list_for_mapping(products) + label_species_atoms(reactants), label_species_atoms(products) + + r_bdes = find_all_breaking_bonds(rxn, r_direction=True, pdi=pdi) + p_bdes = find_all_breaking_bonds(rxn, r_direction=False, pdi=pdi) + r_cuts = cut_species_based_on_atom_indices(reactants, r_bdes) + p_cuts = cut_species_based_on_atom_indices(products, p_bdes) + if r_cuts is None or p_cuts is None: + logger.debug(f'map_rxn_all (rxn={rxn}, pdi={pdi}): could not cut species.') + return list() + + try: + r_label_map = rxn.product_dicts[pdi]['r_label_map'] + p_label_map = rxn.product_dicts[pdi]['p_label_map'] + template_products = rxn.product_dicts[pdi]['products'] + except (IndexError, KeyError) as e: + logger.debug(f'map_rxn_all (rxn={rxn}, pdi={pdi}): no valid template maps. Got:\n{e}') + return list() + try: + template_order = get_template_product_order(rxn, template_products) + except ValueError: + logger.debug(f'map_rxn_all (rxn={rxn}, pdi={pdi}): no valid template order.') + return list() + + updated_p_label_map = reorder_p_label_map(p_label_map=p_label_map, + template_order=template_order, + template_products=template_products, + actual_products=rxn.get_reactants_and_products()[1]) + try: + make_bond_changes(rxn, r_cuts, r_label_map) + except (ValueError, IndexError, ActionError, AtomTypeError) as e: + logger.warning(e) + r_cuts, p_cuts = update_xyz(r_cuts), update_xyz(p_cuts) + pairs = pairing_reactants_and_products_for_mapping(r_cuts, p_cuts) + if p_cuts: + logger.debug(f'map_rxn_all (rxn={rxn}, pdi={pdi}): unpaired scissored products remain.') + return list() + + fragment_map_options = map_pairs_all(pairs, max_per_pair=max_per_pair) + if fragment_map_options is None: + logger.debug(f'map_rxn_all (rxn={rxn}, pdi={pdi}): one or more fragment pairs could not be mapped.') + return list() + + total_atoms = sum(len(sp.mol.atoms) for sp in reactants) + atom_maps, seen = list(), set() + for combination in product(*fragment_map_options): + if len(atom_maps) >= max_maps: + logger.debug(f'map_rxn_all (rxn={rxn}, pdi={pdi}): reached the cap of {max_maps} maps.') + break + try: + atom_map = glue_maps(maps=list(combination), + pairs=pairs, + r_label_map=r_label_map, + p_label_map=updated_p_label_map, + total_atoms=total_atoms, + ) + except (ValueError, IndexError, KeyError) as e: + logger.debug(f'map_rxn_all (rxn={rxn}, pdi={pdi}): gluing a combination failed with {e!r}.') + continue + if atom_map is None: + continue + key = tuple(atom_map) + if key not in seen: + seen.add(key) + atom_maps.append(atom_map) + return atom_maps + + def convert_label_dict(label_dict: dict[str, int], reference_mol_list: list[Molecule], mol_list: list[Molecule], diff --git a/arc/mapping/engine.py b/arc/mapping/engine.py index b67c90a0f4..650654cac7 100644 --- a/arc/mapping/engine.py +++ b/arc/mapping/engine.py @@ -79,6 +79,151 @@ def map_two_species(spc_1: ARCSpecies | Molecule, logger.warning(f'Could not map species {spc_1} and {spc_2}.') return None + trivial_map = trivial_atom_map(spc_1, spc_2, map_type=map_type) + if trivial_map is not None: + return trivial_map + + if backend.lower() not in ['arc']: + raise ValueError(f'The backend {backend} is not supported for 3DAM.') + atom_map = None + + if backend.lower() == 'arc': + candidates = identify_backbone_candidates(spc_1, spc_2, consider_chirality=consider_chirality) + if not candidates: + return None + rmsds, fixed_spcs = list(), list() + for candidate in candidates: + rmsd, fixed_spc_1, fixed_spc_2 = score_backbone_candidate(spc_1, spc_2, candidate) + rmsds.append(rmsd) + fixed_spcs.append((fixed_spc_1, fixed_spc_2)) + lowest_rmsd = min(rmsds) + tied_indices = [i for i, rmsd in enumerate(rmsds) if rmsd - lowest_rmsd <= RMSD_TIE_TOLERANCE] + atom_maps = dict() + for i in tied_indices: + atom_maps[i] = map_hydrogens(fixed_spcs[i][0], fixed_spcs[i][1], candidates[i]) + check_atom_map_is_a_permutation(atom_maps[i], spc_1, spc_2) + if len(tied_indices) > 1: + displacements = {i: fixed_spcs[i][0].kabsch(fixed_spcs[i][1], + [v for k, v in sorted(atom_maps[i].items(), + key=lambda item: item[0])]) + for i in tied_indices} + lowest_displacement = min(displacements.values()) + tied_indices = [i for i in tied_indices + if displacements[i] - lowest_displacement <= RMSD_TIE_TOLERANCE] + chosen_candidate_index = max(tied_indices) + atom_map = atom_maps[chosen_candidate_index] + if map_type == 'list': + atom_map = [v for k, v in sorted(atom_map.items(), key=lambda item: item[0])] + + if inc_vals is not None: + atom_map = [value + inc_vals for value in atom_map] + return atom_map + + +def map_two_species_all(spc_1: ARCSpecies | Molecule, + spc_2: ARCSpecies | Molecule, + map_type: str = 'list', + backend: str = 'ARC', + consider_chirality: bool = True, + inc_vals: int | None = None, + verbose: bool = False, + ) -> list[tuple[list[int] | dict[int, int], float]] | None: + """ + Map the atoms in ``spc_1`` to the atoms in ``spc_2``, returning *every* superimposable candidate map + rather than only the best-scoring one. + + This is the multiplicity-preserving counterpart of :func:`map_two_species`. Where ``map_two_species`` + collapses the candidate list produced by :func:`identify_superimposable_candidates` down to a single + map (lowest backbone RMSD, ties broken by the all-atom Kabsch displacement and then by taking the last + tied candidate), this function scores and returns all of them. It is intended for callers that need to + enumerate distinct reaction channels, see ``arc.mapping.cluster``. + + Note that the first entry of the returned list is not necessarily the map returned by + ``map_two_species``: entries are ordered by ascending backbone RMSD, whereas ``map_two_species`` + resolves a tie by choosing the *last* tied candidate. + + Candidates whose hydrogen mapping fails (a ``ValueError`` raised by :func:`map_hydrogens`, e.g. an + H-count mismatch at a mapped heavy atom) are skipped rather than aborting the whole enumeration. + + Args: + spc_1 (ARCSpecies | Molecule): Species 1. + spc_2 (ARCSpecies | Molecule): Species 2. + map_type (str, optional): Whether to return 'list' or 'dict' type maps. + backend (str, optional): Currently only ``ARC``'s method is implemented as the backend. + consider_chirality (bool, optional): Whether to consider chirality when fingerprinting. + inc_vals (int, optional): An optional integer by which all values in the atom map lists will be incremented. + verbose (bool, optional): Whether to use logging. + + Returns: + list[tuple[list[int] | dict[int, int], float]] | None: + Entries are ``(atom_map, backbone_rmsd)`` tuples ordered by ascending ``backbone_rmsd``. + ``None`` if the two species could not be mapped at all. + """ + spc_1, spc_2 = get_arc_species(spc_1), get_arc_species(spc_2) + if not check_species_before_mapping(spc_1, spc_2, verbose=verbose): + if verbose: + logger.warning(f'Could not map species {spc_1} and {spc_2}.') + return None + + trivial_map = trivial_atom_map(spc_1, spc_2, map_type=map_type) + if trivial_map is not None: + # A trivial map is unique by construction, so the enumeration holds exactly one entry. + return [(trivial_map, 0.0)] + + if backend.lower() not in ['arc']: + raise ValueError(f'The backend {backend} is not supported for 3DAM.') + + candidates = identify_backbone_candidates(spc_1, spc_2, consider_chirality=consider_chirality) + if not candidates: + return None + + scored: list[tuple[list[int] | dict[int, int], float]] = list() + for candidate in candidates: + rmsd, fixed_spc_1, fixed_spc_2 = score_backbone_candidate(spc_1, spc_2, candidate) + try: + atom_map = map_hydrogens(fixed_spc_1, fixed_spc_2, candidate) + # A candidate that does not yield a permutation is discarded rather than raised on, since + # unlike map_two_species this function is enumerating and other candidates may still be good. + check_atom_map_is_a_permutation(atom_map, spc_1, spc_2) + except ValueError as e: + if verbose: + logger.warning(f'Could not map hydrogens for a backbone candidate of ' + f'{spc_1} and {spc_2}, skipping it. Got:\n{e}') + continue + if map_type == 'list': + atom_map = [v for k, v in sorted(atom_map.items(), key=lambda item: item[0])] + if inc_vals is not None: + atom_map = [value + inc_vals for value in atom_map] + elif inc_vals is not None: + atom_map = {k: v + inc_vals for k, v in atom_map.items()} + scored.append((atom_map, rmsd)) + if not scored: + logger.warning(f'Could not map hydrogens for any backbone candidate of {spc_1} and {spc_2}.') + return None + scored.sort(key=lambda entry: entry[1]) + return scored + + +def trivial_atom_map(spc_1: ARCSpecies, + spc_2: ARCSpecies, + map_type: str = 'list', + ) -> list[int] | dict[int, int] | None: + """ + Return the atom map for the trivial cases that do not require the fingerprint/DFS/RMSD pipeline: + mono-atomic species, homonuclear diatomic species, and species in which every atom is a different element. + + Note: Historically these shortcuts returned directly from ``map_two_species`` before the ``inc_vals`` + increment was applied, so a trivial map is *not* incremented by ``inc_vals``. That behaviour is preserved + here. ``inc_vals`` is only exercised by the test suite, it is not used anywhere in the production pipeline. + + Args: + spc_1 (ARCSpecies): Species 1. + spc_2 (ARCSpecies): Species 2. + map_type (str, optional): Whether to return a 'list' or a 'dict' map type. + + Returns: + list[int] | dict[int, int] | None: The trivial atom map, or ``None`` if no shortcut applies. + """ # A shortcut for mono-atomic species. if spc_1.number_of_atoms == spc_2.number_of_atoms == 1: if map_type == 'dict': @@ -104,67 +249,90 @@ def map_two_species(spc_1: ARCSpecies | Molecule, atom_map = [v for k, v in sorted(atom_map.items(), key=lambda item: item[0])] return atom_map - if backend.lower() not in ['arc']: - raise ValueError(f'The backend {backend} is not supported for 3DAM.') - atom_map = None + return None - if backend.lower() == 'arc': + +def identify_backbone_candidates(spc_1: ARCSpecies, + spc_2: ARCSpecies, + consider_chirality: bool = True, + ) -> list[dict[int, int]] | None: + """ + Fingerprint both species and identify all superimposable heavy-atom backbone candidates. + If no candidate is found, the search is retried with the opposite ``consider_chirality`` setting. + + Args: + spc_1 (ARCSpecies): Species 1. + spc_2 (ARCSpecies): Species 2. + consider_chirality (bool, optional): Whether to consider chirality when fingerprinting. + + Returns: + list[dict[int, int]] | None: The backbone candidates, or ``None`` if none could be identified. + """ + fingerprint_1 = fingerprint(spc_1, consider_chirality=consider_chirality) + fingerprint_2 = fingerprint(spc_2, consider_chirality=consider_chirality) + candidates = identify_superimposable_candidates(fingerprint_1, fingerprint_2) + if candidates is None or len(candidates) == 0: + consider_chirality = not consider_chirality fingerprint_1 = fingerprint(spc_1, consider_chirality=consider_chirality) fingerprint_2 = fingerprint(spc_2, consider_chirality=consider_chirality) candidates = identify_superimposable_candidates(fingerprint_1, fingerprint_2) if candidates is None or len(candidates) == 0: - consider_chirality = not consider_chirality - fingerprint_1 = fingerprint(spc_1, consider_chirality=consider_chirality) - fingerprint_2 = fingerprint(spc_2, consider_chirality=consider_chirality) - candidates = identify_superimposable_candidates(fingerprint_1, fingerprint_2) - if candidates is None or len(candidates) == 0: - logger.warning(f'Could not identify superimposable candidates {spc_1} and {spc_2}.') - return None - if not len(candidates): - return None - else: - rmsds, fixed_spcs = list(), list() - for candidate in candidates: - fixed_spc_1, fixed_spc_2 = fix_dihedrals_by_backbone_mapping(spc_1, spc_2, backbone_map=candidate) - fixed_spcs.append((fixed_spc_1, fixed_spc_2)) - backbone_1, backbone_2 = set(list(candidate.keys())), set(list(candidate.values())) - xyz1, xyz2 = fixed_spc_1.get_xyz(), fixed_spc_2.get_xyz() - xyz1 = xyz_from_data(coords=[xyz1['coords'][i] for i in range(fixed_spc_1.number_of_atoms) if i in backbone_1], - symbols=[xyz1['symbols'][i] for i in range(fixed_spc_1.number_of_atoms) if i in backbone_1], - isotopes=[xyz1['isotopes'][i] for i in range(fixed_spc_1.number_of_atoms) if i in backbone_1]) - xyz2 = xyz_from_data(coords=[xyz2['coords'][i] for i in range(fixed_spc_2.number_of_atoms) if i in backbone_2], - symbols=[xyz2['symbols'][i] for i in range(fixed_spc_2.number_of_atoms) if i in backbone_2], - isotopes=[xyz2['isotopes'][i] for i in range(fixed_spc_2.number_of_atoms) if i in backbone_2]) - no_gap_candidate = remove_gaps_from_values(candidate) - xyz2 = sort_xyz_using_indices(xyz2, indices=[v for k, v in sorted(no_gap_candidate.items(), - key=lambda item: item[0])]) - rmsds.append(compare_confs(xyz1=xyz1, xyz2=xyz2, rmsd_score=True)) - lowest_rmsd = min(rmsds) - tied_indices = [i for i, rmsd in enumerate(rmsds) if rmsd - lowest_rmsd <= RMSD_TIE_TOLERANCE] - atom_maps = dict() - for i in tied_indices: - atom_maps[i] = map_hydrogens(fixed_spcs[i][0], fixed_spcs[i][1], candidates[i]) - if sorted(atom_maps[i].keys()) != list(range(spc_1.number_of_atoms)) \ - or sorted(atom_maps[i].values()) != list(range(spc_2.number_of_atoms)): - raise ValueError(f'The atom map of {spc_1.label} and {spc_2.label} is not a permutation of their ' - f'{spc_1.number_of_atoms} and {spc_2.number_of_atoms} atoms, ' - f'got:\n{atom_maps[i]}') - if len(tied_indices) > 1: - displacements = {i: fixed_spcs[i][0].kabsch(fixed_spcs[i][1], - [v for k, v in sorted(atom_maps[i].items(), - key=lambda item: item[0])]) - for i in tied_indices} - lowest_displacement = min(displacements.values()) - tied_indices = [i for i in tied_indices - if displacements[i] - lowest_displacement <= RMSD_TIE_TOLERANCE] - chosen_candidate_index = max(tied_indices) - atom_map = atom_maps[chosen_candidate_index] - if map_type == 'list': - atom_map = [v for k, v in sorted(atom_map.items(), key=lambda item: item[0])] + logger.warning(f'Could not identify superimposable candidates {spc_1} and {spc_2}.') + return None + return candidates - if inc_vals is not None: - atom_map = [value + inc_vals for value in atom_map] - return atom_map + +def check_atom_map_is_a_permutation(atom_map: dict[int, int], + spc_1: ARCSpecies, + spc_2: ARCSpecies, + ) -> None: + """ + Verify that a mapped species pair yields a genuine permutation of both species' atoms. + + Args: + spc_1 (ARCSpecies): Species 1. + spc_2 (ARCSpecies): Species 2. + atom_map (dict[int, int]): The candidate atom map. + + Raises: + ValueError: If the map does not cover every atom of both species exactly once. + """ + if sorted(atom_map.keys()) != list(range(spc_1.number_of_atoms)) \ + or sorted(atom_map.values()) != list(range(spc_2.number_of_atoms)): + raise ValueError(f'The atom map of {spc_1.label} and {spc_2.label} is not a permutation of their ' + f'{spc_1.number_of_atoms} and {spc_2.number_of_atoms} atoms, ' + f'got:\n{atom_map}') + + +def score_backbone_candidate(spc_1: ARCSpecies, + spc_2: ARCSpecies, + candidate: dict[int, int], + ) -> tuple[float, ARCSpecies, ARCSpecies]: + """ + Score a single backbone candidate by the RMSD between the two dihedral-corrected backbone geometries. + + Args: + spc_1 (ARCSpecies): Species 1. + spc_2 (ARCSpecies): Species 2. + candidate (dict[int, int]): The candidate backbone map. + + Returns: + tuple[float, ARCSpecies, ARCSpecies]: + The backbone RMSD, and the two dihedral-corrected species this candidate implies. + """ + fixed_spc_1, fixed_spc_2 = fix_dihedrals_by_backbone_mapping(spc_1, spc_2, backbone_map=candidate) + backbone_1, backbone_2 = set(list(candidate.keys())), set(list(candidate.values())) + xyz1, xyz2 = fixed_spc_1.get_xyz(), fixed_spc_2.get_xyz() + xyz1 = xyz_from_data(coords=[xyz1['coords'][i] for i in range(fixed_spc_1.number_of_atoms) if i in backbone_1], + symbols=[xyz1['symbols'][i] for i in range(fixed_spc_1.number_of_atoms) if i in backbone_1], + isotopes=[xyz1['isotopes'][i] for i in range(fixed_spc_1.number_of_atoms) if i in backbone_1]) + xyz2 = xyz_from_data(coords=[xyz2['coords'][i] for i in range(fixed_spc_2.number_of_atoms) if i in backbone_2], + symbols=[xyz2['symbols'][i] for i in range(fixed_spc_2.number_of_atoms) if i in backbone_2], + isotopes=[xyz2['isotopes'][i] for i in range(fixed_spc_2.number_of_atoms) if i in backbone_2]) + no_gap_candidate = remove_gaps_from_values(candidate) + xyz2 = sort_xyz_using_indices(xyz2, indices=[v for k, v in sorted(no_gap_candidate.items(), + key=lambda item: item[0])]) + return compare_confs(xyz1=xyz1, xyz2=xyz2, rmsd_score=True), fixed_spc_1, fixed_spc_2 def get_arc_species(spc: ARCSpecies | Molecule) -> ARCSpecies: @@ -1571,6 +1739,35 @@ def map_pairs(pairs: list[tuple[ARCSpecies, ARCSpecies]]) -> list[list[int]]: return maps +def map_pairs_all(pairs: list[tuple[ARCSpecies, ARCSpecies]], + max_per_pair: int | None = None, + ) -> list[list[list[int]]] | None: + """ + The multiplicity-preserving counterpart of :func:`map_pairs`: for every matched fragment pair, return + *all* superimposable maps rather than only the best-scoring one. + + Args: + pairs (list[tuple[ARCSpecies, ARCSpecies]]): The matched reactant/product fragment pairs. + max_per_pair (int, optional): Keep at most this many maps per pair, best-scoring first. + + Returns: + list[list[list[int]]] | None: + Per pair, the list of candidate atom maps ordered by ascending backbone RMSD. + ``None`` if any pair could not be mapped at all. + """ + maps = list() + for pair in pairs: + scored = map_two_species_all(pair[0], pair[1]) + if not scored: + return None + candidates = [atom_map for atom_map, _ in scored] + if max_per_pair is not None and len(candidates) > max_per_pair: + logger.debug(f'map_pairs_all: keeping {max_per_pair} of {len(candidates)} maps for a fragment pair.') + candidates = candidates[:max_per_pair] + maps.append(candidates) + return maps + + def label_species_atoms(species: list[ARCSpecies]) -> None: """ Adds the labels to the ``.mol.atoms`` properties of the species object. From 8bd96370b7be7a2c61607c4bc8e97af96d0a6e29 Mon Sep 17 00:00:00 2001 From: kfir4444 Date: Sun, 23 Aug 2026 18:33:45 +0300 Subject: [PATCH 4/8] Cluster atom maps into distinct reaction channels by their Aut(R) orbit --- arc/mapping/cluster.py | 782 ++++++++++++++++++++++++++++++++++++ arc/mapping/cluster_test.py | 562 ++++++++++++++++++++++++++ 2 files changed, 1344 insertions(+) create mode 100644 arc/mapping/cluster.py create mode 100644 arc/mapping/cluster_test.py diff --git a/arc/mapping/cluster.py b/arc/mapping/cluster.py new file mode 100644 index 0000000000..bc1f60fcf1 --- /dev/null +++ b/arc/mapping/cluster.py @@ -0,0 +1,782 @@ +""" +Enumeration and equivalence-clustering of atom maps. + +The standard mapping entry point, ``arc.mapping.driver.map_reaction``, returns a single atom map per +reaction. This module enumerates *many* valid atom maps and groups them into equivalence classes, where +each class corresponds to one chemically distinct reaction channel (one transition state to search for) +and the size of a class is that channel's reaction path degeneracy. + +Theory +------ +An atom map is a bijection ``sigma: atoms(R) -> atoms(P)`` from the reactant complex to the product +complex, stored as a list where ``atom_map[reactant_index] == product_index``. + +If ``sigma`` is valid then so is ``beta . sigma . alpha`` for any ``alpha`` in ``Aut(R)`` and ``beta`` in +``Aut(P)``, since composing with an automorphism only relabels symmetry-equivalent atoms. The equivalence +classes are therefore the double cosets ``Aut(P) . sigma . Aut(R)``. + +That definition can be reduced to something much cheaper. Let ``C(sigma)`` be the *changed-bond set*: the +bonds of ``R`` that are broken, formed, or change order under ``sigma``. Write ``R'`` for ``R`` with those +changes applied. Any valid ``sigma`` is a graph isomorphism ``R' -> P``, so: + +- If ``C(sigma_1) == C(sigma_2)`` then ``sigma_2 . sigma_1^-1`` maps ``P`` to ``P`` preserving all bonds, + i.e. it lies in ``Aut(P)``, hence ``sigma_1 ~ sigma_2``. +- If ``C(sigma_2) == alpha(C(sigma_1))`` for some ``alpha`` in ``Aut(R)``, then ``sigma_1 . alpha^-1`` has + changed-bond set ``C(sigma_2)``, so by the previous point ``sigma_2 = beta . sigma_1 . alpha^-1``. + +Therefore:: + + sigma_1 ~ sigma_2 <=> C(sigma_1) and C(sigma_2) are in the same Aut(R) orbit + +Only ``Aut(R)`` is ever needed; ``Aut(P)`` does not appear. Clustering reduces to canonicalizing an edge +subset of ``R`` under a single group action. + +Hydrogens +--------- +Automorphism groups are computed on the *core skeleton* rather than the full molecular graph, because the +full group is dominated by intra-XHn hydrogen permutations that provably cannot change a cluster (the +hydrogens on a methyl group are one orbit of ``Aut(R)`` by construction). For isobutane the full group has +order 1296 while the core skeleton group has order 6. + +A "core" atom is any non-hydrogen atom, plus any hydrogen that has no unique non-hydrogen neighbor (a free +H radical, H2, ...). Remaining hydrogens are represented in the changed-bond set by their parent core atom, +which is well defined precisely because hydrogens on a common parent are interchangeable. This keeps +hydrogen-transfer reactions - where the migrating atom *is* a hydrogen - fully representable. + +See ``docs/atom_mapping_clustering_design.md`` for the full design, and ``docs/atom_mapping_summary.md`` +for a description of the underlying single-map pipeline. +""" + +from dataclasses import dataclass, field + +from arc.common import logger +from arc.mapping.driver import MAX_PDI, map_rxn_all +from arc.mapping.engine import flip_map +from arc.species import ARCSpecies + +# An upper bound on the number of automorphisms enumerated for one complex. Hitting this bound means the +# clustering is computed under a truncated group and may split classes that are genuinely equivalent, so it +# is reported rather than applied silently. +MAX_AUTOMORPHISMS = 100_000 + +# An upper bound on the number of maps enumerated for one reaction. +MAX_ENUMERATED_MAPS = 5_000 + +# Endpoint descriptor tags used in a changed-bond set. +CORE = 'A' +HYDROGEN = 'H' + + +@dataclass +class ComplexGraph: + """ + A flat graph view of a reactant or product complex, indexed by *running* atom indices across the whole + complex - the same indexing convention used by ``atom_map``. + + Attributes: + n_atoms (int): Total number of atoms in the complex. + symbols (list[str]): Element symbol per atom. + bonds (dict[frozenset[int], float]): Bond orders keyed by the unordered pair of atom indices. + core (list[int]): Indices of core (skeleton) atoms. + parent (dict[int, int]): Maps a non-core hydrogen index to its parent core atom index. + adj (dict[int, dict[int, float]]): Core-core adjacency with bond orders. + invariant (dict[int, tuple]): Per-core-atom initial colour used to seed the automorphism search. + """ + n_atoms: int + symbols: list[str] + bonds: dict[frozenset, float] + core: list[int] + parent: dict[int, int] + adj: dict[int, dict[int, float]] + invariant: dict[int, tuple] + + +@dataclass +class MapCluster: + """ + One equivalence class of atom maps, i.e. one distinct reaction channel. + + Attributes: + representative (list[int]): A representative atom map for this channel. + members (list[list[int]]): All enumerated atom maps belonging to this channel. + centers (set): The distinct changed-bond sets seen in this channel, one per distinct reaction path. + key (tuple): The canonical ``Aut(R)``-invariant form of the changed-bond set. + signature (tuple): A human-readable orbit-level signature of the reaction center. + truncated (bool): Whether the automorphism group used to build ``key`` was truncated. + """ + representative: list[int] + members: list[list[int]] = field(default_factory=list) + centers: set = field(default_factory=set) + key: tuple = () + signature: tuple = () + truncated: bool = False + + @property + def degeneracy(self) -> int: + """ + int: The reaction path degeneracy of this channel, i.e. the number of distinct reaction centers. + + This deliberately counts distinct changed-bond sets rather than distinct atom maps. Two maps with + the same changed-bond set differ only by an element of ``Aut(P)`` - they relabel product atoms + without changing which reactant bonds break and form - so they are the same reaction path and must + not be counted twice. For CH4 + OH the four abstractions give four centers, while a map that only + swaps the two resulting water hydrogens adds a member but no new center. + + Note that this is a count of the paths actually *enumerated*. It is a lower bound on the true + degeneracy whenever the enumeration is incomplete. + """ + return len(self.centers) + + def __repr__(self) -> str: + return f'' + + +def build_complex_graph(species_list: list[ARCSpecies]) -> ComplexGraph: + """ + Build a :class:`ComplexGraph` from a list of species, concatenating their atoms in order so that atom + indices match the running-index convention of ``atom_map``. + + Args: + species_list (list[ARCSpecies]): The species forming the complex. + + Returns: + ComplexGraph: The flat graph view of the complex. + """ + symbols: list[str] = list() + bonds: dict[frozenset, float] = dict() + heavy_neighbors: dict[int, list[int]] = dict() + attributes: dict[int, tuple] = dict() + offset = 0 + for spc in species_list: + atoms = spc.mol.atoms + local_index = {id(atom): i for i, atom in enumerate(atoms)} + for i, atom in enumerate(atoms): + symbols.append(atom.element.symbol) + attributes[offset + i] = (atom.element.symbol, + getattr(atom, 'charge', 0), + getattr(atom, 'radical_electrons', 0), + getattr(atom, 'lone_pairs', 0)) + heavy_neighbors[offset + i] = list() + for bond in spc.mol.get_all_edges(): + i = offset + local_index[id(bond.atom1)] + j = offset + local_index[id(bond.atom2)] + bonds[frozenset((i, j))] = bond.order + if bond.atom2.element.symbol != 'H': + heavy_neighbors[i].append(j) + if bond.atom1.element.symbol != 'H': + heavy_neighbors[j].append(i) + offset += len(atoms) + + # A hydrogen is represented by its parent only when it has exactly one non-hydrogen neighbor. + # Everything else (a free H radical, H2, a bridging H) is promoted to a core atom in its own right. + parent: dict[int, int] = dict() + core: list[int] = list() + for i, symbol in enumerate(symbols): + if symbol == 'H' and len(heavy_neighbors[i]) == 1: + parent[i] = heavy_neighbors[i][0] + else: + core.append(i) + + core_set = set(core) + n_hydrogens = {v: 0 for v in core} + for h, p in parent.items(): + n_hydrogens[p] = n_hydrogens.get(p, 0) + 1 + + adj: dict[int, dict[int, float]] = {v: dict() for v in core} + for pair, order in bonds.items(): + i, j = tuple(pair) + if i in core_set and j in core_set: + adj[i][j] = order + adj[j][i] = order + + invariant = {v: attributes[v] + (n_hydrogens[v],) for v in core} + return ComplexGraph(n_atoms=len(symbols), symbols=symbols, bonds=bonds, + core=core, parent=parent, adj=adj, invariant=invariant) + + +def effective_adjacency(graph: ComplexGraph, ignore_bond_orders: bool = True) -> dict[int, dict[int, float]]: + """ + Return the core-core adjacency used for symmetry detection, optionally with all bond orders normalized + to 1. + + Normalizing matters because RMG stores an aromatic ring as a Kekule structure with alternating single + and double bonds. Only 6 of benzene's 12 graph automorphisms preserve that alternation, so an + order-sensitive group under-counts the symmetry and splits clusters that are genuinely equivalent. + This must agree with the ``ignore_bond_orders`` setting used by :func:`changed_bonds`. + + Note that this normalizes bond orders only. Per-atom radical and lone-pair counts are still part of + ``graph.invariant`` and remain resonance-form dependent, so a delocalized radical can still be seen as + less symmetric than it physically is. + + Args: + graph (ComplexGraph): The complex graph. + ignore_bond_orders (bool, optional): Whether to normalize all bond orders to 1. + + Returns: + dict[int, dict[int, float]]: The adjacency to use. + """ + if not ignore_bond_orders: + return graph.adj + return {v: {u: 1 for u in neighbors} for v, neighbors in graph.adj.items()} + + +def refine_colors(graph: ComplexGraph, ignore_bond_orders: bool = True) -> dict[int, int]: + """ + Compute a stable colour refinement (1-dimensional Weisfeiler-Leman) of the core skeleton, used to prune + the automorphism search. Two core atoms may only be mapped onto each other if they share a colour. + + Args: + graph (ComplexGraph): The complex graph. + ignore_bond_orders (bool, optional): Whether to ignore bond orders, see :func:`effective_adjacency`. + + Returns: + dict[int, int]: The stable colour per core atom. + """ + adj = effective_adjacency(graph, ignore_bond_orders=ignore_bond_orders) + labels = {v: graph.invariant[v] for v in graph.core} + colors = _compress(labels) + while True: + signatures = {v: (colors[v], tuple(sorted((colors[u], order) for u, order in adj[v].items()))) + for v in graph.core} + new_colors = _compress(signatures) + if len(set(new_colors.values())) == len(set(colors.values())): + return new_colors + colors = new_colors + + +def _compress(labels: dict[int, tuple]) -> dict[int, int]: + """ + Relabel arbitrary hashable colours to consecutive integers, deterministically. + + Args: + labels (dict[int, tuple]): Per-node colour. + + Returns: + dict[int, int]: Per-node integer colour. + """ + ranking = {label: i for i, label in enumerate(sorted(set(labels.values()), key=repr))} + return {v: ranking[label] for v, label in labels.items()} + + +def core_automorphisms(graph: ComplexGraph, + max_count: int = MAX_AUTOMORPHISMS, + ignore_bond_orders: bool = True, + ) -> tuple[list[dict[int, int]], bool]: + """ + Enumerate the automorphism group of the core skeleton by colour-pruned backtracking. + + Note that the complex is generally disconnected (it holds several molecules), so the group includes + permutations that exchange identical molecules. That is intended: swapping two identical reactants is a + genuine symmetry of the complex. + + Args: + graph (ComplexGraph): The complex graph. + max_count (int, optional): Stop after this many automorphisms. + ignore_bond_orders (bool, optional): Whether to ignore bond orders, see :func:`effective_adjacency`. + Must agree with the setting used by :func:`changed_bonds`. + + Returns: + tuple[list[dict[int, int]], bool]: + The automorphisms as index->index dicts over core atoms, and whether the search was truncated. + """ + adj = effective_adjacency(graph, ignore_bond_orders=ignore_bond_orders) + colors = refine_colors(graph, ignore_bond_orders=ignore_bond_orders) + buckets: dict[int, list[int]] = dict() + for v, color in colors.items(): + buckets.setdefault(color, list()).append(v) + # Place the most constrained atoms first so that contradictions surface early. + order = sorted(graph.core, key=lambda v: (len(buckets[colors[v]]), v)) + + automorphisms: list[dict[int, int]] = list() + mapping: dict[int, int] = dict() + used: set[int] = set() + truncated = False + + def backtrack(idx: int) -> None: + nonlocal truncated + if truncated: + return + if idx == len(order): + automorphisms.append(dict(mapping)) + if len(automorphisms) >= max_count: + truncated = True + return + v = order[idx] + for w in buckets[colors[v]]: + if w in used: + continue + # Verify both the presence and the absence of every edge to an already-placed atom. + if any(adj[v].get(u) != adj[w].get(image) for u, image in mapping.items()): + continue + mapping[v] = w + used.add(w) + backtrack(idx + 1) + del mapping[v] + used.discard(w) + if truncated: + return + + backtrack(0) + if truncated: + logger.warning(f'Automorphism enumeration hit the cap of {max_count}; atom map clustering may ' + f'report more channels than actually exist.') + return automorphisms, truncated + + +def core_orbits(graph: ComplexGraph, automorphisms: list[dict[int, int]]) -> dict[int, int]: + """ + Compute the orbits of the core atoms under the automorphism group, as a canonical orbit id per atom. + + Args: + graph (ComplexGraph): The complex graph. + automorphisms (list[dict[int, int]]): The automorphism group. + + Returns: + dict[int, int]: Orbit id per core atom, the id being the smallest atom index in that orbit. + """ + representative = {v: v for v in graph.core} + + def find(v: int) -> int: + while representative[v] != v: + representative[v] = representative[representative[v]] + v = representative[v] + return v + + for alpha in automorphisms: + for v, w in alpha.items(): + root_v, root_w = find(v), find(w) + if root_v != root_w: + if root_w < root_v: + root_v, root_w = root_w, root_v + representative[root_w] = root_v + return {v: find(v) for v in graph.core} + + +def changed_bonds(r_graph: ComplexGraph, + p_graph: ComplexGraph, + atom_map: list[int], + ignore_bond_orders: bool = True, + ) -> frozenset: + """ + Compute the changed-bond set ``C(sigma)`` of an atom map, expressed in *reactant-complex* indices. + + Each entry is ``(i, j, order_before, order_after)`` with ``i < j`` reactant atom indices. A broken bond + has ``order_after == 0``, a formed bond has ``order_before == 0``. + + Because the set is expressed entirely in reactant indices, it is invariant under any relabeling of the + product atoms, i.e. under ``Aut(P)``. Two maps sharing a changed-bond set therefore describe the same + reaction path, which is what makes this the right thing to count for reaction path degeneracy. For the + ``Aut(R)``-canonical cluster key the set must first be passed through :func:`collapse_hydrogens`. + + Args: + r_graph (ComplexGraph): The reactant complex. + p_graph (ComplexGraph): The product complex. + atom_map (list[int]): The atom map, ``atom_map[reactant_index] == product_index``. + ignore_bond_orders (bool, optional): If ``True``, only bond breaking and formation are recorded and + pure bond-order changes are ignored. This is the default + because the mapping pipeline mutates bond orders when handling + resonance (see ``make_bond_changes``), which would otherwise + inject spurious entries. + + Returns: + frozenset: The changed-bond set. + """ + changes = set() + inverse_map = {product_index: reactant_index for reactant_index, product_index in enumerate(atom_map)} + pairs = set(r_graph.bonds.keys()) + for pair in p_graph.bonds.keys(): + i, j = tuple(pair) + pairs.add(frozenset((inverse_map[i], inverse_map[j]))) + for pair in pairs: + i, j = tuple(pair) + order_before = r_graph.bonds.get(frozenset((i, j)), 0) + order_after = p_graph.bonds.get(frozenset((atom_map[i], atom_map[j])), 0) + if ignore_bond_orders: + order_before = 1 if order_before else 0 + order_after = 1 if order_after else 0 + if order_before == order_after: + continue + changes.add((min(i, j), max(i, j), order_before, order_after)) + return frozenset(changes) + + +def collapse_hydrogens(center: frozenset, graph: ComplexGraph) -> frozenset: + """ + Rewrite a changed-bond set in terms of core atoms, replacing each non-core hydrogen by its parent core + atom. This is what makes the set act-able by the core-skeleton automorphism group, which is only defined + on core atoms. + + Each entry becomes ``(endpoint_a, endpoint_b, order_before, order_after)`` with the two endpoint + descriptors sorted, where a descriptor is ``('A', core_index)`` for a core atom and + ``('H', parent_core_index)`` for a collapsed hydrogen. + + The collapse is lossy by design: abstracting any of methane's four hydrogens yields the same collapsed + set, which is precisely why all four land in one cluster. Reaction path degeneracy must therefore be + counted on the uncollapsed sets, see :attr:`MapCluster.degeneracy`. + + Args: + center (frozenset): The changed-bond set in reactant indices. + graph (ComplexGraph): The reactant complex graph. + + Returns: + frozenset: The collapsed changed-bond set. + """ + collapsed = set() + for i, j, order_before, order_after in center: + endpoints = tuple(sorted((_endpoint(graph, i), _endpoint(graph, j)))) + collapsed.add((endpoints[0], endpoints[1], order_before, order_after)) + return frozenset(collapsed) + + +def _endpoint(graph: ComplexGraph, index: int) -> tuple: + """ + Describe an atom as a changed-bond endpoint, collapsing a hydrogen onto its parent core atom. + + Args: + graph (ComplexGraph): The complex graph. + index (int): The atom index. + + Returns: + tuple: The endpoint descriptor. + """ + if index in graph.parent: + return HYDROGEN, graph.parent[index] + return CORE, index + + +def canonical_center_key(center: frozenset, automorphisms: list[dict[int, int]]) -> tuple: + """ + Canonicalize a changed-bond set under the reactant automorphism group: the key is the lexicographic + minimum over the group orbit, so two maps share a key exactly when they are equivalent. + + Args: + center (frozenset): The changed-bond set. + automorphisms (list[dict[int, int]]): The reactant automorphism group. + + Returns: + tuple: The canonical, hashable form of the changed-bond set. + """ + if not automorphisms: + return tuple(sorted(center)) + best = None + for alpha in automorphisms: + image = list() + for endpoint_a, endpoint_b, order_before, order_after in center: + mapped_a = (endpoint_a[0], alpha.get(endpoint_a[1], endpoint_a[1])) + mapped_b = (endpoint_b[0], alpha.get(endpoint_b[1], endpoint_b[1])) + if mapped_b < mapped_a: + mapped_a, mapped_b = mapped_b, mapped_a + image.append((mapped_a, mapped_b, order_before, order_after)) + candidate = tuple(sorted(image)) + if best is None or candidate < best: + best = candidate + return best + + +def center_signature(center: frozenset, orbits: dict[int, int], symbols: list[str]) -> tuple: + """ + Build a readable orbit-level signature of a reaction center. Equivalent maps always share a signature, + but distinct maps may also collide, so this is a cheap bucketing invariant rather than a decision + procedure - :func:`canonical_center_key` is the exact test. + + Args: + center (frozenset): The changed-bond set. + orbits (dict[int, int]): Orbit id per reactant core atom. + symbols (list[str]): Element symbol per reactant atom. + + Returns: + tuple: The signature. + """ + entries = list() + for endpoint_a, endpoint_b, order_before, order_after in center: + described = tuple(sorted(f'{"H@" if tag == HYDROGEN else ""}{symbols[index]}{orbits.get(index, index)}' + for tag, index in (endpoint_a, endpoint_b))) + entries.append((described[0], described[1], order_before, order_after)) + return tuple(sorted(entries)) + + +def enumerate_atom_maps(rxn, + backend: str = 'ARC', + include_flipped: bool = True, + max_maps: int = MAX_ENUMERATED_MAPS, + ) -> list[list[int]]: + """ + Enumerate valid atom maps for a reaction by sweeping every RMG template ``product_dict`` and, optionally, + both reaction directions. Duplicate maps are removed while preserving discovery order. + + Multiplicity enters from two independent sources, both of which the single-map path discards: + the template ``product_dict`` sweep, and - within one product dictionary - the several superimposable + backbone maps per scissored fragment pair, combined by ``map_rxn_all``. + + Args: + rxn (ARCReaction): The reaction to map. + backend (str, optional): Currently only ``ARC``'s method is implemented as the backend. + include_flipped (bool, optional): Whether to also map the flipped reaction and un-flip the results. + max_maps (int, optional): Stop after this many distinct maps. + + Returns: + list[list[int]]: The distinct atom maps found. + """ + maps: list[list[int]] = list() + seen: set[tuple] = set() + + def collect(candidates: list[list[int]], flip: bool = False) -> None: + for atom_map in candidates: + if len(maps) >= max_maps: + return + if flip: + atom_map = flip_map(atom_map) + if atom_map is None: + continue + key = tuple(atom_map) + if key not in seen: + seen.add(key) + maps.append(list(atom_map)) + + def sweep(target, flip: bool) -> None: + n_product_dicts = len(target.product_dicts) if getattr(target, 'product_dicts', None) else 1 + for pdi in range(min(n_product_dicts, MAX_PDI)): + try: + collect(map_rxn_all(target, backend=backend, product_dict_index_to_try=pdi), flip=flip) + except (ValueError, IndexError, KeyError) as e: + logger.debug(f'enumerate_atom_maps: product_dict {pdi} of {target} ' + f'(flip={flip}) failed with {e!r}.') + if len(maps) >= max_maps: + return + + sweep(rxn, flip=False) + if include_flipped and len(maps) < max_maps: + sweep(rxn.flip_reaction(), flip=True) + + if len(maps) >= max_maps: + logger.warning(f'enumerate_atom_maps hit the cap of {max_maps} maps for {rxn}; ' + f'reported degeneracies are lower bounds.') + return maps + + +def cluster_atom_maps(atom_maps: list[list[int]], + rxn, + ignore_bond_orders: bool = True, + validate_centers: bool = True, + ) -> list[MapCluster]: + """ + Group atom maps into equivalence classes, one per distinct reaction channel. + + Two maps land in the same cluster exactly when their changed-bond sets lie in the same orbit of the + reactant complex automorphism group - see the module docstring for why this is equivalent to the double + coset ``Aut(P) . sigma . Aut(R)``. + + Args: + atom_maps (list[list[int]]): The atom maps to cluster. + rxn (ARCReaction): The reaction the maps belong to. + ignore_bond_orders (bool, optional): Passed through to :func:`changed_bonds`. + validate_centers (bool, optional): Whether to discard maps whose reaction center is invalid. The + family recipe is used as an absolute reference where available + (:func:`expected_reaction_centers`), otherwise the relative + minimal-center heuristic (:func:`filter_minimal_centers`). + + Returns: + list[MapCluster]: The clusters, ordered by descending degeneracy. + """ + if not atom_maps: + return list() + reactants, products = rxn.get_reactants_and_products(return_copies=True) + r_graph, p_graph = build_complex_graph(reactants), build_complex_graph(products) + # The automorphism group and the changed-bond set must agree on bond orders, otherwise a Kekule + # structure makes the group too small and splits clusters that are genuinely equivalent. + automorphisms, truncated = core_automorphisms(r_graph, ignore_bond_orders=ignore_bond_orders) + orbits = core_orbits(r_graph, automorphisms) + + scored: list[tuple[list[int], frozenset]] = list() + 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))) + if validate_centers and scored: + expected = expected_reaction_centers(rxn, ignore_bond_orders=ignore_bond_orders) + validated = filter_expected_centers(scored, expected, rxn) if expected else list() + if validated: + scored = validated + else: + # Either no recipe reference was available, or it rejected every enumerated map - which would + # mean returning no channels at all. Fall back to the relative filter rather than that. + if expected: + logger.warning(f'The {rxn.family} recipe rejected all {len(scored)} enumerated atom maps ' + f'for {rxn}; falling back to the minimal-center filter.') + scored = filter_minimal_centers(scored, rxn) + + clusters: dict[tuple, MapCluster] = dict() + for atom_map, center in scored: + if not center: + logger.debug(f'An atom map for {rxn} induces no bond changes; clustering it under an empty key.') + collapsed = collapse_hydrogens(center, r_graph) + key = canonical_center_key(collapsed, automorphisms) + if key not in clusters: + clusters[key] = MapCluster(representative=list(atom_map), + key=key, + signature=center_signature(collapsed, orbits, r_graph.symbols), + truncated=truncated) + clusters[key].members.append(list(atom_map)) + clusters[key].centers.add(center) + return sorted(clusters.values(), key=lambda cluster: (-cluster.degeneracy, cluster.key)) + + +def expected_reaction_centers(rxn, ignore_bond_orders: bool = True) -> set[frozenset] | None: + """ + The reaction centers predicted by the RMG family recipe, one per template product dictionary, expressed + in reactant-complex indices. + + This is an *absolute* reference for map validity, unlike :func:`filter_minimal_centers` which can only + compare enumerated maps against each other. Each product dictionary's ``r_label_map`` assigns the + family's labelled atoms (``*1``, ``*2``, ...) to concrete reactant indices, and + ``ARCReaction.get_expected_changing_bonds`` turns the family's ``BREAK_BOND`` and ``FORM_BOND`` actions + into index pairs. For CH4 + OH the four product dictionaries yield exactly the four abstraction centers. + + Note the indices in ``r_label_map`` are 0-indexed, matching ``atom_map``, despite the docstring of + ``find_all_breaking_bonds`` describing its own return value as 1-indexed. + + A product dictionary is skipped, rather than the whole reference being abandoned, when its recipe cannot + be read. Two causes are common: + + - **A missing label.** ``get_expected_changing_bonds`` reads a single family-level ``actions`` list, but + a family such as ``intra_H_migration`` spans several template variants with different label sets, so + some product dictionaries lack a label the actions reference and raise ``KeyError``. Only 6 of the 32 + product dictionaries of one benchmark reaction are readable this way. + - **A degenerate self-bond.** An action naming the same label twice (``R_Recombination``'s + ``* + * -> *-*``) resolves both endpoints to one index. ``find_all_breaking_bonds`` disambiguates + those through suffixed label keys; that logic is not reproduced here, so the affected dictionary is + skipped instead of contributing a wrong center. + + Returns ``None``, meaning "no usable reference at all", when the reaction has no family or no product + dictionaries, when ``ignore_bond_orders`` is ``False`` (the recipe describes only bond breaking and + formation and cannot predict the pure order changes that :func:`changed_bonds` would then report), or + when every product dictionary was skipped. + + Args: + rxn (ARCReaction): The reaction. + ignore_bond_orders (bool, optional): Must match the setting used by :func:`changed_bonds`. + + Returns: + set[frozenset] | None: The predicted centers, or ``None`` if no reliable reference is available. + """ + if not ignore_bond_orders or rxn.family is None or not getattr(rxn, 'product_dicts', None): + return None + centers, skipped = set(), 0 + for product_dict in rxn.product_dicts[:MAX_PDI]: + r_label_map = product_dict.get('r_label_map') + if not r_label_map: + skipped += 1 + continue + try: + breaking, forming = rxn.get_expected_changing_bonds(r_label_dict=r_label_map) + except (KeyError, TypeError, ValueError): + skipped += 1 + continue + if breaking is None and forming is None: + skipped += 1 + continue + center = set() + for pairs, orders in ((breaking or list(), (1, 0)), (forming or list(), (0, 1))): + for i, j in pairs: + if i == j: + center = None + break + center.add((min(i, j), max(i, j)) + orders) + if center is None: + break + if not center: + skipped += 1 + continue + centers.add(frozenset(center)) + if skipped: + logger.debug(f'expected_reaction_centers: skipped {skipped} unreadable product dictionaries of ' + f'{rxn} ({rxn.family}), kept {len(centers)} predicted centers.') + return centers or None + + +def filter_expected_centers(scored: list[tuple[list[int], frozenset]], + expected: set[frozenset], + rxn=None, + ) -> list[tuple[list[int], frozenset]]: + """ + Keep only the maps whose reaction center is one the family recipe actually predicts. + + Args: + scored (list[tuple[list[int], frozenset]]): ``(atom_map, changed_bond_set)`` pairs. + expected (set[frozenset]): The centers from :func:`expected_reaction_centers`. + rxn (ARCReaction, optional): The reaction, used only for logging. + + Returns: + list[tuple[list[int], frozenset]]: The pairs whose center matches the recipe. + """ + kept = [entry for entry in scored if entry[1] in expected] + if len(kept) != len(scored): + logger.debug(f'Discarded {len(scored) - len(kept)} of {len(scored)} atom maps for {rxn} whose ' + f'reaction center is not predicted by the {getattr(rxn, "family", None)} recipe.') + return kept + + +def filter_minimal_centers(scored: list[tuple[list[int], frozenset]], + rxn=None, + ) -> list[tuple[list[int], frozenset]]: + """ + Keep only the maps whose reaction center is as small as the smallest one observed. + + This is a validity filter, and it is needed because :func:`changed_bonds` computes a changed-bond set + for *any* element-preserving bijection - there is nothing in the arithmetic that distinguishes a correct + atom map from a scrambled one. ``map_rxn_all`` deliberately keeps every graph-superimposable backbone + candidate per fragment, including those that ``map_two_species`` would have rejected on RMSD, so + scrambled maps do reach this point. Left unfiltered they register as extra "channels": for the + Diels-Alder of butadiene with ethene the correct map changes 2 bonds while the scrambled ones change 10 + to 14, and each scrambled variant would otherwise open its own cluster. + + An elementary reaction rearranges a minimal set of bonds, so the smallest observed center is the correct + one and anything larger is a mapping error rather than a distinct channel. + + Note this is a heuristic over the enumerated set, not an independent check against the family recipe. It + assumes at least one correct map was enumerated; if every enumerated map is wrong, the least wrong one + survives. ``ARCReaction.get_expected_changing_bonds`` would give an absolute reference and is the natural + next step. + + Args: + scored (list[tuple[list[int], frozenset]]): ``(atom_map, changed_bond_set)`` pairs. + rxn (ARCReaction, optional): The reaction, used only for logging. + + Returns: + list[tuple[list[int], frozenset]]: The pairs whose center is of minimal size. + """ + if not scored: + return scored + smallest = min(len(center) for _, center in scored) + kept = [entry for entry in scored if len(entry[1]) == smallest] + if len(kept) != len(scored): + discarded = sorted({len(center) for _, center in scored if len(center) != smallest}) + logger.debug(f'Discarded {len(scored) - len(kept)} of {len(scored)} atom maps for {rxn} whose ' + f'reaction center ({discarded} changed bonds) exceeds the minimal {smallest}.') + return kept + + +def map_reaction_clusters(rxn, + backend: str = 'ARC', + include_flipped: bool = True, + ignore_bond_orders: bool = True, + ) -> list[MapCluster]: + """ + Convenience wrapper: enumerate every atom map for a reaction and cluster them into distinct channels. + + Args: + rxn (ARCReaction): The reaction to map. + backend (str, optional): Currently only ``ARC``'s method is implemented as the backend. + include_flipped (bool, optional): Whether to also map the flipped reaction. + ignore_bond_orders (bool, optional): Passed through to :func:`changed_bonds`. + + Returns: + list[MapCluster]: The clusters, ordered by descending degeneracy. + """ + return cluster_atom_maps(enumerate_atom_maps(rxn, backend=backend, include_flipped=include_flipped), + rxn, + ignore_bond_orders=ignore_bond_orders) diff --git a/arc/mapping/cluster_test.py b/arc/mapping/cluster_test.py new file mode 100644 index 0000000000..536e74dc5a --- /dev/null +++ b/arc/mapping/cluster_test.py @@ -0,0 +1,562 @@ +#!/usr/bin/env python3 +# encoding: utf-8 + +""" +This module contains unit tests of the arc.mapping.cluster module +""" + +import unittest + +import arc.mapping.cluster as cluster +from arc.reaction import ARCReaction +from arc.species import ARCSpecies + + +class TestComplexGraph(unittest.TestCase): + """ + Contains unit tests for building the flat complex graph. + """ + + def test_build_complex_graph_methane(self): + """Test that a single species is flattened correctly, with hydrogens attached to their parent.""" + graph = cluster.build_complex_graph([ARCSpecies(label='CH4', smiles='C')]) + self.assertEqual(graph.n_atoms, 5) + self.assertEqual(graph.symbols, ['C', 'H', 'H', 'H', 'H']) + self.assertEqual(graph.core, [0]) + self.assertEqual(graph.parent, {1: 0, 2: 0, 3: 0, 4: 0}) + self.assertEqual(len(graph.bonds), 4) + # The single core atom has no core neighbors, but carries four hydrogens in its invariant. + self.assertEqual(graph.adj[0], {}) + self.assertEqual(graph.invariant[0][-1], 4) + + def test_build_complex_graph_uses_running_indices(self): + """Test that atom indices continue across species, matching the atom_map convention.""" + graph = cluster.build_complex_graph([ARCSpecies(label='CH4', smiles='C'), + ARCSpecies(label='OH', smiles='[OH]')]) + self.assertEqual(graph.n_atoms, 7) + self.assertEqual(graph.symbols, ['C', 'H', 'H', 'H', 'H', 'O', 'H']) + self.assertEqual(graph.core, [0, 5]) + self.assertEqual(graph.parent, {1: 0, 2: 0, 3: 0, 4: 0, 6: 5}) + # The O-H bond of the second species must be offset onto the running indices. + self.assertIn(frozenset((5, 6)), graph.bonds) + + def test_build_complex_graph_promotes_parentless_hydrogens(self): + """Test that a hydrogen without a unique heavy neighbor is promoted to a core atom.""" + h2 = cluster.build_complex_graph([ARCSpecies(label='H2', smiles='[H][H]')]) + self.assertEqual(h2.core, [0, 1]) + self.assertEqual(h2.parent, {}) + h_atom = cluster.build_complex_graph([ARCSpecies(label='H', smiles='[H]')]) + self.assertEqual(h_atom.core, [0]) + + def test_build_complex_graph_core_adjacency(self): + """Test that only core-core bonds enter the adjacency, with their bond orders.""" + graph = cluster.build_complex_graph([ARCSpecies(label='ethene', smiles='C=C')]) + self.assertEqual(graph.core, [0, 1]) + self.assertEqual(graph.adj[0], {1: 2}) + self.assertEqual(graph.adj[1], {0: 2}) + + +class TestAutomorphisms(unittest.TestCase): + """ + Contains unit tests for the core-skeleton automorphism machinery. + """ + + @staticmethod + def _automorphism_count(smiles_list, ignore_bond_orders=True): + """Return |Aut| of the core skeleton of a complex built from a list of SMILES.""" + graph = cluster.build_complex_graph([ARCSpecies(label=f's{i}', smiles=smiles) + for i, smiles in enumerate(smiles_list)]) + automorphisms, truncated = cluster.core_automorphisms(graph, ignore_bond_orders=ignore_bond_orders) + return len(automorphisms), truncated + + def test_core_automorphism_counts(self): + """Test |Aut| of the core skeleton against known graph automorphism group orders.""" + for smiles, expected in [('C', 1), # methane, one core atom + ('CC', 2), # ethane, swap the two carbons + ('CCC', 2), # propane, reflect the chain + ('c1ccccc1', 12), # benzene, the dihedral group D6 + ('CC(C)C', 6), # isobutane, permute three methyls + ('CC(C)(C)C', 24), # neopentane, permute four methyls + ('c1ccc2ccccc2c1', 4)]: # naphthalene, Z2 x Z2 + count, truncated = self._automorphism_count([smiles]) + self.assertEqual(count, expected, msg=f'|Aut| of {smiles}') + self.assertFalse(truncated) + + def test_core_automorphisms_include_molecule_swap(self): + """Test that exchanging two identical molecules counts as a symmetry of the complex.""" + self.assertEqual(self._automorphism_count(['O[O]', 'O[O]'])[0], 2) + # Two benzenes: 12 per ring, times the swap. + self.assertEqual(self._automorphism_count(['c1ccccc1', 'c1ccccc1'])[0], 12 * 12 * 2) + + def test_core_automorphisms_bond_order_sensitivity(self): + """Test that honoring Kekule bond orders under-counts an aromatic ring's symmetry.""" + # Only the 6 automorphisms preserving the single/double alternation survive. + self.assertEqual(self._automorphism_count(['c1ccccc1'], ignore_bond_orders=False)[0], 6) + self.assertEqual(self._automorphism_count(['c1ccccc1'], ignore_bond_orders=True)[0], 12) + # A saturated species has no multiple bonds, so the setting cannot matter. + self.assertEqual(self._automorphism_count(['CC(C)(C)C'], ignore_bond_orders=False)[0], + self._automorphism_count(['CC(C)(C)C'], ignore_bond_orders=True)[0]) + + def test_bond_orders_are_ignored_by_default(self): + """Test that the default settings ignore bond orders, everywhere the choice is made. + + The explicit-flag tests above would still pass if a default were wired the wrong way round, so this + pins the defaults themselves. Getting one of them wrong reintroduces the Kekule bug: an aromatic + ring loses half its automorphisms and its clusters over-split. + """ + benzene = cluster.build_complex_graph([ARCSpecies(label='benzene', smiles='c1ccccc1')]) + self.assertEqual(len(cluster.core_automorphisms(benzene)[0]), 12) + ethene = cluster.build_complex_graph([ARCSpecies(label='ethene', smiles='C=C')]) + self.assertEqual(cluster.effective_adjacency(ethene)[0], {1: 1}) + # Colour refinement must see the ring as uniform rather than as alternating bond orders. + self.assertEqual(len(set(cluster.refine_colors(benzene).values())), 1) + + def test_core_automorphisms_are_genuine(self): + """Test that every returned permutation really preserves the core adjacency.""" + graph = cluster.build_complex_graph([ARCSpecies(label='isobutane', smiles='CC(C)C')]) + automorphisms, _ = cluster.core_automorphisms(graph) + adjacency = cluster.effective_adjacency(graph, ignore_bond_orders=True) + for alpha in automorphisms: + self.assertEqual(sorted(alpha.keys()), sorted(graph.core)) + self.assertEqual(sorted(alpha.values()), sorted(graph.core)) + for v in graph.core: + for u in graph.core: + self.assertEqual(adjacency[v].get(u), adjacency[alpha[v]].get(alpha[u])) + + def test_core_automorphisms_truncation_is_reported(self): + """Test that hitting the cap is reported rather than applied silently.""" + graph = cluster.build_complex_graph([ARCSpecies(label='neopentane', smiles='CC(C)(C)C')]) + automorphisms, truncated = cluster.core_automorphisms(graph, max_count=5) + self.assertTrue(truncated) + self.assertLessEqual(len(automorphisms), 5) + + def test_core_orbits(self): + """Test that orbits identify the symmetry-equivalent atoms of propane.""" + graph = cluster.build_complex_graph([ARCSpecies(label='propane', smiles='CCC')]) + automorphisms, _ = cluster.core_automorphisms(graph) + orbits = cluster.core_orbits(graph, automorphisms) + # The two terminal carbons share an orbit, the central one is on its own. + self.assertEqual(orbits[0], orbits[2]) + self.assertNotEqual(orbits[0], orbits[1]) + self.assertEqual(len(set(orbits.values())), 2) + + def test_effective_adjacency_normalizes_orders(self): + """Test that ignoring bond orders replaces every order by 1 without changing connectivity.""" + graph = cluster.build_complex_graph([ARCSpecies(label='ethene', smiles='C=C')]) + honored = cluster.effective_adjacency(graph, ignore_bond_orders=False) + ignored = cluster.effective_adjacency(graph, ignore_bond_orders=True) + self.assertEqual(honored[0], {1: 2}) + self.assertEqual(ignored[0], {1: 1}) + self.assertEqual(set(honored.keys()), set(ignored.keys())) + + def test_refine_colors_separates_inequivalent_atoms(self): + """Test that colour refinement distinguishes propane's primary and secondary carbons.""" + graph = cluster.build_complex_graph([ARCSpecies(label='propane', smiles='CCC')]) + colors = cluster.refine_colors(graph) + self.assertEqual(colors[0], colors[2]) + self.assertNotEqual(colors[0], colors[1]) + + +class TestReactionCenters(unittest.TestCase): + """ + Contains unit tests for changed-bond sets, their hydrogen collapse, and their canonical form. + """ + + @classmethod + def setUpClass(cls): + """A method that is run before all unit tests in this class.""" + cls.maxDiff = None + # Reactants CH4 + OH: 0=C, 1-4=H, 5=O, 6=H. Products CH3 + H2O: 0=C, 1-3=H, 4=O, 5-6=H. + cls.r_graph = cluster.build_complex_graph([ARCSpecies(label='CH4', smiles='C'), + ARCSpecies(label='OH', smiles='[OH]')]) + cls.p_graph = cluster.build_complex_graph([ARCSpecies(label='CH3', smiles='[CH3]'), + ARCSpecies(label='H2O', smiles='O')]) + # Abstract reactant hydrogen 1; it becomes product hydrogen 6 of water. + cls.abstraction_map = [0, 6, 1, 2, 3, 4, 5] + + def test_changed_bonds_h_abstraction(self): + """Test that an abstraction map yields exactly one broken and one formed bond.""" + center = cluster.changed_bonds(self.r_graph, self.p_graph, self.abstraction_map) + self.assertEqual(len(center), 2) + # The C-H bond to reactant hydrogen 1 breaks, and that hydrogen bonds to the oxygen. + self.assertIn((0, 1, 1, 0), center) + self.assertIn((1, 5, 0, 1), center) + + def test_changed_bonds_are_expressed_in_reactant_indices(self): + """Test that every changed-bond endpoint is a valid reactant atom index.""" + center = cluster.changed_bonds(self.r_graph, self.p_graph, self.abstraction_map) + for i, j, _, _ in center: + self.assertLess(i, j) + self.assertLess(j, self.r_graph.n_atoms) + + def test_changed_bonds_invariant_to_product_relabeling(self): + """Test that relabeling product atoms by a product automorphism leaves the center unchanged. + + This is the property that makes the changed-bond set the right unit of degeneracy: swapping the two + equivalent hydrogens of the product water is an element of Aut(P) and is not a different path. + """ + swapped = [5 if index == 6 else 6 if index == 5 else index for index in self.abstraction_map] + self.assertNotEqual(swapped, self.abstraction_map) + self.assertEqual(cluster.changed_bonds(self.r_graph, self.p_graph, self.abstraction_map), + cluster.changed_bonds(self.r_graph, self.p_graph, swapped)) + + def test_changed_bonds_identity_is_empty(self): + """Test that mapping a complex onto an identical complex changes no bonds.""" + identity = list(range(self.r_graph.n_atoms)) + self.assertEqual(cluster.changed_bonds(self.r_graph, self.r_graph, identity), frozenset()) + + def test_collapse_hydrogens(self): + """Test that hydrogens are replaced by their parent core atom, and core atoms are left alone.""" + center = cluster.changed_bonds(self.r_graph, self.p_graph, self.abstraction_map) + collapsed = cluster.collapse_hydrogens(center, self.r_graph) + self.assertEqual(len(collapsed), 2) + # Reactant hydrogen 1 hangs off carbon 0, the oxygen is core index 5. + self.assertIn(((cluster.CORE, 0), (cluster.HYDROGEN, 0), 1, 0), collapsed) + self.assertIn(((cluster.CORE, 5), (cluster.HYDROGEN, 0), 0, 1), collapsed) + + def test_collapse_hydrogens_merges_equivalent_hydrogens(self): + """Test that abstracting any hydrogen of methane gives the same collapsed center. + + This loss is intended - it is what places all four abstractions in one cluster - and it is why + degeneracy is counted on the uncollapsed centers instead. + """ + collapsed = set() + for hydrogen in (1, 2, 3, 4): + center = frozenset({(0, hydrogen, 1, 0)}) + collapsed.add(cluster.collapse_hydrogens(center, self.r_graph)) + self.assertEqual(len(collapsed), 1) + + +class TestCanonicalKey(unittest.TestCase): + """ + Contains unit tests for canonicalizing a reaction center under the reactant automorphism group. + """ + + @classmethod + def setUpClass(cls): + """A method that is run before all unit tests in this class.""" + cls.graph = cluster.build_complex_graph([ARCSpecies(label='propane', smiles='CCC')]) + cls.automorphisms, _ = cluster.core_automorphisms(cls.graph) + cls.hydrogens = dict() + for hydrogen, parent in cls.graph.parent.items(): + cls.hydrogens.setdefault(parent, list()).append(hydrogen) + + def _key_for_abstraction_at(self, carbon): + """Build a minimal abstraction center at ``carbon`` and canonicalize it.""" + hydrogen = self.hydrogens[carbon][0] + center = frozenset({(min(carbon, hydrogen), max(carbon, hydrogen), 1, 0)}) + return cluster.canonical_center_key(cluster.collapse_hydrogens(center, self.graph), + self.automorphisms) + + def test_equivalent_sites_share_a_key(self): + """Test that abstraction at either primary carbon of propane gives the same key.""" + self.assertEqual(self._key_for_abstraction_at(0), self._key_for_abstraction_at(2)) + + def test_inequivalent_sites_have_different_keys(self): + """Test that abstraction at a primary and at the secondary carbon give different keys.""" + self.assertNotEqual(self._key_for_abstraction_at(0), self._key_for_abstraction_at(1)) + + def test_key_is_hashable_and_stable(self): + """Test that the key is hashable and does not depend on which equivalent site produced it.""" + keys = {self._key_for_abstraction_at(0), self._key_for_abstraction_at(2)} + self.assertEqual(len(keys), 1) + + def test_key_without_automorphisms_is_the_sorted_center(self): + """Test that an empty group leaves the center sorted but otherwise untouched.""" + center = frozenset({((cluster.CORE, 1), (cluster.HYDROGEN, 0), 1, 0)}) + self.assertEqual(cluster.canonical_center_key(center, list()), tuple(sorted(center))) + + +class TestMinimalCenterFilter(unittest.TestCase): + """ + Contains unit tests for the validity filter that discards scrambled atom maps. + """ + + def test_keeps_only_the_smallest_centers(self): + """Test that maps with a larger reaction center than the minimum are discarded.""" + small_a = ([0, 1], frozenset({(0, 1, 1, 0)})) + small_b = ([1, 0], frozenset({(0, 2, 1, 0)})) + large = ([0, 1], frozenset({(0, 1, 1, 0), (1, 2, 0, 1), (2, 3, 1, 0)})) + kept = cluster.filter_minimal_centers([small_a, large, small_b]) + self.assertEqual(kept, [small_a, small_b]) + + def test_keeps_everything_when_all_centers_match(self): + """Test that the filter is a no-op when every center is already of minimal size.""" + scored = [([0, 1], frozenset({(0, 1, 1, 0)})), ([1, 0], frozenset({(0, 2, 1, 0)}))] + self.assertEqual(cluster.filter_minimal_centers(scored), scored) + + def test_empty_input(self): + """Test that an empty input is passed through.""" + self.assertEqual(cluster.filter_minimal_centers(list()), list()) + + +class _StubReaction(object): + """A minimal stand-in exposing only what expected_reaction_centers touches.""" + + def __init__(self, family=None, product_dicts=None, breaking=None, forming=None): + self.family = family + self.product_dicts = product_dicts + self._breaking = breaking + self._forming = forming + + def get_expected_changing_bonds(self, r_label_dict, family=None): + """Return the canned recipe bonds.""" + return self._breaking, self._forming + + +class TestExpectedCenters(unittest.TestCase): + """ + Contains unit tests for the family-recipe validity reference. + """ + + def test_returns_none_without_a_family(self): + """Test that a reaction with no family gives no reference.""" + self.assertIsNone(cluster.expected_reaction_centers(_StubReaction(family=None))) + + def test_returns_none_without_product_dicts(self): + """Test that a reaction with no template product dictionaries gives no reference.""" + self.assertIsNone(cluster.expected_reaction_centers(_StubReaction(family='H_Abstraction', + product_dicts=list()))) + + def test_returns_none_when_bond_orders_are_honored(self): + """Test that the recipe is refused when it cannot describe what changed_bonds would report. + + The recipe only names bond breaking and formation, so it cannot predict the pure order changes that + changed_bonds reports when bond orders are honored. + """ + rxn = _StubReaction(family='H_Abstraction', product_dicts=[{'r_label_map': {'*1': 0}}], + breaking=[(0, 1)], forming=[(1, 5)]) + self.assertIsNone(cluster.expected_reaction_centers(rxn, ignore_bond_orders=False)) + self.assertIsNotNone(cluster.expected_reaction_centers(rxn, ignore_bond_orders=True)) + + def test_returns_none_on_a_degenerate_self_bond(self): + """Test that a same-label recipe action is rejected rather than used as a wrong reference.""" + rxn = _StubReaction(family='R_Recombination', product_dicts=[{'r_label_map': {'*': 3}}], + breaking=list(), forming=[(3, 3)]) + self.assertIsNone(cluster.expected_reaction_centers(rxn)) + + def test_unreadable_product_dicts_are_skipped_not_fatal(self): + """Test that one unreadable product dictionary does not discard the whole reference. + + A family such as intra_H_migration spans template variants with different label sets, so some + product dictionaries raise KeyError while others are perfectly readable. Abandoning the reference on + the first failure would silently drop the reaction back to the relative filter. + """ + + class _PartialStub(_StubReaction): + def get_expected_changing_bonds(self, r_label_dict, family=None): + if 'good' not in r_label_dict: + raise KeyError('*6') + return [(0, 1)], [(1, 5)] + + rxn = _PartialStub(family='intra_H_migration', + product_dicts=[{'r_label_map': {'bad': 0}}, {'r_label_map': {'good': 0}}]) + self.assertEqual(cluster.expected_reaction_centers(rxn), {frozenset({(0, 1, 1, 0), (1, 5, 0, 1)})}) + + def test_returns_none_when_every_product_dict_is_unreadable(self): + """Test that the reference is only abandoned once nothing at all could be read.""" + + class _BrokenStub(_StubReaction): + def get_expected_changing_bonds(self, r_label_dict, family=None): + raise KeyError('*6') + + rxn = _BrokenStub(family='intra_H_migration', product_dicts=[{'r_label_map': {'a': 0}}]) + self.assertIsNone(cluster.expected_reaction_centers(rxn)) + + def test_builds_one_center_per_product_dict(self): + """Test that each product dictionary contributes its own predicted center.""" + rxn = _StubReaction(family='H_Abstraction', + product_dicts=[{'r_label_map': {'*1': 0}}, {'r_label_map': {'*1': 0}}], + breaking=[(0, 1)], forming=[(1, 5)]) + centers = cluster.expected_reaction_centers(rxn) + # Both dictionaries yield the same canned bonds here, so they collapse to one center. + self.assertEqual(centers, {frozenset({(0, 1, 1, 0), (1, 5, 0, 1)})}) + + def test_filter_expected_centers(self): + """Test that only maps whose center the recipe predicts are kept.""" + good = ([0, 1], frozenset({(0, 1, 1, 0)})) + bad = ([1, 0], frozenset({(2, 3, 1, 0)})) + kept = cluster.filter_expected_centers([good, bad], {frozenset({(0, 1, 1, 0)})}) + self.assertEqual(kept, [good]) + + def test_filter_expected_centers_can_reject_everything(self): + """Test that the filter reports an empty result rather than silently keeping invalid maps.""" + entry = ([0, 1], frozenset({(9, 9, 1, 0)})) + self.assertEqual(cluster.filter_expected_centers([entry], {frozenset({(0, 1, 1, 0)})}), list()) + + +class TestMapCluster(unittest.TestCase): + """ + Contains unit tests for the MapCluster container. + """ + + def test_degeneracy_counts_distinct_centers_not_maps(self): + """Test that degeneracy counts reaction paths, so an Aut(P) relabeling does not inflate it.""" + map_cluster = cluster.MapCluster(representative=[0, 1, 2]) + map_cluster.members.extend([[0, 1, 2], [0, 2, 1]]) + # Both maps describe the same path, so they share a center. + map_cluster.centers.add(frozenset({(0, 1, 1, 0)})) + self.assertEqual(len(map_cluster.members), 2) + self.assertEqual(map_cluster.degeneracy, 1) + # A genuinely different path adds a center. + map_cluster.centers.add(frozenset({(0, 2, 1, 0)})) + self.assertEqual(map_cluster.degeneracy, 2) + + def test_empty_cluster_has_zero_degeneracy(self): + """Test the degenerate case of a cluster with no recorded center.""" + self.assertEqual(cluster.MapCluster(representative=[0]).degeneracy, 0) + + +class TestClusteringIntegration(unittest.TestCase): + """ + Contains end-to-end unit tests running the full enumerate-and-cluster pipeline on a real reaction. + """ + + @classmethod + def setUpClass(cls): + """A method that is run before all unit tests in this class.""" + cls.maxDiff = None + ch4_xyz = {'symbols': ('C', 'H', 'H', 'H', 'H'), 'isotopes': (12, 1, 1, 1, 1), + 'coords': ((-5.45906343962835e-10, 4.233517924761169e-10, 2.9505240956083194e-10), + (-0.6505520089868748, -0.7742801979689132, -0.4125187934483119), + (-0.34927557824779626, 0.9815958255612931, -0.3276823191685369), + (-0.022337921721882443, -0.04887374527620588, 1.0908766524267022), + (1.0221655095024578, -0.15844188273952128, -0.350675540104908))} + oh_xyz = """O 0.48890387 0.00000000 0.00000000 + H -0.48890387 0.00000000 0.00000000""" + ch3_xyz = """C 0.00000000 0.00000001 -0.00000000 + H 1.06690511 -0.17519582 0.05416493 + H -0.68531716 -0.83753536 -0.02808565 + H -0.38158795 1.01273118 -0.02607927""" + h2o_xyz = """O -0.00032832 0.39781490 0.00000000 + H -0.76330345 -0.19953755 0.00000000 + H 0.76363177 -0.19827735 0.00000000""" + cls.rxn = ARCReaction(r_species=[ARCSpecies(label='CH4', smiles='C', xyz=ch4_xyz), + ARCSpecies(label='OH', smiles='[OH]', xyz=oh_xyz)], + p_species=[ARCSpecies(label='CH3', smiles='[CH3]', xyz=ch3_xyz), + ARCSpecies(label='H2O', smiles='O', xyz=h2o_xyz)]) + cls.atom_maps = cluster.enumerate_atom_maps(cls.rxn) + + def test_enumerate_atom_maps_returns_several_valid_permutations(self): + """Test that enumeration finds more than one map and that each is a permutation of the atoms.""" + self.assertGreater(len(self.atom_maps), 1) + for atom_map in self.atom_maps: + self.assertEqual(sorted(atom_map), list(range(7))) + + def test_enumerate_atom_maps_are_distinct(self): + """Test that the enumeration does not report the same map twice.""" + self.assertEqual(len({tuple(atom_map) for atom_map in self.atom_maps}), len(self.atom_maps)) + + def test_enumerate_atom_maps_honors_the_cap(self): + """Test that the cap on the number of enumerated maps is respected.""" + self.assertLessEqual(len(cluster.enumerate_atom_maps(self.rxn, max_maps=2)), 2) + + def test_cluster_atom_maps_finds_a_single_channel(self): + """Test that abstracting any of methane's four hydrogens is recognized as one channel.""" + clusters = cluster.cluster_atom_maps(self.atom_maps, self.rxn) + self.assertEqual(len(clusters), 1) + + def test_cluster_degeneracy_matches_the_reaction_path_degeneracy(self): + """Test that the reaction path degeneracy of CH4 + OH abstraction is recovered as 4.""" + clusters = cluster.cluster_atom_maps(self.atom_maps, self.rxn) + self.assertEqual(clusters[0].degeneracy, 4) + self.assertFalse(clusters[0].truncated) + + def test_cluster_signature_describes_the_abstraction(self): + """Test that the signature reports one C-H bond breaking and one O-H bond forming.""" + clusters = cluster.cluster_atom_maps(self.atom_maps, self.rxn) + signature = clusters[0].signature + self.assertEqual(len(signature), 2) + broken = [entry for entry in signature if (entry[2], entry[3]) == (1, 0)] + formed = [entry for entry in signature if (entry[2], entry[3]) == (0, 1)] + self.assertEqual(len(broken), 1) + self.assertEqual(len(formed), 1) + self.assertIn('C', broken[0][0] + broken[0][1]) + self.assertIn('O', formed[0][0] + formed[0][1]) + + def test_cluster_representative_is_a_member(self): + """Test that a cluster's representative is one of its own members.""" + for map_cluster in cluster.cluster_atom_maps(self.atom_maps, self.rxn): + self.assertIn(map_cluster.representative, map_cluster.members) + + def test_cluster_atom_maps_with_no_maps(self): + """Test that clustering an empty list returns no clusters.""" + self.assertEqual(cluster.cluster_atom_maps(list(), self.rxn), list()) + + def test_cluster_atom_maps_skips_wrongly_sized_maps(self): + """Test that a map whose length does not match the reactant complex is skipped.""" + self.assertEqual(cluster.cluster_atom_maps([[0, 1, 2]], self.rxn), list()) + + def test_map_reaction_clusters_wrapper(self): + """Test that the convenience wrapper reproduces enumerate followed by cluster.""" + clusters = cluster.map_reaction_clusters(self.rxn) + self.assertEqual(len(clusters), 1) + self.assertEqual(clusters[0].degeneracy, 4) + + def test_expected_reaction_centers_matches_the_product_dicts(self): + """Test that the family recipe predicts one abstraction center per methane hydrogen. + + Reactant frame: 0=C, 1-4=H of methane, 5=O, 6=H of the hydroxyl. Abstracting hydrogen h breaks + (0, h) and forms (h, 5). + """ + centers = cluster.expected_reaction_centers(self.rxn) + self.assertEqual(len(centers), 4) + for hydrogen in (1, 2, 3, 4): + self.assertIn(frozenset({(0, hydrogen, 1, 0), (hydrogen, 5, 0, 1)}), centers) + + def test_clustered_maps_all_have_a_recipe_predicted_center(self): + """Test that every map surviving validation has a center the family recipe predicts. + + This pins that the absolute reference is what actually admitted the maps, rather than the relative + minimal-center fallback silently doing the work. + """ + expected = cluster.expected_reaction_centers(self.rxn) + self.assertTrue(expected) + clusters = cluster.cluster_atom_maps(self.atom_maps, self.rxn) + for map_cluster in clusters: + for center in map_cluster.centers: + self.assertIn(center, expected) + + def test_arc_reaction_atom_map_clusters_property(self): + """Test that ARCReaction exposes the clusters and caches them.""" + clusters = self.rxn.atom_map_clusters + self.assertEqual(len(clusters), 1) + self.assertEqual(clusters[0].degeneracy, 4) + self.assertIn(clusters[0].representative, clusters[0].members) + # A second access must reuse the cached result rather than recomputing. + self.assertIs(self.rxn.atom_map_clusters, clusters) + + def test_arc_reaction_atom_map_does_not_trigger_clustering(self): + """Test that the cheap atom_map property never pays for the expensive enumeration.""" + rxn = ARCReaction(r_species=[spc.copy() for spc in self.rxn.r_species], + p_species=[spc.copy() for spc in self.rxn.p_species]) + self.assertIsNone(rxn._atom_map_clusters) + self.assertIsNotNone(rxn.atom_map) + self.assertIsNone(rxn._atom_map_clusters) + + def test_arc_reaction_atom_map_clusters_is_settable(self): + """Test that the clusters can be set, and reset to None to force a recomputation.""" + rxn = ARCReaction(r_species=[spc.copy() for spc in self.rxn.r_species], + p_species=[spc.copy() for spc in self.rxn.p_species]) + rxn.atom_map_clusters = ['sentinel'] + self.assertEqual(rxn.atom_map_clusters, ['sentinel']) + rxn.atom_map_clusters = None + self.assertIsNone(rxn._atom_map_clusters) + + def test_arc_reaction_atom_map_clusters_not_persisted(self): + """Test that as_dict carries the atom map but never the derived clusters.""" + rxn = ARCReaction(r_species=[spc.copy() for spc in self.rxn.r_species], + p_species=[spc.copy() for spc in self.rxn.p_species]) + self.assertEqual(len(rxn.atom_map_clusters), 1) + # Clustering must not populate the atom map as a side effect, so as_dict still omits it. + self.assertIsNone(rxn._atom_map) + self.assertNotIn('atom_map', rxn.as_dict()) + # Once the atom map is computed it is persisted, but the clusters never are. + self.assertIsNotNone(rxn.atom_map) + reaction_dict = rxn.as_dict() + self.assertIn('atom_map', reaction_dict) + self.assertNotIn('atom_map_clusters', reaction_dict) + + +if __name__ == '__main__': + unittest.main(testRunner=unittest.TextTestRunner(verbosity=2)) From 8ec0af12e57f2ad9ace41f635ef587cfd9630720 Mon Sep 17 00:00:00 2001 From: kfir4444 Date: Sun, 23 Aug 2026 18:33:45 +0300 Subject: [PATCH 5/8] Expose the distinct reaction channels as ARCReaction.atom_map_clusters --- arc/reaction/reaction.py | 32 ++++++++++++++++++++++++++++++++ 1 file changed, 32 insertions(+) diff --git a/arc/reaction/reaction.py b/arc/reaction/reaction.py index 3b21baa3ff..29b6faebb7 100644 --- a/arc/reaction/reaction.py +++ b/arc/reaction/reaction.py @@ -12,6 +12,7 @@ translate_xyz, xyz_to_str, ) +from arc.mapping.cluster import map_reaction_clusters from arc.mapping.driver import map_reaction from arc.species.species import ARCSpecies, check_atom_balance, check_label @@ -123,6 +124,7 @@ def __init__(self, self.preserve_param_in_scan = preserve_param_in_scan self._product_dicts = None self._atom_map = None + self._atom_map_clusters = None self._charge = charge self._multiplicity = multiplicity if reaction_dict is not None: @@ -169,6 +171,36 @@ def atom_map(self, value): """Allow setting the atom map""" self._atom_map = value + @property + def atom_map_clusters(self): + """The distinct reaction channels of this reaction, as clustered atom maps. + + Each entry is an ``arc.mapping.cluster.MapCluster``: one chemically distinct channel, hence one + transition state to search for, whose ``degeneracy`` is that channel's reaction path degeneracy. + More than one cluster is uncommon - it was seen for 1 of 293 mapped reactions of the benchmark + corpus - and marks a reaction for which ``atom_map`` picked one channel by an arbitrary tie-break. + + This is much more expensive than ``atom_map``: it sweeps every RMG template product dictionary, in + both reaction directions, and every superimposable backbone candidate per scissored fragment. It is + therefore computed lazily, on first access only, and is never triggered by ``atom_map``. The result + is cached even when empty, so a failed attempt is not silently repaid on every access; assign + ``None`` to force a recomputation. + + Not persisted by ``as_dict``, unlike ``atom_map`` - it is a derived quantity that can always be + recomputed from the species. + """ + 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 + + @atom_map_clusters.setter + def atom_map_clusters(self, value): + """Allow setting or resetting the atom map clusters""" + self._atom_map_clusters = value + @property def product_dicts(self): """The RMG reaction family product dictionaries""" From 09237cdcd9ece5a47a4fabdab45fc9a7b133353d Mon Sep 17 00:00:00 2001 From: kfir4444 Date: Mon, 24 Aug 2026 10:32:06 +0300 Subject: [PATCH 6/8] Seed the flipped reaction with a family when enumerating atom maps --- arc/mapping/cluster.py | 8 ++++++-- arc/mapping/cluster_test.py | 21 +++++++++++++++++++++ 2 files changed, 27 insertions(+), 2 deletions(-) diff --git a/arc/mapping/cluster.py b/arc/mapping/cluster.py index bc1f60fcf1..34531849ff 100644 --- a/arc/mapping/cluster.py +++ b/arc/mapping/cluster.py @@ -50,7 +50,7 @@ from dataclasses import dataclass, field from arc.common import logger -from arc.mapping.driver import MAX_PDI, map_rxn_all +from arc.mapping.driver import MAX_PDI, map_rxn_all, prepare_flipped_reaction from arc.mapping.engine import flip_map from arc.species import ARCSpecies @@ -547,7 +547,11 @@ def sweep(target, flip: bool) -> None: sweep(rxn, flip=False) if include_flipped and len(maps) < max_maps: - sweep(rxn.flip_reaction(), flip=True) + # Use prepare_flipped_reaction rather than flip_reaction: the latter resets the family, so the + # flipped copy re-derives its product dictionaries with the default family set and comes back + # empty. That is the whole enumeration for a reaction whose template was only discovered in + # reverse, since every forward product dictionary then fails at get_template_product_order. + sweep(prepare_flipped_reaction(rxn), flip=True) if len(maps) >= max_maps: logger.warning(f'enumerate_atom_maps hit the cap of {max_maps} maps for {rxn}; ' diff --git a/arc/mapping/cluster_test.py b/arc/mapping/cluster_test.py index 536e74dc5a..26300d3b75 100644 --- a/arc/mapping/cluster_test.py +++ b/arc/mapping/cluster_test.py @@ -517,6 +517,27 @@ def test_clustered_maps_all_have_a_recipe_predicted_center(self): for center in map_cluster.centers: self.assertIn(center, expected) + def test_enumerates_a_reverse_discovered_reaction(self): + """Test that a reaction whose template was only discovered in reverse still yields maps. + + Every forward product dictionary of such a reaction fails at ``get_template_product_order``, so the + whole enumeration rests on the flipped sweep. That sweep is useless unless the flipped reaction is + seeded with a family, since ``flip_reaction`` resets it and the default family set finds nothing. + """ + rxn = ARCReaction(r_species=[ARCSpecies(label='C2H5Cl', smiles='CCCl')], + p_species=[ARCSpecies(label='C2H4', smiles='C=C'), + ARCSpecies(label='HCl', smiles='Cl')]) + rxn.product_dicts = rxn.get_product_dicts(rmg_family_set='all') + rxn.family = rxn.product_dicts[0]['family'] + rxn.family_own_reverse = rxn.product_dicts[0]['own_reverse'] + self.assertTrue(rxn.product_dicts[0]['discovered_in_reverse']) + atom_maps = cluster.enumerate_atom_maps(rxn, max_maps=20) + self.assertTrue(atom_maps) + for atom_map in atom_maps: + self.assertEqual(sorted(atom_map), list(range(len(atom_map)))) + clusters = cluster.cluster_atom_maps(atom_maps, rxn) + self.assertEqual(len(clusters), 1) + def test_arc_reaction_atom_map_clusters_property(self): """Test that ARCReaction exposes the clusters and caches them.""" clusters = self.rxn.atom_map_clusters From 7fe06b877760ecb14146147d46969775c3849bd1 Mon Sep 17 00:00:00 2001 From: kfir4444 Date: Mon, 24 Aug 2026 12:47:35 +0300 Subject: [PATCH 7/8] Compute reaction path degeneracy exactly from the reaction center orbit --- arc/mapping/cluster.py | 96 +++++++++++++++++++++++++++++++------ arc/mapping/cluster_test.py | 79 ++++++++++++++++++++++++++++-- 2 files changed, 156 insertions(+), 19 deletions(-) diff --git a/arc/mapping/cluster.py b/arc/mapping/cluster.py index 34531849ff..890bc51a16 100644 --- a/arc/mapping/cluster.py +++ b/arc/mapping/cluster.py @@ -43,8 +43,9 @@ which is well defined precisely because hydrogens on a common parent are interchangeable. This keeps hydrogen-transfer reactions - where the migrating atom *is* a hydrogen - fully representable. -See ``docs/atom_mapping_clustering_design.md`` for the full design, and ``docs/atom_mapping_summary.md`` -for a description of the underlying single-map pipeline. +Reaction path degeneracy is not counted from the enumeration, which only ever sees the maps RMG's templates +happened to produce. It is the orbit size of the reaction center under the *full* automorphism group - the +core group extended by permutations of the hydrogens on each core atom - see :func:`center_degeneracy`. """ from dataclasses import dataclass, field @@ -107,29 +108,27 @@ class MapCluster: representative: list[int] members: list[list[int]] = field(default_factory=list) centers: set = field(default_factory=set) + degeneracy: int = 0 key: tuple = () signature: tuple = () truncated: bool = False @property - def degeneracy(self) -> int: + def enumerated_degeneracy(self) -> int: """ - int: The reaction path degeneracy of this channel, i.e. the number of distinct reaction centers. + int: How many distinct reaction centers the enumeration actually found. - This deliberately counts distinct changed-bond sets rather than distinct atom maps. Two maps with - the same changed-bond set differ only by an element of ``Aut(P)`` - they relabel product atoms - without changing which reactant bonds break and form - so they are the same reaction path and must - not be counted twice. For CH4 + OH the four abstractions give four centers, while a map that only - swaps the two resulting water hydrogens adds a member but no new center. - - Note that this is a count of the paths actually *enumerated*. It is a lower bound on the true - degeneracy whenever the enumeration is incomplete. + This counts distinct changed-bond sets rather than distinct atom maps, since two maps sharing a + changed-bond set differ only by an element of ``Aut(P)`` and are the same reaction path. It is a + lower bound on :attr:`degeneracy`, and can be a poor one: RMG generates only two template matches + for ``C2H5Cl -> C2H4 + HCl``, so the enumeration sees at most two of the three equivalent methyl + hydrogens even though the true degeneracy is 3. Compare the two to gauge enumeration coverage. """ return len(self.centers) def __repr__(self) -> str: - return f'' + return f'' def build_complex_graph(species_list: list[ARCSpecies]) -> ComplexGraph: @@ -445,6 +444,69 @@ def _endpoint(graph: ComplexGraph, index: int) -> tuple: return CORE, index +def center_degeneracy(center: frozenset, + graph: ComplexGraph, + automorphisms: list[dict[int, int]], + ) -> int: + """ + The exact reaction path degeneracy of a reaction center: the size of its orbit under the *full* + automorphism group of the reactant complex. + + Counting the distinct centers actually enumerated is a lower bound, and often a poor one - RMG generates + only two template matches for ``C2H5Cl -> C2H4 + HCl``, so enumeration sees at most two of the three + equivalent methyl hydrogens. The orbit size is exact and does not depend on the enumeration at all. + + The full automorphism group is the core-skeleton group extended by arbitrary permutations of the + hydrogens attached to each core atom, and the orbit factorises accordingly:: + + degeneracy = (number of distinct images of the collapsed center under the core group) + * product over core atoms p of n_p! / (n_p - k_p)! + + where ``n_p`` is how many hydrogens hang off core atom ``p`` and ``k_p`` how many distinct ones the + center actually names. The falling factorial counts the ordered choices of which hydrogens play the + named roles. For CH4 + OH the core group is trivial and one of methane's four hydrogens is named, giving + 4; for C2H6 + OH the core group swaps the two carbons and one of three hydrogens is named, giving 6. + + Args: + center (frozenset): The changed-bond set, in reactant indices and *not* hydrogen-collapsed. + graph (ComplexGraph): The reactant complex. + automorphisms (list[dict[int, int]]): The core-skeleton automorphism group. + + Returns: + int: The reaction path degeneracy. + """ + if not center: + return 0 + collapsed = collapse_hydrogens(center, graph) + images = set() + for alpha in automorphisms: + image = list() + for endpoint_a, endpoint_b, order_before, order_after in collapsed: + mapped_a = (endpoint_a[0], alpha.get(endpoint_a[1], endpoint_a[1])) + mapped_b = (endpoint_b[0], alpha.get(endpoint_b[1], endpoint_b[1])) + if mapped_b < mapped_a: + mapped_a, mapped_b = mapped_b, mapped_a + image.append((mapped_a, mapped_b, order_before, order_after)) + images.add(tuple(sorted(image))) + core_orbit = len(images) or 1 + + hydrogen_counts = {v: 0 for v in graph.core} + for hydrogen, parent in graph.parent.items(): + hydrogen_counts[parent] = hydrogen_counts.get(parent, 0) + 1 + named: dict[int, set] = dict() + for i, j, _, _ in center: + for index in (i, j): + if index in graph.parent: + named.setdefault(graph.parent[index], set()).add(index) + + hydrogen_factor = 1 + for parent, hydrogens in named.items(): + available = hydrogen_counts.get(parent, 0) + for offset in range(len(hydrogens)): + hydrogen_factor *= max(available - offset, 1) + return core_orbit * hydrogen_factor + + def canonical_center_key(center: frozenset, automorphisms: list[dict[int, int]]) -> tuple: """ Canonicalize a changed-bond set under the reactant automorphism group: the key is the lexicographic @@ -622,6 +684,7 @@ def cluster_atom_maps(atom_maps: list[list[int]], if key not in clusters: clusters[key] = MapCluster(representative=list(atom_map), key=key, + degeneracy=center_degeneracy(center, r_graph, automorphisms), signature=center_signature(collapsed, orbits, r_graph.symbols), truncated=truncated) clusters[key].members.append(list(atom_map)) @@ -672,7 +735,10 @@ def expected_reaction_centers(rxn, ignore_bond_orders: bool = True) -> set[froze centers, skipped = set(), 0 for product_dict in rxn.product_dicts[:MAX_PDI]: r_label_map = product_dict.get('r_label_map') - if not r_label_map: + if not r_label_map or product_dict.get('discovered_in_reverse'): + # A reverse-discovered dictionary's label maps describe the flipped reaction, so the indices + # they carry are in the flipped frame and cannot be compared against a changed-bond set built + # in this reaction's own frame. Using them anyway predicts centers that match nothing. skipped += 1 continue try: diff --git a/arc/mapping/cluster_test.py b/arc/mapping/cluster_test.py index 26300d3b75..ccdc4eee22 100644 --- a/arc/mapping/cluster_test.py +++ b/arc/mapping/cluster_test.py @@ -363,6 +363,22 @@ def get_expected_changing_bonds(self, r_label_dict, family=None): rxn = _BrokenStub(family='intra_H_migration', product_dicts=[{'r_label_map': {'a': 0}}]) self.assertIsNone(cluster.expected_reaction_centers(rxn)) + def test_reverse_discovered_product_dicts_are_skipped(self): + """Test that a reverse-discovered dictionary is not used to predict centers. + + Its label maps describe the flipped reaction, so the indices they carry are in the flipped frame. + Building a reference from them predicts centers that match nothing, which silently rejects every + enumerated map and drops the reaction back to the relative filter. + """ + rxn = _StubReaction(family='XY_Addition_MultipleBond', + product_dicts=[{'r_label_map': {'*1': 0}, 'discovered_in_reverse': True}], + breaking=[(0, 1)], forming=[(1, 5)]) + self.assertIsNone(cluster.expected_reaction_centers(rxn)) + # The same dictionary discovered forward is usable. + rxn.product_dicts = [{'r_label_map': {'*1': 0}, 'discovered_in_reverse': False}] + self.assertEqual(cluster.expected_reaction_centers(rxn), + {frozenset({(0, 1, 1, 0), (1, 5, 0, 1)})}) + def test_builds_one_center_per_product_dict(self): """Test that each product dictionary contributes its own predicted center.""" rxn = _StubReaction(family='H_Abstraction', @@ -390,21 +406,76 @@ class TestMapCluster(unittest.TestCase): Contains unit tests for the MapCluster container. """ - def test_degeneracy_counts_distinct_centers_not_maps(self): - """Test that degeneracy counts reaction paths, so an Aut(P) relabeling does not inflate it.""" + def test_enumerated_degeneracy_counts_centers_not_maps(self): + """Test that the enumerated count tracks reaction paths, so an Aut(P) relabeling does not inflate it.""" map_cluster = cluster.MapCluster(representative=[0, 1, 2]) map_cluster.members.extend([[0, 1, 2], [0, 2, 1]]) # Both maps describe the same path, so they share a center. map_cluster.centers.add(frozenset({(0, 1, 1, 0)})) self.assertEqual(len(map_cluster.members), 2) - self.assertEqual(map_cluster.degeneracy, 1) + self.assertEqual(map_cluster.enumerated_degeneracy, 1) # A genuinely different path adds a center. map_cluster.centers.add(frozenset({(0, 2, 1, 0)})) - self.assertEqual(map_cluster.degeneracy, 2) + self.assertEqual(map_cluster.enumerated_degeneracy, 2) + + def test_degeneracy_is_stored_not_derived_from_the_enumeration(self): + """Test that the reported degeneracy is the exact orbit size, independent of what was enumerated.""" + map_cluster = cluster.MapCluster(representative=[0, 1, 2], degeneracy=3) + map_cluster.centers.add(frozenset({(0, 1, 1, 0)})) + # Only one path was enumerated, but three exist. + self.assertEqual(map_cluster.enumerated_degeneracy, 1) + self.assertEqual(map_cluster.degeneracy, 3) def test_empty_cluster_has_zero_degeneracy(self): """Test the degenerate case of a cluster with no recorded center.""" self.assertEqual(cluster.MapCluster(representative=[0]).degeneracy, 0) + self.assertEqual(cluster.MapCluster(representative=[0]).enumerated_degeneracy, 0) + + +class TestCenterDegeneracy(unittest.TestCase): + """ + Contains unit tests for the exact, orbit-based reaction path degeneracy. + """ + + @staticmethod + def _degeneracy(smiles_list, center): + """Compute the exact degeneracy of a center on a complex built from a list of SMILES.""" + graph = cluster.build_complex_graph([ARCSpecies(label=f's{i}', smiles=smiles) + for i, smiles in enumerate(smiles_list)]) + automorphisms, _ = cluster.core_automorphisms(graph) + return cluster.center_degeneracy(center, graph, automorphisms), graph + + def test_methane_abstraction(self): + """Test that abstracting one of methane's four equivalent hydrogens has degeneracy 4.""" + # CH4 + OH: 0=C, 1-4=H, 5=O, 6=H. Break C-H(1), form H(1)-O. + degeneracy, _ = self._degeneracy(['C', '[OH]'], frozenset({(0, 1, 1, 0), (1, 5, 0, 1)})) + self.assertEqual(degeneracy, 4) + + def test_ethane_abstraction_combines_core_symmetry_and_hydrogen_choice(self): + """Test that ethane's degeneracy of 6 is the 2 equivalent carbons times 3 hydrogens each.""" + # CC: 0,1=C, 2-7=H. Break C(0)-H(2). + graph = cluster.build_complex_graph([ARCSpecies(label='ethane', smiles='CC')]) + automorphisms, _ = cluster.core_automorphisms(graph) + hydrogen = min(h for h, parent in graph.parent.items() if parent == graph.core[0]) + center = frozenset({(graph.core[0], hydrogen, 1, 0)}) + self.assertEqual(cluster.center_degeneracy(center, graph, automorphisms), 6) + + def test_propane_primary_and_secondary_differ(self): + """Test that propane's primary sites give 6 and its secondary site 2.""" + graph = cluster.build_complex_graph([ARCSpecies(label='propane', smiles='CCC')]) + automorphisms, _ = cluster.core_automorphisms(graph) + hydrogens = dict() + for hydrogen, parent in graph.parent.items(): + hydrogens.setdefault(parent, list()).append(hydrogen) + primary = frozenset({(0, min(hydrogens[0]), 1, 0)}) + secondary = frozenset({(1, min(hydrogens[1]), 1, 0)}) + self.assertEqual(cluster.center_degeneracy(primary, graph, automorphisms), 6) + self.assertEqual(cluster.center_degeneracy(secondary, graph, automorphisms), 2) + + def test_empty_center(self): + """Test that a center with no changed bonds has zero degeneracy.""" + degeneracy, _ = self._degeneracy(['C'], frozenset()) + self.assertEqual(degeneracy, 0) class TestClusteringIntegration(unittest.TestCase): From c8b75f7fcec42ffab3fe5160870eefb4b3d22b43 Mon Sep 17 00:00:00 2001 From: kfir4444 Date: Mon, 24 Aug 2026 13:26:50 +0300 Subject: [PATCH 8/8] Test that sec-butyl to n-butyl resolves into its two migration channels --- arc/mapping/cluster_test.py | 41 +++++++++++++++++++++++++++++++++++++ 1 file changed, 41 insertions(+) diff --git a/arc/mapping/cluster_test.py b/arc/mapping/cluster_test.py index ccdc4eee22..a673c65c1a 100644 --- a/arc/mapping/cluster_test.py +++ b/arc/mapping/cluster_test.py @@ -609,6 +609,47 @@ def test_enumerates_a_reverse_discovered_reaction(self): clusters = cluster.cluster_atom_maps(atom_maps, rxn) self.assertEqual(len(clusters), 1) + def test_sec_butyl_to_n_butyl_has_two_channels(self): + """Test that sec-butyl -> n-butyl is recognized as two distinct channels, three paths each. + + ``CC[CH]C <=> CCC[CH2]`` is C0-C1-C2(.)-C3 losing a hydrogen from a terminal methyl onto the radical + carbon C2. Either methyl works, and both give n-butyl, because n-butyl is the same molecule numbered + from either end: + + C3 -> C2 is a 1,2 shift through a three-membered TS + C0 -> C2 is a 1,3 shift through a four-membered TS + + The two are not symmetry-equivalent - sec-butyl has no automorphism exchanging its ends, and the two + transition states are structurally different - so they must not merge into one channel of degeneracy + 6. Each carries the three hydrogens of its methyl. + + Both transition states have been located and optimized separately, with distinct barriers, by + supplying the two atom maps below to ARC by hand. This test pins that the clustering derives the same + two channels without them. + """ + rxn = ARCReaction(r_species=[ARCSpecies(label='sec-butyl', smiles='CC[CH]C')], + p_species=[ARCSpecies(label='n-butyl', smiles='CCC[CH2]')]) + self.assertEqual(rxn.family, 'intra_H_migration') + + # The two hand-written maps, one per located transition state. + one_two_shift = [0, 1, 2, 3, 6, 4, 5, 8, 7, 10, 9, 12, 11] + one_three_shift = [3, 2, 1, 0, 7, 12, 11, 10, 9, 8, 4, 5, 6] + manual = cluster.cluster_atom_maps([one_two_shift, one_three_shift], rxn) + self.assertEqual(len(manual), 2) + for map_cluster in manual: + self.assertEqual(map_cluster.degeneracy, 3) + + # The enumeration finds the same two channels unaided. + enumerated = cluster.cluster_atom_maps(cluster.enumerate_atom_maps(rxn), rxn) + self.assertEqual(len(enumerated), 2) + self.assertEqual(sorted(c.degeneracy for c in enumerated), [3, 3]) + self.assertEqual({c.key for c in manual}, {c.key for c in enumerated}) + + # One channel moves a hydrogen off C0, the other off C3; both onto the radical carbon C2. + donors = {tuple(sorted(entry[:2])) for c in enumerated for entry in c.signature + if (entry[2], entry[3]) == (1, 0)} + self.assertEqual(len(donors), 2) + def test_arc_reaction_atom_map_clusters_property(self): """Test that ARCReaction exposes the clusters and caches them.""" clusters = self.rxn.atom_map_clusters