Merge pull request #499 from avcopan/dev

WIP Implementing generic scan functions

automol/geom/__init__.py

@@ -58,6 +58,7 @@
 from automol.geom.base._0core import distance
 from automol.geom.base._0core import central_angle
 from automol.geom.base._0core import dihedral_angle
+from automol.geom.base._0core import measure
 from automol.geom.base._0core import zmatrix_row_values
 # # binary functions
 from automol.geom.base._0core import join
@@ -191,6 +192,7 @@
     'distance',
     'central_angle',
     'dihedral_angle',
+    'measure',
     'zmatrix_row_values',
     # # binary functions
     'join',

automol/geom/base/_0core.py

@@ -706,6 +706,28 @@ def dihedral_angle(geo, idx1, idx2, idx3, idx4, degree=False):
     return dih
 
 
+def measure(geo, coo, angstrom=False, degree=False):
+    """Measure a coordinate value for distance, central angle, or dihedral angle
+
+    :param geo: A molecular geometry
+    :param coo: The indices defining the coordinate
+        2 for distance, 3 for central angle, 4 for dihedral angle
+    :param angstrom: Give distances in angstrom?, defaults to False
+    :param degree: Give angles in degrees?, defaults to False
+    :return: The measured value
+    """
+    assert 2 <= len(coo) <= 4, f"Invalid coordinate: {coo}"
+
+    if len(coo) == 2:
+        return distance(geo, *coo, angstrom=angstrom)
+    if len(coo) == 3:
+        return central_angle(geo, *coo, degree=degree)
+    if len(coo) == 4:
+        return dihedral_angle(geo, *coo, degree=degree)
+
+    return None
+
+
 def zmatrix_row_values(
     geo, idx, idx1=None, idx2=None, idx3=None, angstrom=True, degree=True
 ):

automol/geom/base/__init__.py

@@ -50,6 +50,7 @@
 from automol.geom.base._0core import distance
 from automol.geom.base._0core import central_angle
 from automol.geom.base._0core import dihedral_angle
+from automol.geom.base._0core import measure
 from automol.geom.base._0core import zmatrix_row_values
 # # binary functions
 from automol.geom.base._0core import join
@@ -134,6 +135,7 @@
     'distance',
     'central_angle',
     'dihedral_angle',
+    'measure',
     'zmatrix_row_values',
     # # binary functions
     'join',

automol/reac/_6scan.py

@@ -2,6 +2,7 @@
 """
 
 import numpy
+import more_itertools as mit
 
 from automol import geom
 from automol.const import ReactionClass
@@ -10,6 +11,7 @@
 from automol.reac._1util import (
     hydrogen_migration_atom_keys,
     hydrogen_migration_might_dissociate,
+    ring_forming_scission_chain,
 )
 
 
@@ -23,6 +25,10 @@ def scan_coordinates(rxn: Reaction) -> tuple[int, ...]:
         (frm_bkey,) = ts.forming_bond_keys(ts_graph(rxn))
         scan_coo = tuple(sorted(frm_bkey))
         return (scan_coo,)
+    if class_(rxn) in (ReactionClass.BETA_SCISSION, ReactionClass.RING_FORM_SCISSION):
+        (frm_bkey,) = ts.breaking_bond_keys(ts_graph(rxn))
+        scan_coo = tuple(sorted(frm_bkey))
+        return (scan_coo,)
 
     return ()
 
@@ -39,11 +45,13 @@ def scan_values(rxn: Reaction) -> tuple[numpy.ndarray, ...]:
 
     # 1. Form a grid for the first coordinate
     dist1 = dists[0]
-    npoints1, start1, end1 = {
-        # ReactionClass.HYDROGEN_MIGRATION: (18, dist1 - 0.1, dist1 + 0.3),
-        # ReactionClass.HYDROGEN_MIGRATION: (9, dist1 - 0.35, dist1 + 0.05),
-        ReactionClass.HYDROGEN_MIGRATION: (16, dist1 + 0.05, dist1 - 0.35),
-    }.get(class_(rxn), (None, None, None))
+    npoints1 = start1 = end1 = None
+    if class_(rxn) == ReactionClass.HYDROGEN_MIGRATION:
+        npoints1, start1, end1 = (16, dist1 + 0.05, dist1 - 0.55)
+    if class_(rxn) == ReactionClass.BETA_SCISSION:
+        npoints1, start1, end1 = (14, dist1 + 0.1, dist1 + 0.8)
+    if class_(rxn) == ReactionClass.RING_FORM_SCISSION:
+        npoints1, start1, end1 = (7, dist1 + 0.1, dist1 + 0.7)
 
     grid = numpy.linspace(start1, end1, npoints1)
 
@@ -61,5 +69,14 @@ def constraint_coordinates(rxn: Reaction) -> tuple[tuple[int, ...]]:
         diss_coo1 = tuple(sorted([att_key, att_nkey]))
         diss_coo2 = hydrogen_migration_might_dissociate(rxn, att_key, att_nkey, don_key)
         return (diss_coo1,) if diss_coo2 is None else (diss_coo1, diss_coo2)
+    if class_(rxn) == ReactionClass.RING_FORM_SCISSION:
+        chain_keys = ring_forming_scission_chain(rxn)
+        ang_keys_lst = []
+        ang_keys_lst.extend(sorted(mit.windowed(chain_keys[1:], 3)))
+        ang_keys_lst.extend(sorted(mit.windowed(chain_keys, 4)))
+        # Use standard coordinate keys
+        rev_ang_keys_lst = list(map(tuple, map(reversed, ang_keys_lst)))
+        ang_keys_lst = tuple(map(min, zip(ang_keys_lst, rev_ang_keys_lst)))
+        return ang_keys_lst
 
     return ()