Fix: Stereo-dependent reaction mapping from geometries

Problem case:
    [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.

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)