Support reactions where a species participates more than once (A + A)

arc/checks/nmd.py

@@ -34,6 +34,9 @@ def analyze_ts_normal_mode_displacement(reaction: ARCReaction,
     Analyze the normal mode displacement by identifying bonds that break and form
     and comparing them to the expected given reaction.
     Note that the TS geometry must be in the standard orientation for the normal mode displacement to be relevant.
+    The forming, breaking and changed bonds of the reaction are indexed into the concatenated reactant atoms,
+    where a species participating more than once (e.g. ``A + A``, for which ``r_species`` holds a single entry)
+    contributes its atoms once per occurrence. The TS geometry must span that same set of atoms.
 
     Args:
         reaction (ARCReaction): The reaction for which the TS is checked.
@@ -46,9 +49,23 @@ def analyze_ts_normal_mode_displacement(reaction: ARCReaction,
 
     Returns:
         bool | None: Whether the TS normal mode displacement is consistent with the desired reaction.
+                     ``None`` if the analysis could not be performed, either because no frequency job was
+                     given, because the TS geometry does not span the reactant atoms the bond indices refer
+                     to, or because no normal mode displacements could be parsed from the job's output file.
+                     Only ``False`` marks the TS as inconsistent with the reaction.
     """
     if job is None:
         return None
+    ts_xyz = reaction.ts_species.get_xyz()
+    n_ts = len(ts_xyz['symbols'])
+    reactants, _ = reaction.get_reactants_and_products(return_copies=False)
+    n_expected = sum(spc.number_of_atoms for spc in reactants)
+    if n_ts != n_expected:
+        logger.warning(f'The geometry of TS {reaction.ts_species.label} has {n_ts} atoms, while the reactants of '
+                       f'reaction {reaction.label} have {n_expected} atoms in total. The forming, breaking and '
+                       f'changed bonds of the reaction are indexed into the reactant atoms, so they cannot be '
+                       f'applied to this TS geometry. Skipping the normal mode displacement analysis.')
+        return None
     if reaction.atom_map is None:
         # Without an atom map the formed/broken/changed bonds cannot be determined; skip the NMD
         # check for this reaction rather than crashing the whole run (e.g. reactions ARC cannot map).
@@ -58,11 +75,6 @@ def analyze_ts_normal_mode_displacement(reaction: ARCReaction,
                 not in reaction.ts_species.ts_checks['warnings']:
             reaction.ts_species.ts_checks['warnings'] += 'Atom map is None; skipped the TS normal mode displacement check; '
         return None
-    ts_xyz = reaction.ts_species.get_xyz()
-    n_ts = len(ts_xyz['symbols'])
-    n_expected = sum(spc.number_of_atoms for spc in reaction.r_species)
-    if n_ts != n_expected:
-        return False
     try:
         freqs, normal_mode_disp = parser.parse_normal_mode_displacement(log_file_path=job.local_path_to_output_file)
     except NotImplementedError:
@@ -486,13 +498,46 @@ def get_bond_length_in_reaction(bond: tuple[int, int] | list[int],
     return float(distance)
 
 
+def get_repeated_species_atom_equivalences(reaction: ARCReaction,
+                                           well: int = 0,
+                                           ) -> list[list[int]]:
+    """
+    Build atom-equivalence groups across repeated identical species in a reaction well.
+
+    When a species participates more than once (e.g. OH + OH), ``r_species`` is deduplicated, so the
+    atom map may assign a reactive atom to one copy while the located TS uses the equivalent atom of
+    another copy. For each atom position within such a species, the atoms at that position across all
+    copies are equivalent. Indices follow the expanded (per-occurrence) atom ordering used by the atom
+    map and by ``get_reactants_and_products``.
+
+    Args:
+        reaction (ARCReaction): The reaction.
+        well (int): ``0`` for the reactants well, ``1`` for the products well.
+
+    Returns:
+        list[list[int]]: Equivalence groups (each a list of global atom indices) for repeated species.
+    """
+    species = reaction.r_species if well == 0 else reaction.p_species
+    groups, offset = list(), 0
+    for spc in species:
+        count = reaction.get_species_count(species=spc, well=well)
+        n_atoms = spc.number_of_atoms
+        if count > 1:
+            for pos in range(n_atoms):
+                groups.append([offset + copy_i * n_atoms + pos for copy_i in range(count)])
+        offset += count * n_atoms
+    return groups
+
+
 def find_equivalent_atoms(reaction: ARCReaction,
                           reactant_only: bool = True,
                           ) -> tuple[list[list[int]], list[list[int]]]:
     """
     Find equivalent atoms in the reactants and products of a reaction.
     This is a tentative function that should be replaced when atom mapping returns a list.
     It is meant to suggest additional atoms that can move instead of the ones selected by the atom map.
+    Reactant equivalences cover both atoms made equivalent within a molecule and atoms at the same
+    position across the copies of a species that participates more than once.
 
     Args:
         reaction (ARCReaction): The reaction for which equivalent atoms are searched.
@@ -510,6 +555,7 @@ def find_equivalent_atoms(reaction: ARCReaction,
                                                                 inc=sum([len(r.mol.atoms) for r in reactants[:i]]),
                                                                 atom_map=None,
                                                                 ))
+    r_eq_atoms.extend(get_repeated_species_atom_equivalences(reaction, well=0))
     if not reactant_only:
         for i, product in enumerate(products):
             p_eq_atoms.extend(identify_equivalent_atoms_in_molecule(molecule=product.mol,

arc/checks/nmd_test.py

@@ -8,6 +8,7 @@
 import unittest
 import os
 import shutil
+from unittest.mock import patch
 
 import numpy as np
 
@@ -56,19 +57,19 @@ def setUpClass(cls):
                          H      -0.34927558    0.98159583   -0.32768232
                          H      -0.02233792   -0.04887375    1.09087665
                          H       1.02216551   -0.15844188   -0.35067554"""
-        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.oh_xyz = """O       0.48890387    0.00000000    0.00000000
+                        H      -0.48890387    0.00000000    0.00000000"""
+        cls.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"""
+        cls.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_1 = ARCReaction(r_species=[ARCSpecies(label='CH4', smiles='C', xyz=cls.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)])
+                                           ARCSpecies(label='OH', smiles='[OH]', xyz=cls.oh_xyz)],
+                                p_species=[ARCSpecies(label='CH3', smiles='[CH3]', xyz=cls.ch3_xyz),
+                                           ARCSpecies(label='H2O', smiles='O', xyz=cls.h2o_xyz)])
         cls.ts_1_xyz = check_xyz_dict("""C                    -1.212192   -0.010161    0.000000
                                          H                     0.010122    0.150115    0.000001
                                          H                    -1.460491   -0.555461   -0.907884
@@ -225,6 +226,53 @@ def setUpClass(cls):
                                          H       -1.951439    0.465285   -1.158262""")
         cls.rxn_3.ts_species = ARCSpecies(label='TS3', is_ts=True, xyz=cls.ts_3_xyz)
 
+        cls.ts_oh_oh_xyz = check_xyz_dict("""O      -1.10000000    0.00000000    0.00000000
+                                             H      -0.15000000    0.00000000    0.00000000
+                                             O       1.20000000    0.00000000    0.00000000
+                                             H       1.55000000    0.85000000    0.00000000""")
+
+    def test_get_repeated_species_atom_equivalences(self):
+        """Repeated identical reactants (e.g. OH + OH) must yield cross-copy atom-equivalence groups
+        (the two O's, the two H's) so NMD can try the alternative mapping when the atom map picks one
+        copy's atom but the TS uses the other's. Non-repeated reactions must yield none."""
+        rxn = ARCReaction(r_species=[ARCSpecies(label='OH', smiles='[OH]'),
+                                     ARCSpecies(label='OH', smiles='[OH]')],
+                          p_species=[ARCSpecies(label='H2O', smiles='O'),
+                                     ARCSpecies(label='O', smiles='[O]', multiplicity=3)],
+                          reactants=['OH', 'OH'], products=['H2O', 'O'])
+        groups = nmd.get_repeated_species_atom_equivalences(rxn, well=0)
+        self.assertEqual(sorted(sorted(g) for g in groups), [[0, 2], [1, 3]])
+        self.assertEqual(nmd.get_repeated_species_atom_equivalences(self.rxn_1, well=0), [])
+
+    def test_analyze_ts_normal_mode_displacement_repeated_reactant(self):
+        """Test that a TS of an A + A reaction is not reported as inconsistent with the reaction."""
+        rxn = ARCReaction(r_species=[ARCSpecies(label='OH', smiles='[OH]', xyz=self.oh_xyz),
+                                     ARCSpecies(label='OH', smiles='[OH]', xyz=self.oh_xyz)],
+                          p_species=[ARCSpecies(label='H2O', smiles='O', xyz=self.h2o_xyz),
+                                     ARCSpecies(label='O', smiles='[O]', multiplicity=3,
+                                                xyz='O 0.00000000 0.00000000 0.00000000')])
+        rxn.ts_species = ARCSpecies(label='TS_OH_OH', is_ts=True, xyz=self.ts_oh_oh_xyz)
+        self.assertEqual(len(rxn.r_species), 1)
+        reactants, _ = rxn.get_reactants_and_products(return_copies=False)
+        self.assertEqual(sum(spc.number_of_atoms for spc in reactants), 4)
+        self.assertEqual(len(rxn.ts_species.get_xyz()['symbols']), 4)
+        self.generic_job.local_path_to_output_file = os.path.join(ARC_TESTING_PATH, 'freq', 'CHO_neg_freq.out')
+        with patch.object(nmd.parser, 'parse_normal_mode_displacement', side_effect=NotImplementedError):
+            valid = nmd.analyze_ts_normal_mode_displacement(reaction=rxn, job=self.generic_job, amplitude=0.25)
+        self.assertIsNone(valid)
+
+    def test_analyze_ts_normal_mode_displacement_atom_count_mismatch(self):
+        """Test that a TS spanning a different set of atoms than the reactants is not reported as consistent."""
+        rxn = ARCReaction(r_species=[ARCSpecies(label='CH4', smiles='C', xyz=self.ch4_xyz),
+                                     ARCSpecies(label='OH', smiles='[OH]', xyz=self.oh_xyz)],
+                          p_species=[ARCSpecies(label='CH3', smiles='[CH3]', xyz=self.ch3_xyz),
+                                     ARCSpecies(label='H2O', smiles='O', xyz=self.h2o_xyz)])
+        rxn.ts_species = ARCSpecies(label='TS_short', is_ts=True, xyz=self.ch4_xyz)
+        self.generic_job.local_path_to_output_file = os.path.join(ARC_TESTING_PATH, 'freq', 'TS_CH4_OH.log')
+        valid = nmd.analyze_ts_normal_mode_displacement(reaction=rxn, job=self.generic_job, amplitude=0.25)
+        self.assertIsNot(valid, True)
+        self.assertIsNone(valid)
+
     def test_analyze_ts_normal_mode_displacement_simple_rxns(self):
         """Test the analyze_ts_normal_mode_displacement() function with simple reactions."""
         # CH4 + OH <=> CH3 + H2O

arc/job/adapters/ts/autotst_ts.py

@@ -211,6 +211,8 @@ def set_input_file_memory(self) -> None:
     def execute_incore(self):
         """
         Execute a job incore.
+        The reverse-direction reaction label repeats each species by the number of times it
+        participates in its well, so a well with a repeated species stays atom-balanced.
         """
         if not AUTOTST_PYTHON or not os.path.isfile(AUTOTST_PYTHON):
             raise FileNotFoundError('AutoTST python executable was not found. '
@@ -230,10 +232,14 @@ def execute_incore(self):
                                                 multiplicity=rxn.multiplicity,
                                                 )
                 reaction_label_fwd = get_autotst_reaction_string(rxn)
+                rev_reactant_labels = [lbl for lbl in rxn.products
+                                       for _ in range(rxn.get_species_count(label=lbl, well=1))]
+                rev_product_labels = [lbl for lbl in rxn.reactants
+                                      for _ in range(rxn.get_species_count(label=lbl, well=0))]
                 reaction_label_rev = get_autotst_reaction_string(ARCReaction(r_species=rxn.p_species,
                                                                              p_species=rxn.r_species,
-                                                                             reactants=rxn.products,
-                                                                             products=rxn.reactants))
+                                                                             reactants=rev_reactant_labels,
+                                                                             products=rev_product_labels))
 
                 i = 0
                 for reaction_label, direction in zip([reaction_label_fwd, reaction_label_rev], ['F', 'R']):

arc/processor.py

@@ -32,6 +32,30 @@ def resolve_neb_level(ts_adapters: list) -> Level | None:
     return None
 
 
+def _thermo_lib_has_isomorph(spc: 'ARCSpecies', species_list: list) -> bool:
+    """
+    Whether ``species_list`` already contains a species with the same molecular identity as ``spc``.
+
+    Arkane keys thermo-library entries by adjacency list + multiplicity and rejects duplicates, so a
+    reaction with identical participants (e.g. OH + OH) must contribute each unique species only once.
+
+    Args:
+        spc (ARCSpecies): The candidate species.
+        species_list (list): Species already collected for the thermo library.
+
+    Returns:
+        bool: ``True`` if an isomorphic, same-multiplicity species is already present.
+    """
+    if spc.mol is None:
+        return False
+    for existing in species_list:
+        if existing.mol is not None \
+                and existing.multiplicity == spc.multiplicity \
+                and existing.mol.is_isomorphic(spc.mol):
+            return True
+    return False
+
+
 def process_arc_project(thermo_adapter: str,
                         kinetics_adapter: str,
                         project: str,
@@ -196,7 +220,7 @@ def process_arc_project(thermo_adapter: str,
                                             )
         statmech_adapter.compute_thermo(e0_only=True)
     for spc in converged_species:
-        if spc.thermo is not None:
+        if spc.thermo is not None and not _thermo_lib_has_isomorph(spc, species_for_thermo_lib):
             species_for_thermo_lib.append(spc)
         plotter.augment_arkane_yml_file_with_mol_repr(spc, output_directory)
     if species_for_thermo_lib:

arc/reaction/reaction.py

@@ -873,6 +873,8 @@ def get_single_mapped_product_xyz(self) -> ARCSpecies | None:
     def get_reactants_xyz(self, return_format='str') -> dict | str:
         """
         Get a combined string/dict representation of the cartesian coordinates of all reactant species.
+        A species participating more than once in the reactant well contributes its atoms once per
+        occurrence.
 
         Args:
             return_format (str): Either ``'dict'`` to return a dict format or ``'str'`` to return a string format.
@@ -885,15 +887,17 @@ def get_reactants_xyz(self, return_format='str') -> dict | str:
             Orient the fragments according to the reactive site.
         """
         xyz_dict = dict()
-        if len(self.r_species) == 1:
-            xyz_dict = self.r_species[0].get_xyz()
-        elif len(self.r_species) >= 2:
+        reactants = [spc for spc in self.r_species
+                     for _ in range(self.get_species_count(species=spc, well=0))]
+        if len(reactants) == 1:
+            xyz_dict = reactants[0].get_xyz()
+        elif len(reactants) >= 2:
             xyz_dict = {'symbols': tuple(), 'isotopes': tuple(), 'coords': tuple()}
-            for i, reactant in enumerate(self.r_species):
+            for i, reactant in enumerate(reactants):
                 xyz = translate_to_center_of_mass(reactant.get_xyz())
                 if i:
                     xyz = translate_xyz(xyz_dict=xyz,
-                                        translation=(sum(spc.radius for spc in self.r_species[:i]) * 1.1 * i, 0, 0))
+                                        translation=(sum(spc.radius for spc in reactants[:i]) * 1.1 * i, 0, 0))
                 xyz_dict['symbols'] += xyz['symbols']
                 xyz_dict['isotopes'] += xyz['isotopes']
                 xyz_dict['coords'] += xyz['coords']
@@ -906,6 +910,8 @@ def get_products_xyz(self, return_format='str') -> dict | str:
         """
         Get a combined string/dict representation of the cartesian coordinates of all product species.
         The resulting coordinates are ordered as the reactants using an atom map.
+        A species participating more than once in the product well contributes its atoms once per
+        occurrence.
 
         Args:
             return_format (str): Either ``'dict'`` to return a dict format or ``'str'`` to return a string format.
@@ -917,15 +923,17 @@ def get_products_xyz(self, return_format='str') -> dict | str:
         Todo:
             Orient the fragments according to the reactive site.
         """
-        if len(self.p_species) == 1:
-            xyz_dict = self.p_species[0].get_xyz()
+        products = [spc for spc in self.p_species
+                    for _ in range(self.get_species_count(species=spc, well=1))]
+        if len(products) == 1:
+            xyz_dict = products[0].get_xyz()
         else:
             xyz_dict = {'symbols': tuple(), 'isotopes': tuple(), 'coords': tuple()}
-            for i, product in enumerate(self.p_species):
+            for i, product in enumerate(products):
                 xyz = translate_to_center_of_mass(product.get_xyz())
                 if i:
                     xyz = translate_xyz(xyz_dict=xyz,
-                                        translation=(sum(spc.radius for spc in self.p_species[:i]) * 1.1 * i, 0, 0))
+                                        translation=(sum(spc.radius for spc in products[:i]) * 1.1 * i, 0, 0))
                 xyz_dict['symbols'] += xyz['symbols']
                 xyz_dict['isotopes'] += xyz['isotopes']
                 xyz_dict['coords'] += xyz['coords']