Merge pull request #593 from avcopan/dev

Fix: Stop including strained structures by default

automol/graph/base/_02algo.py

@@ -707,9 +707,16 @@ def branch_dict(gra, atm_key, keep_root=False, stereo=False):
     return bnch_dct
 
 
-def branch(gra, atm_key, branch_key, keep_root=False, stereo=False):
+def branch(
+    gra,
+    atm_key,
+    branch_key,
+    keep_root: bool = False,
+    stereo: bool = False,
+    return_missing_neighbor: bool = False,
+):
     """Get the atom keys for a branch off of an atom extending a long a
-    specific neighboring atom or bond
+    specific neighboring atom or bond.
 
     :param gra: molecular graph
     :type gra: automol graph data structure
@@ -718,9 +725,9 @@ def branch(gra, atm_key, branch_key, keep_root=False, stereo=False):
     :param branch_key: the adjacent atom or bond key that begins the branch
     :type branch_key: int or frozenset[int]
     :param keep_root: Keep root atom as part of branch? defaults to False
-    :type keep_root: bool, optional
     :param stereo: Keep stereochemistry in the subgraph? defaults to False
-    :type stereo: bool, optional
+    :param return_missing_neighbor: If the branch key is not a neighbor, return it as a
+        one-atom branch?
     :return: The branch atom keys
     :rtype: frozenset[int]
     """
@@ -734,10 +741,19 @@ def branch(gra, atm_key, branch_key, keep_root=False, stereo=False):
         f"Input graph:\n{string(gra)}"
     )
     bnch_dct = branch_dict(gra, atm_key, keep_root=keep_root, stereo=stereo)
+    if return_missing_neighbor:
+        return bnch_dct.get(branch_key, subgraph_(gra, {branch_key}, stereo=stereo))
+
     return bnch_dct[branch_key]
 
 
-def branch_atom_keys(gra, atm_key, branch_key, keep_root=False):
+def branch_atom_keys(
+    gra,
+    atm_key,
+    branch_key,
+    keep_root: bool = False,
+    return_missing_neighbor: bool = False,
+):
     """Get the atom keys for a branch off of an atom extending a long a
     specific neighboring atom or bond
 
@@ -748,11 +764,20 @@ def branch_atom_keys(gra, atm_key, branch_key, keep_root=False):
     :param branch_key: the adjacent atom or bond key that begins the branch
     :type branch_key: int or frozenset[int]
     :param keep_root: Keep root atom as part of branch? defaults to False
-    :type keep_root: bool, optional
+    :param return_missing_neighbor: If the branch key is not a neighbor, return it as a
+        one-atom branch?
     :return: The branch atom keys
     :rtype: frozenset[int]
     """
-    return atom_keys(branch(gra, atm_key, branch_key, keep_root=keep_root))
+    return atom_keys(
+        branch(
+            gra,
+            atm_key,
+            branch_key,
+            keep_root=keep_root,
+            return_missing_neighbor=return_missing_neighbor,
+        )
+    )
 
 
 def is_branched(gra):

automol/graph/base/_05stereo.py

@@ -31,6 +31,7 @@
     rigid_planar_bonds_with_ring_constraints,
 )
 from ._04class import (
+    atom_transfers,
     insertions,
     substitutions,
     vinyl_addition_candidates,
@@ -208,10 +209,8 @@ def _bond_is_stereogenic(nkeys):
 
 
 def stereoatom_bridgehead_pairs(
-    gra, cand_dct: CenterNeighborDict | None = None
-) -> dict[
-    AtomKey, tuple[tuple[AtomKey, AtomKey, AtomKey], tuple[AtomKey, AtomKey, AtomKey]]
-]:
+    gra, cand_dct: CenterNeighborDict | None = None, migration_only: bool = False
+) -> dict[tuple[AtomKey, AtomKey], tuple[AtomKeys, AtomKeys]]:
     r"""Identify pairs of interdependent bridgehead stereoatoms, if any.
 
     Bridgehead stereoatoms sharing the same three bridges are interdependent -- the
@@ -226,6 +225,7 @@ def stereoatom_bridgehead_pairs(
     :param gra: A molecular graph
     :param cand_dct: A mapping of stereocenter candidates onto their neighbor keys
         (To avoid recalculating)
+    :param migration_only: Only detect bridgehead pairs for atom migrations?
     :return: A dictionary mapping pairs of bridgehead atoms onto the connected neighbors
         of each, respectively
     """
@@ -244,6 +244,10 @@ def stereoatom_bridgehead_pairs(
             conn_nkeys1, conn_nkeys2 = zip(*sorted(conn_dct.items()), strict=True)
             bhp_dct[(key1, key2)] = (conn_nkeys1, conn_nkeys2)
 
+    if migration_only:
+        mig_pairs = list(map(sorted, atom_transfers(gra).values()))
+        bhp_dct = dict_.filter_by_key(bhp_dct, lambda k: sorted(k) in mig_pairs)
+
     return bhp_dct
 
 

automol/graph/base/_11stereo.py

@@ -9,6 +9,7 @@
 
 import more_itertools as mit
 import numpy
+
 from phydat import phycon
 
 from ... import util
@@ -98,8 +99,11 @@ def expand_stereo(
     )
 
     # 2. If requested, filter out strained stereoisomers
-    if not strained:
-        gps = _remove_strained_stereoisomers_from_expansion(gps, cand_dct=cand_dct)
+    # (If not, filter out only strained ring H-migrations, which are guaranteed to cause
+    # problems)
+    gps = _remove_strained_stereoisomers_from_expansion(
+        gps, cand_dct=cand_dct, migration_only=strained
+    )
 
     # 3. If requested, filter out non-canonical enantiomers
     if not enant:
@@ -196,14 +200,19 @@ def _select_ts_canonical_direction_priorities(gprs, fcand_dct: dict, rcand_dct:
     return gps
 
 
-def _remove_strained_stereoisomers_from_expansion(gps, cand_dct):
-    """Remove strained stereoisomers from an expansion"""
+def _remove_strained_stereoisomers_from_expansion(
+    gps, cand_dct, migration_only: bool = False
+):
+    """Remove strained stereoisomers from an expansion.
+
+    :param migration_only: Only remove strained H-migration bridgeheads?
+    """
     if not gps:
         return gps
 
     gps0 = list(gps)
     gra = without_stereo(gps0[0][0])
-    bhp_dct = stereoatom_bridgehead_pairs(gra, cand_dct)
+    bhp_dct = stereoatom_bridgehead_pairs(gra, cand_dct, migration_only=migration_only)
 
     gps = gps0.copy()
     for gra, pri_dct in gps0:
@@ -646,7 +655,12 @@ def geometry_pseudorotate_atom(
     # forming ring
     if nkey1 in rot_keys or nkey2 in rot_keys:
         rot_keys = set(
-            itertools.chain(*(branch_atom_keys(gra_reac, key, k) for k in rot_nkeys))
+            itertools.chain(
+                *(
+                    branch_atom_keys(gra_reac, key, k, return_missing_neighbor=True)
+                    for k in rot_nkeys
+                )
+            )
         )
 
     geo = geom_base.rotate(geo, rot_axis, ang, orig_xyz=xyz, idxs=rot_keys)

automol/reac/_5conv.py

@@ -30,7 +30,7 @@ def from_graphs(
     stereo: bool = True,
     struc_typ: str | None = None,
     enant: bool = True,
-    strained: bool = True,
+    strained: bool = False,
 ) -> tuple[Reaction, ...]:
     """Get reaction objects from graphs.
 
@@ -58,7 +58,7 @@ def from_amchis(
     stereo: bool = True,
     struc_typ: str | None = None,
     enant: bool = True,
-    strained: bool = True,
+    strained: bool = False,
 ) -> tuple[Reaction, ...]:
     """Get reaction objects from AMChIs.
 
@@ -86,7 +86,7 @@ def from_inchis(
     stereo: bool = True,
     struc_typ: str | None = None,
     enant: bool = True,
-    strained: bool = True,
+    strained: bool = False,
 ) -> tuple[Reaction, ...]:
     """Get reaction objects from InChIs.
 
@@ -114,7 +114,7 @@ def from_chis(
     stereo: bool = True,
     struc_typ: str | None = None,
     enant: bool = True,
-    strained: bool = True,
+    strained: bool = False,
 ) -> tuple[Reaction, ...]:
     """Get reaction objects from ChIs.
 
@@ -142,7 +142,7 @@ def from_smiles(
     stereo: bool = True,
     struc_typ: str | None = None,
     enant: bool = True,
-    strained: bool = True,
+    strained: bool = False,
 ) -> tuple[Reaction, ...]:
     """Get reaction objects from SMILES.
 
@@ -170,7 +170,7 @@ def from_geometries(
     stereo: bool = True,
     struc_typ: str | None = "geom",
     enant: bool = True,
-    strained: bool = True,
+    strained: bool = False,
 ) -> tuple[Reaction, ...]:
     """Get reaction objects from geometries.
 
@@ -202,7 +202,7 @@ def from_zmatrices(
     stereo: bool = True,
     struc_typ: str | None = "zmat",
     enant: bool = True,
-    strained: bool = True,
+    strained: bool = False,
 ) -> tuple[Reaction, ...]:
     """Get reaction objects from z-matrices.
 

automol/tests/test_graph_ts.py

@@ -1,9 +1,10 @@
 """ test graph.ts
 """
 
-import automol
 import numpy
 import pytest
+
+import automol
 from automol import graph
 
 # Sn2 Atom Stereo
@@ -906,6 +907,54 @@
 )
 
 
+# Simple Bridgehead Transfer
+C5H7_TSG = (
+    {
+        0: ("C", 0, True),
+        1: ("C", 0, None),
+        2: ("H", 0, None),
+        3: ("C", 0, False),
+        4: ("H", 0, None),
+        5: ("C", 0, None),
+        6: ("H", 0, None),
+        7: ("H", 0, None),
+        8: ("C", 0, None),
+        9: ("H", 0, None),
+        10: ("H", 0, None),
+        11: ("H", 0, None),
+    },
+    {
+        frozenset({1, 4}): (1, None),
+        frozenset({8, 11}): (1, None),
+        frozenset({0, 6}): (0.1, None),
+        frozenset({3, 7}): (1, None),
+        frozenset({0, 1}): (1, None),
+        frozenset({0, 2}): (1, None),
+        frozenset({3, 6}): (0.9, None),
+        frozenset({3, 5}): (1, None),
+        frozenset({10, 5}): (1, None),
+        frozenset({8, 5}): (1, None),
+        frozenset({1, 3}): (1, None),
+        frozenset({9, 5}): (1, None),
+        frozenset({0, 8}): (1, None),
+    },
+)
+C5H7_TS_GEO = (
+    ("C", (2.411517, 0.003628, 0.0)),
+    ("C", (0.9272, -2.162191, 0.0)),
+    ("H", (4.474794, 0.007151, -0.0)),
+    ("C", (-1.830041, -1.468223, -0.0)),
+    ("H", (1.626336, -4.101193, 0.0)),
+    ("C", (-1.834561, 1.463071, 0.0)),
+    ("H", (-2.817187, -2.252453, -1.662617)),
+    ("H", (-2.81719, -2.252455, 1.662611)),
+    ("C", (0.920338, 2.16457, -0.0)),
+    ("H", (-2.82419, 2.244254, 1.66257)),
+    ("H", (-2.824192, 2.244256, -1.662566)),
+    ("H", (1.614892, 4.105309, -2e-06)),
+)
+
+
 @pytest.mark.parametrize(
     "formula,tsg,geo,npars1,npars2",
     [
@@ -1014,6 +1063,7 @@ def _test(formula, tsg):
         ("C8H14", C8H14_TSG, [1, 1, 1, 1]),
         ("C5H7O2", C5H7O2_TSG, [1, 1]),
         ("C5H7O3", C5H7O3_TSG, [1, 1, 1, 1]),
+        ("C5H7", C5H7_TSG, [2]),
     ],
 )
 def test__ts__expand_reaction_stereo(formula, ts_gra, ts_counts):
@@ -1527,6 +1577,16 @@ def test__stereoatom_bridgehead_pairs(tsg, bh_dct0):
     assert bh_dct == bh_dct0, f"{bh_dct} != {bh_dct0}"
 
 
+@pytest.mark.parametrize("tsg0, ts_geo0", [(C5H7_TSG, C5H7_TS_GEO)])
+def test__stereo_corrected_geometry(tsg0, ts_geo0):
+    # Invert the stereochemistry so that correction of atoms is required
+    ts_geo0 = automol.geom.reflect_coordinates(ts_geo0)
+
+    ts_geo = automol.graph.stereo_corrected_geometry(tsg0, ts_geo0)
+    tsg = automol.graph.set_stereo_from_geometry(tsg0, ts_geo)
+    assert tsg == tsg0, f"   {tsg}\n!= {tsg0}"
+
+
 if __name__ == "__main__":
     # test__set_stereo_from_geometry()
     # test__to_local_stereo()
@@ -1545,6 +1605,8 @@ def test__stereoatom_bridgehead_pairs(tsg, bh_dct0):
     # test__ts__reagents_graph(
     #     "C5H7O_B", C5H7O_B_TSG, {3: True, 4: False}, {3: False, 4: True}
     # )
-    test__stereoatom_bridgehead_pairs(
-        C5H7O_B_TSG, {(1, 0): ((2, 3, 8), (4, 3, 8)), (3, 4): ((0, 1, 5), (0, 2, 5))}
-    )
+    # test__stereoatom_bridgehead_pairs(
+    #     C5H7O_B_TSG, {(1, 0): ((2, 3, 8), (4, 3, 8)), (3, 4): ((0, 1, 5), (0, 2, 5))}
+    # )
+    # test__stereo_corrected_geometry(C5H7_TSG, C5H7_TS_GEO)
+    test__ts__expand_reaction_stereo("C5H7", C5H7_TSG, [2])