TMP

Author: alongd

arc/job/adapters/ts/linear.py

@@ -221,40 +221,41 @@ def execute_incore(self):
                                                           )
 
             t0_0 = datetime.datetime.now()
-            xyz_0 = interpolate(rxn=rxn, use_weights=False)
+            xyzs_0 = interpolate(rxn=rxn, use_weights=False)
             t_ex_0 = datetime.datetime.now() - t0_0
 
             t0_1 = datetime.datetime.now()
-            xyz_1 = interpolate(rxn=rxn, use_weights=True)
+            xyzs_1 = interpolate(rxn=rxn, use_weights=True)
             t_ex_1 = datetime.datetime.now() - t0_1
 
-            for method_index, (xyz, t0, t_ex) in enumerate(zip([xyz_0, xyz_1], [t0_0, t0_1], [t_ex_0, t_ex_1])):
-                if colliding_atoms(xyz):
-                    continue
-                unique = True
-                for other_tsg in rxn.ts_species.ts_guesses:
-                    if almost_equal_coords(xyz, other_tsg.initial_xyz):
-                        if 'linear' not in other_tsg.method.lower():
-                            other_tsg.method += f' and Linear {method_index}'
-                        unique = False
-                        break
-                if unique:
-                    ts_guess = TSGuess(method=f'linear {method_index}',
-                                       index=len(rxn.ts_species.ts_guesses),
-                                       method_index=method_index,
-                                       t0=t0,
-                                       execution_time=t_ex,
-                                       success=True,
-                                       family=family_label,
-                                       xyz=xyz,
-                                       )
-                    rxn.ts_species.ts_guesses.append(ts_guess)
-                    save_geo(xyz=xyz,
-                             path=self.local_path,
-                             filename=f'Linear {method_index}',
-                             format_='xyz',
-                             comment=f'Linear {method_index}, family: {family_label}',
-                             )
+            for method_index, (xyzs, t0, t_ex) in enumerate(zip([xyzs_0, xyzs_1], [t0_0, t0_1], [t_ex_0, t_ex_1])):
+                for xyz in xyzs:
+                    if colliding_atoms(xyz):
+                        continue
+                    unique = True
+                    for other_tsg in rxn.ts_species.ts_guesses:
+                        if almost_equal_coords(xyz, other_tsg.initial_xyz):
+                            if 'linear' not in other_tsg.method.lower():
+                                other_tsg.method += f' and Linear {method_index}'
+                            unique = False
+                            break
+                    if unique:
+                        ts_guess = TSGuess(method=f'linear {method_index}',
+                                           index=len(rxn.ts_species.ts_guesses),
+                                           method_index=method_index,
+                                           t0=t0,
+                                           execution_time=t_ex,
+                                           success=True,
+                                           family=family_label,
+                                           xyz=xyz,
+                                           )
+                        rxn.ts_species.ts_guesses.append(ts_guess)
+                        save_geo(xyz=xyz,
+                                 path=self.local_path,
+                                 filename=f'Linear {method_index}',
+                                 format_='xyz',
+                                 comment=f'Linear {method_index}, family: {family_label}',
+                                 )
 
             if len(self.reactions) < 5:
                 successes = len([tsg for tsg in rxn.ts_species.ts_guesses if tsg.success and 'linear' in tsg.method])
@@ -275,7 +276,7 @@ def execute_queue(self):
 
 def interpolate(rxn: 'ARCReaction',
                 use_weights: bool = False,
-                ) -> Optional[dict]:
+                ) -> Optional[List[dict]]:
     """
     Search for a TS by interpolating internal coords.
 
@@ -284,7 +285,7 @@ def interpolate(rxn: 'ARCReaction',
         use_weights (bool, optional): Whether to use the well energies to determine relative interpolation weights.
 
     Returns:
-        Optional[dict]: The XYZ coordinates guess.
+        Optional[List[dict]]: Entries are the XYZ coordinate guesses.
     """
     if rxn.is_isomerization():
         return interpolate_isomerization(rxn=rxn, use_weights=use_weights)
@@ -295,7 +296,7 @@ def interpolate(rxn: 'ARCReaction',
 
 def interpolate_isomerization(rxn: 'ARCReaction',
                               use_weights: bool = False,
-                              ) -> Optional[dict]:
+                              ) -> Optional[List[dict]]:
     """
     Search for a TS of an isomerization reaction by interpolating internal coords.
 
@@ -304,23 +305,25 @@ def interpolate_isomerization(rxn: 'ARCReaction',
         use_weights (bool, optional): Whether to use the well energies to determine relative interpolation weights.
 
     Returns:
-        Optional[dict]: The XYZ coordinates guess.
+        Optional[List[dict]]: Entries are the XYZ coordinate guesses.
     """
     rxn.r_species[0].get_xyz()
     rxn.p_species[0].get_xyz()
     r_zmat = xyz_to_zmat(xyz=rxn.r_species[0].get_xyz(),
                          consolidate=False,
                          atom_order=list(range(sum(r.number_of_atoms for r in rxn.r_species))))
-    p_zmat = xyz_to_zmat(xyz=rxn.p_species[0].get_xyz(),
-                         consolidate=False,
-                         atom_order=rxn.atom_map)
     weight = get_rxn_weight(rxn) if use_weights else 0.5
-    if weight is None:
-        return None
-    ts_zmat = average_zmat_params(zmat_1=r_zmat, zmat_2=p_zmat, weight=weight)
-    if ts_zmat is None:
+    if weight is None or rxn.atom_maps is None:
         return None
-    return zmat_to_xyz(ts_zmat)
+    ts_xyzs = list()
+    for atom_map in rxn.atom_maps:
+        p_zmat = xyz_to_zmat(xyz=rxn.p_species[0].get_xyz(),
+                             consolidate=False,
+                             atom_order=atom_map)
+        ts_zmat = average_zmat_params(zmat_1=r_zmat, zmat_2=p_zmat, weight=weight)
+        if ts_zmat is not None:
+            ts_xyzs.append(zmat_to_xyz(ts_zmat))
+    return ts_xyzs
 
 
 def average_zmat_params(zmat_1: dict,

arc/job/adapters/ts/linear_test.py

@@ -22,9 +22,9 @@
 from arc.species.zmat import _compare_zmats
 
 
-class TestHeuristicsAdapter(unittest.TestCase):
+class TestLinearAdapter(unittest.TestCase):
     """
-    Contains unit tests for the HeuristicsAdapter class.
+    Contains unit tests for the LinearAdapter class.
     """
 
     @classmethod
@@ -129,7 +129,7 @@ def test_get_rxn_weight(self):
         rxn_1.ts_species.e0 = 391.6
         self.assertAlmostEquals(get_rxn_weight(rxn_1), 0.3417832)
 
-    def test_interpolate_isomerization(self):
+    def test_interpolate_intra_h_migration(self):
         """Test the interpolate_isomerization() function."""
         nc3h7_xyz = """C                  0.00375165   -0.48895802   -1.20586379
                        C                  0.00375165   -0.48895802    0.28487510
@@ -164,8 +164,9 @@ def test_interpolate_isomerization(self):
                                         H       0.76979130    1.33747945    0.33815513
                                         H      -0.04544494    0.70455273    1.77835334
                                         H      -1.00071642    1.24557408    0.38839197""")
-        ts_xyz = interpolate_isomerization(rxn, use_weights=False)
-        self.assertTrue(almost_equal_coords(ts_xyz, expected_ts_xyz))
+        ts_xyzs = interpolate_isomerization(rxn, use_weights=False)
+        self.assertEqual(len(ts_xyzs), 2)
+        self.assertTrue(almost_equal_coords(ts_xyzs[0], expected_ts_xyz))
 
         nc3h7.e0 = 101.55
         ic3h7.e0 = 88.91
@@ -182,8 +183,51 @@ def test_interpolate_isomerization(self):
                                         H       0.79813573    1.38347069    0.43772483
                                         H      -0.03897336    0.76031233    1.86961141
                                         H      -0.97425159    1.33180895    0.47825005""")
-        ts_xyz = interpolate_isomerization(rxn, use_weights=True)
-        self.assertTrue(almost_equal_coords(ts_xyz, expected_ts_xyz))
+        ts_xyzs = interpolate_isomerization(rxn, use_weights=True)
+        self.assertTrue(almost_equal_coords(ts_xyzs[0], expected_ts_xyz))
+
+    def test_interpolate_no2_ono_isomerization(self):
+        """Test the interpolate_isomerization() function."""
+        c2h5no2_xyz = """C      -1.12362739   -0.04664655   -0.08575959
+                         C       0.24488022   -0.51587553    0.36119196
+                         N       0.57726975   -1.77875156   -0.37104243
+                         O       1.16476543   -1.66382529   -1.45384186
+                         O       0.24561669   -2.84385320    0.16410116
+                         H      -1.87655344   -0.80826847    0.13962125
+                         H      -1.14729169    0.14493421   -1.16405294
+                         H      -1.41423043    0.87863077    0.42354512
+                         H       1.02430791    0.21530309    0.12674144
+                         H       0.27058353   -0.73979548    1.43184405"""
+        c2h5ono_xyz = """C      -1.36894499    0.07118059   -0.24801399
+                         C      -0.01369535    0.17184136    0.42591278
+                         O      -0.03967083   -0.62462610    1.60609048
+                         N       1.23538512   -0.53558048    2.24863846
+                         O       1.25629155   -1.21389295    3.27993827
+                         H      -2.16063255    0.41812452    0.42429392
+                         H      -1.39509985    0.66980796   -1.16284741
+                         H      -1.59800183   -0.96960842   -0.49986392
+                         H       0.19191326    1.21800574    0.68271847
+                         H       0.76371340   -0.19234475   -0.25650067"""
+        c2h5no2 = ARCSpecies(label='C2H5NO2', smiles='CC[N+](=O)[O-]', xyz=c2h5no2_xyz)
+        c2h5ono = ARCSpecies(label='C2H5ONO', smiles='CCON=O', xyz=c2h5ono_xyz)
+
+        rxn = ARCReaction(r_species=[c2h5no2], p_species=[c2h5ono])
+        expected_ts_xyz = str_to_xyz("""C       0.01099731   -0.46789926   -1.15958911
+                                        C       0.01099731   -0.46789926    0.33114978
+                                        C       0.01099731    0.94103865    0.90031155
+                                        H       0.57795661   -1.24174248   -1.65467180
+                                        H      -0.39690222    0.34527841   -1.69240298
+                                        H      -1.19440431   -1.28933062   -0.47327539
+                                        H       0.89689057   -1.16420498    0.45967951
+                                        H       0.76979130    1.33747945    0.33815513
+                                        H      -0.04544494    0.70455273    1.77835334
+                                        H      -1.00071642    1.24557408    0.38839197""")
+        ts_xyzs = interpolate_isomerization(rxn, use_weights=False)
+        for xyz in ts_xyzs:
+            print('\n\n\n')
+            print(xyz_to_str(xyz))
+        self.assertEqual(len(ts_xyzs), 1)
+        self.assertTrue(almost_equal_coords(ts_xyzs[0], expected_ts_xyz))
 
     @classmethod
     def tearDownClass(cls):

arc/mapping/driver.py

@@ -257,11 +257,13 @@ def map_rxn(rxn: 'ARCReaction',
             logger.error(f"Could not find isomorphism for scissored species: {[cut.mol.smiles for cut in p_cuts]}")
             return None
 
-        atom_map = map_pairs(pairs_of_reactant_and_products)
-        if atom_map is not None:
-            atom_map = glue_maps(atom_map, pairs_of_reactant_and_products)
+        atom_map_list = map_pairs(pairs_of_reactant_and_products)
+        if isinstance(atom_map_list, list) and all(e is not None for e in atom_map_list):
+            atom_map = glue_maps(atom_map_list, pairs_of_reactant_and_products)
             atom_maps.append(atom_map)
 
+    print(f'atom_maps: {atom_maps}')
+
     return atom_maps
 
 

arc/mapping/driver_test.py

@@ -1047,6 +1047,32 @@ def test_map_intra_h_migration(self):
         self.assertIn(atom_map[7], [3, 4, 5, 8])
         self.assertIn(atom_map[8], [3, 4, 5, 8])
 
+        ic3h7_xyz = """C                 -0.40735690   -0.74240205   -0.34312948
+                       C                  0.38155377   -0.25604705    0.82450968
+                       C                  0.54634593    1.25448345    0.81064511
+                       H                  0.00637731   -1.58836501   -0.88041673
+                       H                 -0.98617584   -0.01198912   -0.89732723
+                       H                 -1.29710684   -1.29092340    0.08598983
+                       H                  1.36955428   -0.72869684    0.81102246
+                       H                  1.06044877    1.58846788   -0.09702437
+                       H                  1.13774084    1.57830484    1.67308862
+                       H                 -0.42424546    1.75989927    0.85794283"""
+        nc3h7_xyz = """C                  0.00375165   -0.48895802   -1.20586379
+                       C                  0.00375165   -0.48895802    0.28487510
+                       C                  0.00375165    0.91997987    0.85403684
+                       H                  0.41748586   -1.33492098   -1.74315104
+                       H                 -0.57506729    0.24145491   -1.76006154
+                       H                 -0.87717095   -1.03203740    0.64280162
+                       H                  0.88948616   -1.02465371    0.64296621
+                       H                  0.88512433    1.48038223    0.52412379
+                       H                  0.01450405    0.88584135    1.94817394
+                       H                 -0.88837301    1.47376959    0.54233121"""
+        rxn_1 = ARCReaction(r_species=[ARCSpecies(label='iC3H7', smiles='C[CH]C', xyz=ic3h7_xyz)],
+                            p_species=[ARCSpecies(label='nC3H7', smiles='[CH2]CC', xyz=nc3h7_xyz)])
+        print('---------')
+        print(rxn_1.atom_maps[0])
+        self.assertEqual(rxn_1.atom_maps[0], [0, 1, 2, 6, 3, 4, 5, 7, 8, 9])
+
     def test_map_isomerization_reaction(self):
         """Test the map_isomerization_reaction() function."""
         reactant_xyz = """C  -1.3087    0.0068    0.0318

arc/mapping/engine.py

@@ -1088,23 +1088,18 @@ def make_bond_changes(rxn: 'ARCReaction',
 def assign_labels_to_products(rxn: 'ARCReaction',
                               p_label_dict: dict):
     """
-    Add the indices to the reactants and products.
+    Add labels to the atoms of the reactants and products.
 
     Args:
         rxn: ARCReaction object to be mapped
         p_label_dict: the labels of the products
-        Consider changing in rmgpy.
-
-    Returns:
-        Adding labels to the atoms of the reactants and products, to be identified later.
     """
-
     atom_index = 0
-    for product in rxn.p_species:
-        for atom in product.mol.atoms:
+    for product_ in rxn.p_species:
+        for atom in product_.mol.atoms:
             if atom_index in p_label_dict.values() and (atom.label is str or atom.label is None):
                 atom.label = key_by_val(p_label_dict,atom_index)
-            atom_index+=1
+            atom_index += 1
 
 
 def update_xyz(spcs: List[ARCSpecies]) -> List[ARCSpecies]:
@@ -1169,7 +1164,7 @@ def pairing_reactants_and_products_for_mapping(r_cuts: List[ARCSpecies],
 
 def map_pairs(pairs):
     """
-    A function that maps the mached species together
+    A function that maps the matched species together.
 
     Args:
         pairs: A list of the pairs of reactants and species
@@ -1181,7 +1176,6 @@ def map_pairs(pairs):
     maps = list()
     for pair in pairs:
         maps.append(map_two_species(pair[0], pair[1]))
-
     return maps
 
 
@@ -1196,20 +1190,23 @@ def label_species_atoms(spcs):
     for spc in spcs:
         for atom in spc.mol.atoms:
             atom.label = str(index)
-            index+=1
+            index += 1
 
 
-def glue_maps(maps, pairs_of_reactant_and_products):
+def glue_maps(maps: List[List[int]],
+              pairs_of_reactant_and_products: List[Tuple[ARCSpecies, ARCSpecies]]
+              ) -> List[int]:
     """
     a function that joins together the maps from the parts of the reaction.
 
     Args:
-        rxn: ARCReaction that requires atom mapping
-        maps: The list of all maps of the isomorphic cuts.
+        maps (List[List[int]]): The list of all maps of the isomorphic cuts.
+        pairs_of_reactant_and_products (List[Tuple[ARCSpecies, ARCSpecies]]): The list of all pairs of reactants and products.
 
     Returns:
-        an Atom Map of the compleate reaction.
+        List[int]: An Atom Map of the complete reaction.
     """
+    # print(f'maps: {maps}, pairs_of_reactant_and_products: {pairs_of_reactant_and_products}')
     am_dict = dict()
     for _map, pair in zip(maps, pairs_of_reactant_and_products):
         r_atoms = pair[0].mol.atoms
@@ -1221,16 +1218,16 @@ def glue_maps(maps, pairs_of_reactant_and_products):
 
 def determine_bdes_on_spc_based_on_atom_labels(spc: "ARCSpecies", bde: Tuple[int, int]) -> bool:
     """
-    A function for determining whether or not the species in question containt the bond specified by the bond dissociation indices.
+    A function for determining whether the species in question contains the bond specified by the bond dissociation indices.
     Also, assigns the correct BDE to the species.
     
     Args:
         spc (ARCSpecies): The species in question, with labels atom indices.
         bde (Tuple[int, int]): The bde in question.
-        add_bdes (bool): Whether or not to add the bde to the species.
+        add_bdes (bool): Whether to add the bde to the species.
     
     Returns:
-        bool: Whether or not the bde is based on the atom labels.
+        bool: Whether the bde is based on the atom labels.
     """
     bde = convert_list_index_0_to_1(bde, direction=-1)
     index1, index2 = bde[0], bde[1]

arc/mapping/engine_test.py

@@ -698,7 +698,7 @@ def test_glue_maps(self):
         rxn_1_test.determine_family(self.db)
         rmg_reactions = get_rmg_reactions_from_arc_reaction(arc_reaction=rxn_1_test, backend="ARC")
         r_label_dict, p_label_dict = get_atom_indices_of_labeled_atoms_in_an_rmg_reaction(arc_reaction=rxn_1_test,
-                                                                                      rmg_reaction=rmg_reactions[0])
+                                                                                          rmg_reaction=rmg_reactions[0])
         assign_labels_to_products(rxn_1_test, p_label_dict)
         reactants, products = copy_species_list_for_mapping(rxn_1_test.r_species), copy_species_list_for_mapping(rxn_1_test.p_species)
         label_species_atoms(reactants), label_species_atoms(products)
@@ -709,7 +709,7 @@ def test_glue_maps(self):
         p_cuts = cut_species_based_on_atom_indices(products, p_bdes)
         pairs_of_reactant_and_products = pairing_reactants_and_products_for_mapping(r_cuts, p_cuts)
         maps = map_pairs(pairs_of_reactant_and_products)
-        atom_map = glue_maps(maps,pairs_of_reactant_and_products)
+        atom_map = glue_maps(maps, pairs_of_reactant_and_products)
         self.assertEqual(len(atom_map), self.r_1.mol.get_num_atoms() + self.r_2.mol.get_num_atoms())
         atoms_r = [atom for atom in self.r_1.mol.copy(deep=True).atoms] + [atom for atom in self.r_2.mol.copy(deep=True).atoms]
         atoms_p = [atom for atom in self.p_1.mol.copy(deep=True).atoms] + [atom for atom in self.p_2.mol.copy(deep=True).atoms]