Several more atom mapping improvements

arc/job/adapters/ts/heuristics_test.py

@@ -36,6 +36,10 @@
 from arc.species.species import ARCSpecies
 from arc.species.zmat import _compare_zmats, get_parameter_from_atom_indices
 
+from arc.species.species import check_isomorphism
+from arc.species.zmat import remove_zmat_atom_0
+from arc.species.converter import relocate_zmat_dummy_atoms_to_the_end
+
 
 class TestHeuristicsAdapter(unittest.TestCase):
     """
@@ -1409,11 +1413,32 @@ def test_get_new_zmat2_map(self):
         # expected_new_map = {0: 12, 1: 13, 2: 'X24', 3: 14, 4: 15, 5: 16, 6: 'X25', 7: 17, 8: 'X26', 9: 18, 10: 19,
         #                     11: 20, 12: 21, 13: 22, 14: 'X27', 15: 23, 16: 'X28', 17: 2, 18: 3, 19: 1, 21: 4, 23: 0,
         #                     25: 7, 26: 6, 28: 5, 20: 'X8', 22: 'X9', 24: 'X10', 27: 'X11'}
-        expected_new_map = {0: 12, 1: 13, 2: 'X24', 3: 14, 4: 15, 5: 16, 6: 'X25', 7: 17, 8: 'X26', 9: 18, 10: 19,
-                            11: 20, 12: 21, 13: 22, 14: 'X27', 15: 23, 16: 'X28', 17: 2, 18: 1, 19: 3, 21: 0, 23: 4,
-                            25: 5, 26: 6, 28: 7, 20: 'X8', 22: 'X9', 24: 'X10', 27: 'X11'}
-
-        self.assertEqual(new_map, expected_new_map)
+
+        # Test isomorphism of the mapped reactant_2 part
+        zmat_2_mod = remove_zmat_atom_0(self.zmat_6)
+        zmat_2_mod['map'] = relocate_zmat_dummy_atoms_to_the_end(zmat_2_mod['map'])
+        spc_from_zmat_2 = ARCSpecies(label='spc_from_zmat_2', xyz=zmat_2_mod, multiplicity=reactant_2.multiplicity,
+                                     number_of_radicals=reactant_2.number_of_radicals, charge=reactant_2.charge)
+
+        # Verify that all physical atom indices in new_map that came from zmat_2 correctly map to reactant_2
+        # Atom indices in new_map are for the combined species.
+        # Atoms 0-16 in self.zmat_5, atoms 1-12 in self.zmat_6 (13 atoms total, index 0 removed).
+        # In get_new_zmat_2_map, zmat_2 atoms are mapped to indices in new_map.
+
+        num_atoms_1 = len(self.zmat_5['symbols'])
+        atom_map = dict()
+        for i in range(1, len(self.zmat_6['symbols'])):
+            if not isinstance(self.zmat_6['symbols'][i], str) or self.zmat_6['symbols'][i] != 'X':
+                # This is a physical atom in zmat_2 (at index i)
+                # Its index in the combined Z-Matrix is num_atoms_1 + i - 1
+                combined_idx = num_atoms_1 + i - 1
+                if combined_idx in new_map:
+                    # new_map[combined_idx] is the index in reactant_2
+                    # i-1 is the index in spc_from_zmat_2
+                    atom_map[i-1] = new_map[combined_idx]
+
+        # Verify the atom_map is a valid isomorphism
+        self.assertTrue(check_isomorphism(spc_from_zmat_2.mol, reactant_2.mol, atom_map))
 
     def test_get_new_map_based_on_zmat_1(self):
         """Test the get_new_map_based_on_zmat_1() function."""

arc/mapping/driver.py

@@ -272,6 +272,13 @@ def map_rxn(rxn: ARCReaction,
     r_bdes, p_bdes = find_all_breaking_bonds(rxn, r_direction=True, pdi=pdi), find_all_breaking_bonds(rxn, r_direction=False, pdi=pdi)
     r_cuts, p_cuts = cut_species_based_on_atom_indices(reactants, r_bdes), cut_species_based_on_atom_indices(products, p_bdes)
 
+    if r_cuts is None or p_cuts is None:
+        if rxn.product_dicts is not None and len(rxn.product_dicts) - 1 > pdi < MAX_PDI:
+            return map_rxn(rxn, backend=backend, product_dict_index_to_try=pdi + 1)
+        else:
+            logger.error(f'Could not cut species for reaction {rxn}')
+            return None
+
     try:
         r_label_map = rxn.product_dicts[pdi]['r_label_map']
         p_label_map = rxn.product_dicts[pdi]['p_label_map']

arc/mapping/driver_test.py

@@ -1120,7 +1120,7 @@ def test_get_atom_map_2(self):
                                                                     11 H u0 p0 c0 {4,S}
                                                                     12 H u0 p0 c0 {5,S}""")
         rxn = ARCReaction(reactants=['C6H6_a'], products=['C6H6_b'], r_species=[r_1], p_species=[p_1])
-        self.assertEqual(rxn.atom_map, [3, 2, 1, 0, 5, 4, 10, 9, 8, 7, 6, 11])
+        self.assertIn(rxn.atom_map, [[3, 2, 1, 0, 5, 4, 10, 9, 8, 7, 6, 11], [3, 2, 1, 0, 5, 4, 10, 9, 8, 6, 7, 11]])
         self.assertTrue(check_atom_map(rxn))
 
         # Disproportionation: HO2 + NHOH <=> NH2OH + O2
@@ -1373,7 +1373,7 @@ def test_get_atom_map_6(self):
                                                                     11 H u0 p0 c0 {4,S}
                                                                     12 H u0 p0 c0 {5,S}""")
         rxn = ARCReaction(reactants=['C6H6_1'], products=['C6H6_b'], r_species=[r_1], p_species=[p_1])
-        self.assertEqual(rxn.atom_map, [3, 2, 1, 0, 5, 4, 10, 9, 8, 7, 6, 11])
+        self.assertIn(rxn.atom_map, [[3, 2, 1, 0, 5, 4, 10, 9, 8, 7, 6, 11], [3, 2, 1, 0, 5, 4, 10, 9, 8, 6, 7, 11]])
         self.assertTrue(check_atom_map(rxn))
 
     def test_get_atom_map_7(self):
@@ -1420,10 +1420,12 @@ def test_get_atom_map_8(self):
         rxn = ARCReaction(r_species=[ARCSpecies(label="r1", smiles="F[C]F", xyz=r1_xyz),
                                      ARCSpecies(label="r2", smiles="[CH3]", xyz=r2_xyz)],
                           p_species=[ARCSpecies(label="p1", smiles="F[C](F)C", xyz=p1_xyz)])
-        self.assertIn(rxn.atom_map[:2], [[0, 1], [1, 0]])
-        self.assertEqual(rxn.atom_map[2], 2)
-        self.assertEqual(rxn.atom_map[3], 3)
-        self.assertIn(tuple(rxn.atom_map[4:]), list(permutations([4, 5, 6])))
+        if rxn.atom_map[0] == 0:
+            self.assertEqual(rxn.atom_map[:4], [0, 1, 2, 3])
+        else: # Only other F can be in position 0.
+            self.assertEqual(rxn.atom_map[:4], [2, 1, 0, 3])
+        self.assertIn(tuple(rxn.atom_map[4:]), tuple(permutations([4, 5, 6])))
+
         self.assertTrue(check_atom_map(rxn))
 
     def test_get_atom_map_9(self):
@@ -1542,22 +1544,42 @@ def test_get_atom_map_11(self):
         rxn = ARCReaction(reactants=['C4H10', 'CO'], products=['C5H10O'],
                           r_species=[r_1, r_2], p_species=[p_1])
         atom_map = rxn.atom_map
-        self.assertEqual(atom_map[:4], [0, 1, 2, 3])
-        self.assertIn(tuple(rxn.atom_map[4:7]), permutations([6, 7, 8]))
-        self.assertEqual(atom_map[7], 15)
-        self.assertIn(tuple(rxn.atom_map[8:11]), permutations([9, 10, 11]))
-        self.assertIn(tuple(rxn.atom_map[11:14]), permutations([12, 13, 14]))
-        self.assertEqual(atom_map[14:], [4, 5])
         self.assertTrue(check_atom_map(rxn))
+        # Set all anchor- atoms uneffected by symmetry.
+        self.assertEqual(atom_map[1], 1) # Middle Carbon
+        self.assertEqual(atom_map[7], 15) # Middel Hydrogen
+        self.assertEqual(atom_map[-2:], [4, 5]) # CO (In that order!)
+        # Check the symmetric carbons:
+        symm_carbon_hydrogens_r = {
+                0: [4, 5, 6],
+                2: [8, 9, 10],
+                3: [11, 12, 13]
+        }
+        symm_carbon_hydrogens_p = {
+                0: [6, 7, 8],
+                2: [9, 10, 11],
+                3: [12, 13, 14]
+        }
+        for r_atom, p_atom in enumerate(atom_map[:4]):
+            if r_atom == 1:
+                continue # anchor carbon.
+            self.assertIn(p_atom, [0, 2, 3])
+            for r_h in symm_carbon_hydrogens_r[r_atom]:
+                self.assertIn(atom_map[r_h], symm_carbon_hydrogens_p[p_atom])
+
+
         # same reaction in reverse:
         rxn_rev = ARCReaction(r_species=[p_1], p_species=[r_1, r_2])
         atom_map = rxn_rev.atom_map
-        for index in [0, 2, 3]:
-            self.assertIn(atom_map[index], [0, 2, 3])
-        self.assertEqual(atom_map[1], 1)
-        self.assertEqual(atom_map[4], 14)
-        self.assertEqual(atom_map[5], 15)
-        self.assertEqual(atom_map[15], 7)
+        self.assertEqual(atom_map[1], 1) # Middle Carbon
+        self.assertEqual(atom_map[15], 7) # Middel Hydrogen
+        self.assertEqual(atom_map[4:6], [14, 15]) # CO (In that order!)
+        for r_atom, p_atom in enumerate(atom_map[:4]):
+            if r_atom == 1:
+                continue # anchor carbon.
+            self.assertIn(p_atom, [0, 2, 3])
+            for p_h in symm_carbon_hydrogens_p[r_atom]:
+                self.assertIn(atom_map[p_h], symm_carbon_hydrogens_r[p_atom])
         self.assertTrue(check_atom_map(rxn_rev))
 
     def test_get_atom_map_12(self):

arc/mapping/engine.py

@@ -285,12 +285,14 @@ def identify_superimposable_candidates(fingerprint_1: dict[int, dict[str, str |
                               of species 1, values are potentially mapped atom indices of species 2.
     """
     candidates = list()
-    for key_1 in fingerprint_1.keys():
-        for key_2 in fingerprint_2.keys():
-            # Try all combinations of heavy atoms.
-            result = iterative_dfs(fingerprint_1, fingerprint_2, key_1, key_2)
-            if result is not None:
-                candidates.append(result)
+    if not fingerprint_1:
+        return []
+    key_1 = list(fingerprint_1.keys())[0]
+    for key_2 in fingerprint_2.keys():
+        # Try all combinations of heavy atoms.
+        result = iterative_dfs(fingerprint_1, fingerprint_2, key_1, key_2)
+        if result is not None:
+            candidates.append(result)
     return prune_identical_dicts(candidates)
 
 
@@ -382,14 +384,7 @@ def prune_identical_dicts(dicts_list: list[dict]) -> list[dict]:
     """
     new_dicts_list = list()
     for new_dict in dicts_list:
-        unique_ = True
-        for existing_dict in new_dicts_list:
-            if unique_:
-                for new_key, new_val in new_dict.items():
-                    if new_key not in existing_dict.keys() or new_val == existing_dict[new_key]:
-                        unique_ = False
-                        break
-        if unique_:
+        if new_dict not in new_dicts_list:
             new_dicts_list.append(new_dict)
     return new_dicts_list
 
@@ -1195,11 +1190,18 @@ def pairing_reactants_and_products_for_mapping(r_cuts: list[ARCSpecies],
         list[tuple[ARCSpecies,ARCSpecies]]: A list of paired reactant and products, to be sent to map_two_species.
     """
     pairs: list[tuple[ARCSpecies, ARCSpecies]] = list()
-    for react in r_cuts:
+    r_res = [generate_resonance_structures_safely(react.mol, save_order=True) or [react.mol] for react in r_cuts]
+    for i, react in enumerate(r_cuts):
+        res1 = r_res[i]
         for idx, prod in enumerate(p_cuts):
-            if r_cut_p_cut_isomorphic(react, prod):
-                pairs.append((react, prod))
-                p_cuts.pop(idx)
+            found = False
+            for res in res1:
+                if res.fingerprint == prod.mol.fingerprint or prod.mol.is_isomorphic(res, save_order=True):
+                    pairs.append((react, prod))
+                    p_cuts.pop(idx)
+                    found = True
+                    break
+            if found:
                 break
     return pairs
 
@@ -1460,7 +1462,10 @@ def copy_species_list_for_mapping(species: list["ARCSpecies"]) -> list["ARCSpeci
     Returns:
         list[ARCSpecies]: The copied species list.
     """
-    copies = [spc.copy() for spc in species]
+    copies = list()
+    for spc in species:
+        new_spc = ARCSpecies(label=spc.label, mol=spc.mol.copy(deep=True), xyz=spc.get_xyz(), keep_mol=True)
+        copies.append(new_spc)
     for copy, spc in zip(copies, species):
         for atom1, atom2 in zip(copy.mol.atoms, spc.mol.atoms):
             atom1.label = atom2.label

arc/mapping/engine_test.py

@@ -826,7 +826,7 @@ def test_identify_superimposable_candidates(self):
 
         candidates = engine.identify_superimposable_candidates(fingerprint_1=self.butenylnebzene_fingerprint,
                                                         fingerprint_2=self.butenylnebzene_fingerprint)
-        self.assertEqual(candidates, [{0: 0, 5: 5, 4: 4, 3: 3, 2: 2, 1: 1, 6: 6, 7: 7, 8: 8, 9: 9}])
+        self.assertEqual(candidates[0], {0: 0, 5: 5, 4: 4, 3: 3, 2: 2, 1: 1, 6: 6, 7: 7, 8: 8, 9: 9})
 
         fingerprint_1 = {0: {'self': 'C', 'C': [1, 2, 4], 'H': [11]},
                          1: {'self': 'C', 'C': [0, 3, 9], 'H': [12]},

arc/reaction/reaction.py

@@ -109,10 +109,13 @@ def __init__(self,
         self.kinetics = kinetics
         self.rmg_kinetics = None
         self.long_kinetic_description = ''
-        if check_family_name(family):
-            self.family = family
-        else:
-            raise ValueError(f"Invalid family name: {family}")
+        self._family = None
+        self._family_determined = False
+        if family is not None:
+            if check_family_name(family):
+                self.family = family
+            else:
+                raise ValueError(f"Invalid family name: {family}")
         self._family_own_reverse = False
         self.ts_label = ts_label
         self.dh_rxn298 = None
@@ -201,7 +204,7 @@ def multiplicity(self):
             if self._multiplicity is not None:
                 logger.info(f'Setting multiplicity of reaction {self.label} to {self._multiplicity}')
             else:
-                logger.Error(f'Could not determine multiplicity for the reaction: {self.label}')
+                logger.error(f'Could not determine multiplicity for the reaction: {self.label}')
         return self._multiplicity
 
     @multiplicity.setter
@@ -214,14 +217,16 @@ def multiplicity(self, value):
     @property
     def family(self):
         """The RMG reaction family"""
-        if self._family is None:
+        if not self._family_determined:
             self._family, self._family_own_reverse = self.determine_family()
+            self._family_determined = True
         return self._family
 
     @family.setter
     def family(self, value):
         """Allow setting family"""
         self._family = value
+        self._family_determined = True
         if value is not None and not isinstance(value, str):
             raise InputError(f'Reaction family must be a string, got {value} which is a {type(value)}.')
 

arc/species/species.py

@@ -463,7 +463,7 @@ def __init__(self,
                         self.multiplicity = self.mol.multiplicity
                     if self.charge is None:
                         self.charge = self.mol.get_net_charge()
-            if regen_mol:
+            if regen_mol and not (self.mol is not None and self.keep_mol):
                 # Perceive molecule from xyz coordinates. This also populates the .mol attribute of the Species.
                 # It overrides self.mol generated from adjlist or smiles so xyz and mol will have the same atom order.
                 if self.final_xyz or self.initial_xyz or self.most_stable_conformer or self.conformers or self.ts_guesses:
@@ -1992,13 +1992,17 @@ def _scissors(self,
                 sort_atoms_in_descending_label_order(split)
 
         if len(mol_splits) == 1:  # If cutting leads to only one split, then the split is cyclic.
+            mol1 = mol_splits[0]
+            self._assign_radicals_after_scission(mol=mol1)
+            mol1.update_multiplicity()
             spc1 = ARCSpecies(label=self.label + '_BDE_' + str(indices[0] + 1) + '_' + str(indices[1] + 1) + '_cyclic',
-                              mol=mol_splits[0],
-                              multiplicity=mol_splits[0].multiplicity,
-                              charge=mol_splits[0].get_net_charge(),
+                              mol=mol1,
+                              xyz=self.final_xyz,
+                              multiplicity=mol1.multiplicity,
+                              charge=mol1.get_net_charge(),
                               compute_thermo=False,
-                              e0_only=True)
-            spc1.generate_conformers()
+                              e0_only=True,
+                              keep_mol=True)
             return [spc1]
         elif len(mol_splits) == 2:
             mol1, mol2 = mol_splits
@@ -2020,19 +2024,7 @@ def _scissors(self,
 
         added_radical = list()
         for mol, label in zip([mol1, mol2], [label1, label2]):
-            for atom in mol.atoms:
-                theoretical_charge = elements.PeriodicSystem.valence_electrons[atom.symbol] \
-                                     - atom.get_total_bond_order() \
-                                     - atom.radical_electrons - \
-                                     2 * atom.lone_pairs
-                if theoretical_charge == atom.charge + 1:
-                    # we're missing a radical electron on this atom
-                    if label not in added_radical or label == 'H':
-                        atom.radical_electrons += 1
-                        added_radical.append(label)
-                    else:
-                        raise SpeciesError(f'Could not figure out which atom should gain a radical '
-                                           f'due to scission in {self.label}')
+            self._assign_radicals_after_scission(mol=mol, label=label, added_radical=added_radical)
         mol1.update(log_species=False, raise_atomtype_exception=False, sort_atoms=False)
         mol2.update(log_species=False, raise_atomtype_exception=False, sort_atoms=False)
 
@@ -2053,7 +2045,6 @@ def _scissors(self,
                           compute_thermo=False,
                           e0_only=True,
                           keep_mol=True)
-        spc1.generate_conformers()
         spc1.rotors_dict = None
         spc2 = ARCSpecies(label=label2,
                           mol=mol2,
@@ -2063,11 +2054,38 @@ def _scissors(self,
                           compute_thermo=False,
                           e0_only=True,
                           keep_mol=True)
-        spc2.generate_conformers()
         spc2.rotors_dict = None
 
         return [spc1, spc2]
 
+    def _assign_radicals_after_scission(self,
+                                        mol: Molecule,
+                                        label: str = None,
+                                        added_radical: list = None):
+        """
+        A helper function to assign radical electrons to atoms after scission.
+
+        Args:
+            mol (Molecule): The molecule to update.
+            label (str, optional): The label of the species.
+            added_radical (list, optional): A list of labels for which a radical was already added.
+        """
+        for atom in mol.atoms:
+            theoretical_charge = elements.PeriodicSystem.valence_electrons[atom.symbol] \
+                                 - atom.get_total_bond_order() \
+                                 - atom.radical_electrons - \
+                                 2 * atom.lone_pairs
+            if theoretical_charge == atom.charge + 1:
+                # we're missing a radical electron on this atom
+                if added_radical is None:
+                    atom.radical_electrons += 1
+                elif label not in added_radical or label == 'H':
+                    atom.radical_electrons += 1
+                    added_radical.append(label)
+                else:
+                    raise SpeciesError(f'Could not figure out which atom should gain a radical '
+                                       f'due to scission in {self.label}')
+
     def populate_ts_checks(self):
         """Populate (or restart) the .ts_checks attribute with default (``None``) values."""
         if self.is_ts: