Merge pull request #571 from avcopan/dev

Fix: Stereo-dependent reaction mapping from geometries

automol/reac/_3find.py

@@ -466,19 +466,28 @@ def additions(rct_gras, prd_gras):
         iso_dct = graph.isomorphism(rcts_gra_, prds_gra, stereo=False)
         if iso_dct:
             f_frm_bnd_key = (atm1_key, atm2_key)
-            b_brk_bnd_key = (iso_dct[atm1_key], iso_dct[atm2_key])
-            forw_tsg = ts.graph(rcts_gra, frm_bnd_keys=[f_frm_bnd_key], brk_bnd_keys=[])
-            back_tsg = ts.graph(prds_gra, frm_bnd_keys=[], brk_bnd_keys=[b_brk_bnd_key])
-
-            # Create the reaction object
-            rxn = from_forward_reverse(
-                cla=ReactionClass.ADDITION,
-                ftsg=forw_tsg,
-                rtsg=back_tsg,
-                rcts_keys=list(map(graph.atom_keys, rct_gras)),
-                prds_keys=list(map(graph.atom_keys, prd_gras)),
+            b_brk_bnd_key0 = (iso_dct[atm1_key], iso_dct[atm2_key])
+            b_brk_bnd_keys = graph.equivalent_bonds(
+                prds_gra, b_brk_bnd_key0, stereo=False, dummy=False
             )
-            rxns.append(rxn)
+
+            for b_brk_bnd_key in b_brk_bnd_keys:
+                forw_tsg = ts.graph(
+                    rcts_gra, frm_bnd_keys=[f_frm_bnd_key], brk_bnd_keys=[]
+                )
+                back_tsg = ts.graph(
+                    prds_gra, frm_bnd_keys=[], brk_bnd_keys=[b_brk_bnd_key]
+                )
+
+                # Create the reaction object
+                rxn = from_forward_reverse(
+                    cla=ReactionClass.ADDITION,
+                    ftsg=forw_tsg,
+                    rtsg=back_tsg,
+                    rcts_keys=list(map(graph.atom_keys, rct_gras)),
+                    prds_keys=list(map(graph.atom_keys, prd_gras)),
+                )
+                rxns.append(rxn)
 
     return rxns
 
@@ -767,7 +776,6 @@ def find(rct_gras, prd_gras, stereo=False):
     all_rxns = []
     for finder_ in finders_:
         rxns = finder_(rct_gras, prd_gras)
-        rxns = unique(rxns)
         for rxn in rxns:
             if not stereo:
                 all_rxns.append(rxn)
@@ -776,7 +784,8 @@ def find(rct_gras, prd_gras, stereo=False):
                     rxn, rct_gras0, prd_gras0, shift_keys=True
                 )
                 all_rxns.extend(srxns)
-
+    # Check for uniqueness *after* stereochemistry is assigned
+    all_rxns = unique(all_rxns)
     return tuple(all_rxns)
 
 

automol/tests/test_reac.py

@@ -1,10 +1,56 @@
-"""Test reac
-"""
+"""Test reac."""
 
 import pytest
+
 from automol import chi as chi_
 from automol import geom, graph, reac, smiles, zmat
 
+# Stereo-dependent reaction ID
+#   [CH]1[C@H]2CC[C@@H]1O2 => C1=C2CC[C@@H]1O2 + [H]
+#         *       ^              *    ^
+#
+# Notes:
+#   1. Without stereochemistry, * and ^ are symmetrically equivalent, so that the
+#      C-H bond of either one could be broken to yield the product.
+#   2. With stereochemistry, this is not the case. Only one of them can yield the product.
+#   3. This was causing a bug in our reaction mapping strategy, as only one of
+#      the (apparently equivalent) bonds was enumerated for the beta scission,
+#      but it was for the wrong (impossible) stereoisomer.
+C5H7O_RGEOS = (
+    (
+        ("C", (-1.818776, 0.000365, -1.870351)),
+        ("C", (-0.728737, 1.849057, 0.019341)),
+        ("H", (-3.815204, 0.00202, -2.406224)),
+        ("O", (-1.43218, -0.00102, 1.941119)),
+        ("C", (2.14661, 1.467934, -0.165888)),
+        ("H", (-1.480623, 3.753408, 0.344455)),
+        ("C", (-0.7282, -1.849316, 0.017276)),
+        ("H", (-1.479758, -3.753888, 0.342001)),
+        ("C", (2.147088, -1.467251, -0.165725)),
+        ("H", (2.945064, 2.30276, -1.89216)),
+        ("H", (3.116856, 2.276812, 1.482364)),
+        ("H", (3.11582, -2.275557, 1.483743)),
+        ("H", (2.947381, -2.302127, -1.891058)),
+    ),
+)
+C5H7O_PGEOS = (
+    (
+        ("C", (-1.656404, -1.057905, 1.788813)),
+        ("C", (-1.32713, 1.452595, 0.301729)),
+        ("H", (-1.047337, -1.533683, 3.704721)),
+        ("O", (-1.448161, 0.042876, -2.01624)),
+        ("C", (1.639596, 1.769402, 0.474459)),
+        ("H", (-2.548583, 3.132631, 0.445664)),
+        ("C", (-0.411648, -1.765102, -0.368446)),
+        ("C", (2.371036, -1.082399, -0.302428)),
+        ("H", (2.362823, 3.238507, -0.810887)),
+        ("H", (2.320699, 2.182828, 2.397452)),
+        ("H", (3.469656, -2.054244, 1.162259)),
+        ("H", (3.335331, -1.208571, -2.134068)),
+    ),
+    (("H", (0.0, 0.0, 0.0)),),
+)
+
 
 def test__reactant_graphs():
     """Test reac.reactant_graphs"""
@@ -231,6 +277,19 @@ def _test(rct_smis, prd_smis):
     _test([r"F\N=[C]/F", "[C]#C"], [r"F\N=C(C#C)/F"])
 
 
+@pytest.mark.parametrize(
+    "fml,rgeos,pgeos,nrxns",
+    [
+        ("C5H7O", C5H7O_RGEOS, C5H7O_PGEOS, 1),
+    ],
+)
+def test__from_geometries(fml, rgeos, pgeos, nrxns):
+    print(f"Testing for {fml}")
+    rxns = reac.from_geometries(rct_geos=rgeos, prd_geos=pgeos, stereo=True)
+    print(len(rxns))
+    assert len(rxns) == nrxns, f"{len(rxns)} != {nrxns}"
+
+
 @pytest.mark.parametrize(
     "rct_smi,prd_smi",
     [
@@ -456,5 +515,6 @@ def test__canonical_enantiomer():
     # test__end_to_end("[C@H](O)(C)F.[Cl]", "[C@@H](O)(C)Cl.[F]")
     # test__end_to_end("CCCO[O]", "[CH2]CCOO")
     # test__end_to_end("C[C]1O[C@H]1COO", r"C/C([O])=C\COO")
-    test__end_to_end("CC1[C](O1)COO", "CC1C2(O1)CO2.[OH]")
+    # test__end_to_end("CC1[C](O1)COO", "CC1C2(O1)CO2.[OH]")
     # test__canonical_enantiomer()
+    test__from_geometries("C5H7O", C5H7O_RGEOS, C5H7O_PGEOS, 1)