Skip to content

Place the retro-ene migrating H on the recipe acceptor #1036

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
calvinp0 wants to merge 1 commit into
base: main
Choose a base branch
from
Open

Place the retro-ene migrating H on the recipe acceptor #1036

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
61 changes: 22 additions & 39 deletions arc/job/adapters/ts/linear_test.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -65,6 +65,7 @@
)
from arc.job.adapters.ts.linear_utils.postprocess import (FAMILY_POSTPROCESSORS,
FAMILY_VALIDATORS,
PAULING_DELTA,
has_inward_migrating_group_h,
postprocess_generic,
postprocess_ts_guess,
Expand All @@ -76,6 +77,7 @@
from arc.reaction import ARCReaction
from arc.species.converter import compare_zmats, order_mol_by_atom_map, order_xyz_by_atom_map, str_to_xyz, xyz_to_str, zmat_from_xyz
from arc.species.species import ARCSpecies, colliding_atoms
from arc.species.vectors import calculate_angle


def _angle_deg(p1: np.ndarray, vertex: np.ndarray, p2: np.ndarray) -> float:
Expand Down Expand Up @@ -4372,51 +4374,32 @@ def test_interpolate_retroene(self):
p2 = ARCSpecies(label='P2', smiles='CC(=O)O', xyz=p2_xyz)
rxn = ARCReaction(r_species=[r], p_species=[p1, p2])
ts_xyzs = interpolate(rxn)
found_ring = False
d_ch = get_single_bond_length('C', 'H') + PAULING_DELTA
d_oh = get_single_bond_length('O', 'H') + PAULING_DELTA
ring_guess = None
for ts_xyz in ts_xyzs:
coords = np.array(ts_xyz['coords'], dtype=float)
d_o3c4 = float(np.linalg.norm(coords[3] - coords[4]))
d_o2c4 = float(np.linalg.norm(coords[2] - coords[4]))
d_c5h13 = float(np.linalg.norm(coords[5] - coords[13]))
d_h13o3 = float(np.linalg.norm(coords[13] - coords[3]))
# Check 6-membered ring characteristics:
# O3-C4 stretched (breaking)
if d_o3c4 < 1.8:
if float(np.linalg.norm(coords[3] - coords[4])) < 1.8:
continue
# O2-C4 contracted (ring contact, both O's participate)
if d_o2c4 > 2.5:
if abs(float(np.linalg.norm(coords[5] - coords[13])) - d_ch) > 0.1:
continue
# H13 migrating from C5 (near Pauling TS distance)
if d_c5h13 > 1.7:
if abs(float(np.linalg.norm(coords[2] - coords[13])) - d_oh) > 0.1:
continue
# H13 approaching O3
if d_h13o3 > 2.5:
continue
found_ring = True
ring_guess = ts_xyz
break
self.assertTrue(found_ring, msg='No TS guess has 6-membered ring retroene characteristics: O3-C4 stretched '
'(breaking), O2-C4 contracted (ring), H13 migrating from C5 toward O3')
expected_ts = """C 3.35667786 -0.45750645 0.53734155
C 2.24637997 0.53978750 0.40948895
O 1.25975689 0.57306185 1.13089404
O 1.90798055 1.60139464 -0.37185519
C -0.25465664 2.12528094 0.76768134
C -0.03131089 3.04259438 -0.26607527
C -1.40160640 3.23262677 0.38080580
C 0.80777963 4.31367617 -0.11759400
H 3.43291583 -1.04518301 -0.37749397
H 4.29490197 0.05808884 0.74227330
H 3.14278867 -1.13176712 1.36663279
H -0.88551430 1.24137878 0.67368347
H -0.03778772 2.29994313 1.82151292
H 0.87654037 2.16543158 -1.09462311
H -1.31181604 3.41929928 1.45094180
H -1.92716397 4.07668244 -0.06580283
H -2.02058432 2.34596864 0.24367923
H 1.76588869 4.21544196 -0.62796354
H 0.29073745 5.16969147 -0.55118965
H 1.00524898 4.53585118 0.93109285"""
self.assertTrue(any(almost_equal_coords(ts, str_to_xyz(expected_ts)) for ts in ts_xyzs))
self.assertIsNotNone(ring_guess,
msg='No TS guess has 6-membered ring retroene characteristics: O3-C4 stretched '
'(breaking), H13 migrating from C5 toward the carbonyl O2')
ring_coords = np.array(ring_guess['coords'], dtype=float)
d_o3h13 = float(np.linalg.norm(ring_coords[3] - ring_coords[13]))
self.assertGreater(d_o3h13, 2.0,
msg=f'H13 was placed on the ester O3 instead of the carbonyl O2: {d_o3h13:.3f}')
d_o2c4 = float(np.linalg.norm(ring_coords[2] - ring_coords[4]))
self.assertGreater(d_o2c4, 2.4,
msg=f'a spurious O2-C4 contact was fabricated: {d_o2c4:.3f}')
angle = calculate_angle(coords=ring_guess, atoms=[5, 13, 2])
self.assertAlmostEqual(angle, 126.0, delta=5.0,
msg=f'C5-H13...O2 bridge angle is {angle:.1f}, expected ~126')

def test_interpolate_singlet_carbene_intra_disproportionation(self):
"""Test the interpolate_isomerization() function for Singlet_Carbene_Intra_Disproportionation: C=C1C=C[C]C1 <=> C=C1C=CC=C1"""
Expand Down
121 changes: 64 additions & 57 deletions arc/job/adapters/ts/linear_utils/families.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -42,6 +42,7 @@
xyz_with_coords,
)
from arc.job.adapters.ts.linear_utils.isomerization import ring_closure_xyz
from arc.job.adapters.ts.linear_utils.local_geometry import clean_migrating_h
from arc.job.adapters.ts.linear_utils.postprocess import PAULING_DELTA
from arc.species.species import colliding_atoms

Expand All @@ -64,6 +65,10 @@

_OH_MIGRATION_EARLY_TS_CO_DISTANCE: float = 2.08

# Retro-ene bridging-H interior angle (degrees); the six-ring strain minimum
# for this builder's seed distances.
_RETROENE_H_BRIDGE_ANGLE: float = 126.0


# ---------------------------------------------------------------------------
# Geometry helpers (thin wrappers over geom_utils where possible)
Expand Down Expand Up @@ -1215,20 +1220,36 @@ def build_retroene_ts(r_xyz: dict,
"""
Build a 6-membered-ring TS for retro-ene fragmentation.

In a Retroene reaction, an ester fragments through a concerted 6-membered ring TS::

H(mig)
/ \\
C(donor) O(ester)
| |
C(leaving) C(carbonyl)
\\ /
O(C=O)

The migrating H leaves the donor C and arrives at the ester O,
while the C=O oxygen develops a contact with the leaving C.

Returns ``None`` if the motif cannot be identified.
The RMG ``Retroene`` recipe labels the ring ``*1`` … ``*6``. Five of its
six bonds are already present in the reactant; the ring is closed by the
forming ``*1``···``*6`` bond::

*6 (migrating H)
/ \\
*5 (donor) *1 (acceptor)
| |
*4 (leaving) *2
\\ /
*3

``*5``-``*6`` and ``*3``-``*4`` break, ``*1``-``*6`` forms, ``*1``-``*2``
loses a bond order and ``*2``-``*3`` / ``*4``-``*5`` each gain one. For an
ester ``*1`` is the carbonyl O, ``*2`` the carbonyl C and ``*3`` the ester
O; for an all-carbon retro-ene ``*1`` is a carbon.

The acceptor ``*1`` is read from ``forming_bonds``; the donor ``*5``,
the migrating H ``*6``, the leaving atom ``*4`` and its sigma partner
``*3`` are read from ``breaking_bonds``. ``*4`` is first repositioned by
two-sphere triangulation so the ``*3``-``*4`` sigma bond stretches while
the ``*4``-``*5`` pi bond contracts, the ring is then folded until the
donor-acceptor separation matches a relaxed hydrogen bridge, and the
migrating H is finally triangulated between donor and acceptor at
element-specific Pauling-elongated distances, on the side of the
donor-acceptor axis opposite the rest of the ring.

Returns ``None`` if the motif cannot be identified, if ``*1`` is not two
bonds from ``*3``, or if the fold leaves the donor and acceptor too far
apart to place the migrating H between them.
"""
bb = list(breaking_bonds or [])
fb = list(forming_bonds or [])
Expand All @@ -1255,53 +1276,39 @@ def build_retroene_ts(r_xyz: dict,
if mig_h is None or other_bb_idx is None:
return None

acceptor = fb[0][0] if fb[0][1] == mig_h else fb[0][1]
if acceptor == mig_h or not 0 <= acceptor < len(symbols) or symbols[acceptor] == 'H':
return None

# The other breaking bond is the sigma bond (e.g. O3-C4).
sigma_bb = bb[other_bb_idx]

leaving_c = None
ester_o = None
sigma_partner = None
for candidate in sigma_bb:
for nbr in r_mol.atoms[candidate].bonds:
if atom_to_idx[nbr] == donor:
leaving_c = candidate
ester_o = sigma_bb[0] if sigma_bb[1] == candidate else sigma_bb[1]
sigma_partner = sigma_bb[0] if sigma_bb[1] == candidate else sigma_bb[1]
break
if leaving_c is not None:
break
if leaving_c is None or ester_o is None:
return None

# carbonyl_c = non-H, non-leaving-C neighbour of ester_o.
carbonyl_c = None
for nbr in r_mol.atoms[ester_o].bonds:
ni = atom_to_idx[nbr]
if ni != leaving_c and symbols[ni] != 'H':
carbonyl_c = ni
break
if carbonyl_c is None:
if leaving_c is None or sigma_partner is None:
return None

# carbonyl_o = double-bonded O on carbonyl_c (not ester_o).
carbonyl_o = None
for nbr, bond in r_mol.atoms[carbonyl_c].bonds.items():
ni = atom_to_idx[nbr]
if symbols[ni] == 'O' and ni != ester_o and bond.order >= 1.5:
carbonyl_o = ni
break
if carbonyl_o is None:
adj = mol_to_adjacency(r_mol)
ring_bridge = sorted(adj[acceptor] & adj[sigma_partner])
if not ring_bridge:
logger.debug(f'Retroene builder: acceptor {acceptor} is not two bonds from {sigma_partner}, '
f'so no *1-*2-*3 ring path exists; skipping.')
return None

# The acid fragment and the donor stay at reactant positions; only C4 is
# repositioned via two-sphere triangulation so d(O3,C4) stretches and
# d(C5,C4) contracts simultaneously. Then ring_closure_xyz folds O2
# toward C4, and finally H13 is placed between C5 and O3.

# Step 1: Place C4 at the two-sphere intersection of O3 and C5.
d_break = 2.5 # O3-C4 (breaking sigma)
d_pi = 1.40 # C5-C4 (forming pi bond)

new_c4 = two_sphere_intersection(
coords[ester_o], d_break, coords[donor], d_pi, coords[leaving_c])
coords[sigma_partner], d_break, coords[donor], d_pi, coords[leaving_c])
if new_c4 is None:
return None
c4_displacement = new_c4 - coords[leaving_c]
Expand All @@ -1315,28 +1322,28 @@ def build_retroene_ts(r_xyz: dict,
for k in c4_moving:
coords[k] += c4_displacement

# Step 2: Fold O2 toward C4's new position.
intermediate_xyz = xyz_with_coords(r_xyz, coords)
d_co_target = 2.2
d_donor_h = get_single_bond_length(symbols[donor], 'H') + PAULING_DELTA
d_acceptor_h = get_single_bond_length(symbols[acceptor], 'H') + PAULING_DELTA
d_donor_acceptor = float(np.sqrt(
d_donor_h ** 2 + d_acceptor_h ** 2
- 2.0 * d_donor_h * d_acceptor_h * np.cos(np.radians(_RETROENE_H_BRIDGE_ANGLE))))

rc_xyz = ring_closure_xyz(
intermediate_xyz, r_mol,
forming_bond=(carbonyl_o, leaving_c),
target_distance=d_co_target)
xyz_with_coords(r_xyz, coords), r_mol,
forming_bond=(acceptor, donor),
target_distance=d_donor_acceptor)
if rc_xyz is not None:
coords = np.array(rc_xyz['coords'], dtype=float)

# Step 3: Place migrating H between C5 and O3.
sbl_ch = get_single_bond_length('C', 'H') or 1.09
sbl_oh = get_single_bond_length('O', 'H') or 0.97
d_ch = sbl_ch + PAULING_DELTA # ~1.51 Å
d_oh = sbl_oh + PAULING_DELTA # ~1.39 Å
new_h_pos = two_sphere_intersection(
coords[donor], d_ch, coords[ester_o], d_oh, coords[mig_h])
if new_h_pos is None:
gap = float(np.linalg.norm(coords[acceptor] - coords[donor]))
if not abs(d_donor_h - d_acceptor_h) <= gap <= d_donor_h + d_acceptor_h:
logger.debug(f'Retroene builder: the ring did not close -- donor {donor} and acceptor {acceptor} '
f'are {gap:.3f} A apart, outside the triangulable window; skipping.')
return None
coords[mig_h] = new_h_pos

ts_xyz = xyz_with_coords(r_xyz, coords)
ring_interior = coords[[ring_bridge[0], sigma_partner, leaving_c]].mean(axis=0)
ts_xyz = clean_migrating_h(xyz_with_coords(r_xyz, coords), donor, acceptor, mig_h,
ref_pos=coords[donor] + coords[acceptor] - ring_interior)
if colliding_atoms(ts_xyz):
return None
return ts_xyz
Loading
Loading