From 9ed1554055b77b07a1d466a6608d0d4027e9ea51 Mon Sep 17 00:00:00 2001 From: avcopan Date: Thu, 12 Dec 2024 11:38:25 -0600 Subject: [PATCH 1/3] Fix: Stereo correction for ring H-migrations Includes a test case that fails with the previous version of the code. --- automol/graph/base/_02algo.py | 39 +++++++++++++++---- automol/graph/base/_11stereo.py | 8 +++- automol/tests/test_graph_ts.py | 68 +++++++++++++++++++++++++++++++-- 3 files changed, 103 insertions(+), 12 deletions(-) diff --git a/automol/graph/base/_02algo.py b/automol/graph/base/_02algo.py index 57c37148..3513281a 100644 --- a/automol/graph/base/_02algo.py +++ b/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): diff --git a/automol/graph/base/_11stereo.py b/automol/graph/base/_11stereo.py index 13b3836d..93e5ecde 100644 --- a/automol/graph/base/_11stereo.py +++ b/automol/graph/base/_11stereo.py @@ -9,6 +9,7 @@ import more_itertools as mit import numpy + from phydat import phycon from ... import util @@ -646,7 +647,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) diff --git a/automol/tests/test_graph_ts.py b/automol/tests/test_graph_ts.py index 35b4f2af..3bbb3034 100644 --- a/automol/tests/test_graph_ts.py +++ b/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", [ @@ -1527,6 +1576,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 +1604,7 @@ 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) From d86e8727ba85deac74434adf40fd8faa8a96ab3a Mon Sep 17 00:00:00 2001 From: avcopan Date: Thu, 12 Dec 2024 12:23:20 -0600 Subject: [PATCH 2/3] Fix: Never include bad ring migrations H-migrations can't traverse to the other side of a ring, so these should be filtered out in all cases. --- automol/graph/base/_05stereo.py | 12 ++++++++---- automol/graph/base/_11stereo.py | 18 +++++++++++++----- automol/tests/test_graph_ts.py | 4 +++- 3 files changed, 24 insertions(+), 10 deletions(-) diff --git a/automol/graph/base/_05stereo.py b/automol/graph/base/_05stereo.py index 6de2f47e..16d4a160 100644 --- a/automol/graph/base/_05stereo.py +++ b/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 diff --git a/automol/graph/base/_11stereo.py b/automol/graph/base/_11stereo.py index 93e5ecde..0e837fd7 100644 --- a/automol/graph/base/_11stereo.py +++ b/automol/graph/base/_11stereo.py @@ -99,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: @@ -197,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: diff --git a/automol/tests/test_graph_ts.py b/automol/tests/test_graph_ts.py index 3bbb3034..baa40cec 100644 --- a/automol/tests/test_graph_ts.py +++ b/automol/tests/test_graph_ts.py @@ -1063,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): @@ -1607,4 +1608,5 @@ def test__stereo_corrected_geometry(tsg0, ts_geo0): # 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__stereo_corrected_geometry(C5H7_TSG, C5H7_TS_GEO) + test__ts__expand_reaction_stereo("C5H7", C5H7_TSG, [2]) From 16683a4fa52cf6d30ed9b84201f0cf93ffcde33b Mon Sep 17 00:00:00 2001 From: avcopan Date: Thu, 12 Dec 2024 13:23:00 -0600 Subject: [PATCH 3/3] Fix: Stop including strained structures by default This is causing no end of headaches for me and leads to many issues. --- automol/reac/_5conv.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/automol/reac/_5conv.py b/automol/reac/_5conv.py index 0ea14e9b..af174a71 100644 --- a/automol/reac/_5conv.py +++ b/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.