From 4fb70cdae9c44ac942a622206d2683403de5f037 Mon Sep 17 00:00:00 2001 From: Andreas Copan Date: Mon, 10 Jun 2024 17:00:10 -0400 Subject: [PATCH] WIP Implementing generic scan functions --- automol/geom/__init__.py | 2 ++ automol/geom/base/_0core.py | 22 ++++++++++++++++++++++ automol/geom/base/__init__.py | 2 ++ automol/reac/_6scan.py | 27 ++++++++++++++++++++++----- 4 files changed, 48 insertions(+), 5 deletions(-) diff --git a/automol/geom/__init__.py b/automol/geom/__init__.py index c05f60c0..08b5be38 100644 --- a/automol/geom/__init__.py +++ b/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', diff --git a/automol/geom/base/_0core.py b/automol/geom/base/_0core.py index 129e7102..3c128818 100644 --- a/automol/geom/base/_0core.py +++ b/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 ): diff --git a/automol/geom/base/__init__.py b/automol/geom/base/__init__.py index 98b4c5d8..29c07b21 100644 --- a/automol/geom/base/__init__.py +++ b/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', diff --git a/automol/reac/_6scan.py b/automol/reac/_6scan.py index 2b9fe353..5b58e628 100644 --- a/automol/reac/_6scan.py +++ b/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 ()