New: Projection of torsions / internal rotations from Hessian

The resulting frequencies look reasonable to me, at a glance.

automol/geom/_2vib.py

@@ -1,51 +1,48 @@
 """Vibrational analysis."""
 
-from collections.abc import Sequence
+from collections.abc import Callable, Collection, Sequence
 
 import numpy
 from qcelemental import constants as qcc
 
-from .base import masses, rotational_normal_modes, translational_normal_modes
+from .. import util
+from ..graph import base as graph_base
+from ._1conv import graph
+from .base import (
+    coordinates,
+    count,
+    mass_weight_vector,
+    masses,
+    rotational_normal_modes,
+    translational_normal_modes,
+)
 
 MatrixLike = Sequence[Sequence[float]] | numpy.ndarray
 
 
-def normal_mode_projection(
-    geo, trans: bool = False, rot: bool = False
-) -> numpy.ndarray:
-    """Get the matrix for projecting onto a subset of normal modes.
-
-    Note: This must be applied to the *mass-weighted* normal modes!
-
-    Uses SVD to calculate an orthonormal basis for the null space of the external normal
-    modes, which is the space of internal normal modes.
-
-    :param geo: The geometry
-    :param trans: Keep translations, instead of removing them?
-    :param rot: Keep rotations, instead of removing them?
-    :return: The projection onto a subset normal modes
-    """
-    coos_lst = []
-    if not trans:
-        coos_lst.append(translational_normal_modes(geo, mass_weight=True))
-    if not rot:
-        coos_lst.append(rotational_normal_modes(geo, mass_weight=True))
-    coos = numpy.hstack(coos_lst)
-    dim = numpy.shape(coos)[-1]
-    coo_basis, *_ = numpy.linalg.svd(coos, full_matrices=True)
-    proj = coo_basis[:, dim:]
-    return proj
-
-
 def vibrational_analysis(
-    geo, hess: MatrixLike, trans: bool = False, rot: bool = False, wavenum: bool = True
+    geo,
+    hess: MatrixLike,
+    trans: bool = False,
+    rot: bool = False,
+    tors: bool = True,
+    gra: object | None = None,
+    bkeys: Sequence[Collection[int]] | None = None,
+    with_h_rotors: bool = True,
+    with_ch_rotors: bool = True,
+    wavenum: bool = True,
 ) -> tuple[tuple[float, ...], numpy.ndarray]:
     """Get vibrational modes and frequencies from the Hessian matrix.
 
     :param geo: The molecular geometry
     :param hess: The Hessian matrix
     :param trans: Keep translations, instead of removing them?
     :param rot: Keep rotations, instead of removing them?
+    :param tors: Keep torsions, instead of removing them?
+    :param gra: For torsions, a graph specifying connectivity
+    :param bkeys: For torsions, specify the rotational bonds by index
+    :param with_h_rotors: For torsions, include XH rotors?
+    :param with_ch_rotors: For torsions, include CH rotors?
     :param wavenum: Return frequencies (cm^-1), instead of force constants (a.u.)?
     :return: The vibrational frequencies (or force constants) and normal modes
     """
@@ -55,7 +52,16 @@ def vibrational_analysis(
 
     # 2. Project onto the space of internal motions
     #       K = Qmwt Hmw Qmw = (PI)t Hmw (PI) = It Hint I
-    proj = normal_mode_projection(geo, trans=trans, rot=rot)
+    proj = normal_mode_projection(
+        geo,
+        trans=trans,
+        rot=rot,
+        tors=tors,
+        gra=gra,
+        bkeys=bkeys,
+        with_h_rotors=with_h_rotors,
+        with_ch_rotors=with_ch_rotors,
+    )
     hess_proj = proj.T @ hess_mw @ proj
 
     # 2. Compute eigenvalues and eigenvectors of the mass-weighted Hessian matrix
@@ -78,3 +84,137 @@ def vibrational_analysis(
     freqs = tuple(map(float, numpy.real(freqs) - numpy.imag(freqs)))
 
     return freqs, norm_coos
+
+
+def normal_mode_projection(
+    geo,
+    trans: bool = False,
+    rot: bool = False,
+    tors: bool = True,
+    gra: object | None = None,
+    bkeys: Sequence[Collection[int]] | None = None,
+    with_h_rotors: bool = True,
+    with_ch_rotors: bool = True,
+) -> numpy.ndarray:
+    """Get the matrix for projecting onto a subset of normal modes.
+
+    Note: This must be applied to the *mass-weighted* normal modes!
+
+    Uses SVD to calculate an orthonormal basis for the null space of the external normal
+    modes, which is the space of internal normal modes.
+
+    :param geo: The geometry
+    :param trans: Keep translations, instead of removing them?
+    :param rot: Keep rotations, instead of removing them?
+    :param tors: Keep torsions, instead of removing them?
+    :param gra: For torsions, a graph specifying connectivity
+    :param bkeys: For torsions, specify the rotational bonds by index
+    :param with_h_rotors: For torsions, include XH rotors?
+    :param with_ch_rotors: For torsions, include CH rotors?
+    :return: The projection onto a subset normal modes
+    """
+    coos_lst = []
+    if not trans:
+        coos_lst.append(translational_normal_modes(geo, mass_weight=True))
+    if not rot:
+        coos_lst.append(rotational_normal_modes(geo, mass_weight=True))
+    if not tors:
+        coos_lst.append(
+            torsional_normal_modes(
+                geo,
+                gra=gra,
+                bkeys=bkeys,
+                with_h_rotors=with_h_rotors,
+                with_ch_rotors=with_ch_rotors,
+            )
+        )
+
+    # If no modes are being project, return an identity matrix
+    if not coos_lst:
+        return numpy.eye(count(geo) * 3)
+
+    coos = numpy.hstack(coos_lst)
+    dim = numpy.shape(coos)[-1]
+    coo_basis, *_ = numpy.linalg.svd(coos, full_matrices=True)
+    proj = coo_basis[:, dim:]
+    return proj
+
+
+# Torsions
+def torsional_normal_modes(
+    geo,
+    gra: object | None = None,
+    bkeys: Sequence[Collection[int]] | None = None,
+    mass_weight: bool = True,
+    with_h_rotors: bool = True,
+    with_ch_rotors: bool = True,
+) -> numpy.ndarray:
+    """Calculate the torsional normal modes at rotational bonds.
+
+    :param geo: The geometry
+    :param gra: A graph specifying the connectivity of the geometry
+    :param bkeys: Specify the rotational bonds by index (ignores `with_x_rotors` flags);
+        If `None`, they will be determined automatically
+    :param mass_weight: Return mass-weighted normal modes?
+    :param with_h_rotors: Include rotors neighbored by only hydrogens on one side?
+    :param with_ch_rotors: Include rotors with C neighbored by only hydrogens?
+    :return: The torsional normal modes, as a 3N x n_tors matrix
+    """
+    gra = graph(geo, stereo=False) if gra is None else gra
+    bkeys = (
+        graph_base.rotational_bond_keys(
+            gra,
+            with_h_rotors=with_h_rotors,
+            with_ch_rotors=with_ch_rotors,
+            with_rings_rotors=False,
+        )
+        if bkeys is None
+        else bkeys
+    )
+
+    all_idxs = list(range(count(geo)))
+    tors_coos = []
+    for bkey in bkeys:
+        axis_idxs = sorted(bkey)
+        groups = graph_base.rotational_groups(gra, *axis_idxs)
+        tors_ = torsional_motion_calculator_(geo, axis_idxs, groups)
+        tors_coo = numpy.concatenate([tors_(idx) for idx in all_idxs])
+        tors_coos.append(tors_coo)
+    tors_coos = numpy.transpose(tors_coos)
+
+    if mass_weight:
+        tors_coos *= mass_weight_vector(geo)[:, numpy.newaxis]
+    return tors_coos
+
+
+def torsional_motion_calculator_(
+    geo, axis_idxs: Sequence[int], groups: Sequence[Sequence[int]]
+) -> Callable[[int], numpy.ndarray]:
+    """Generate a torsional motion calculator.
+
+    :param geo: The geometry
+    :param axis_idxs: The indices of the rotational axis
+    :param groups: The groups of the rotational axis
+    """
+    group1, group2 = groups
+    xyz1, xyz2 = coordinates(geo, idxs=axis_idxs)
+    axis1 = util.vector.unit_norm(numpy.subtract(xyz1, xyz2))
+    axis2 = numpy.negative(axis1)
+
+    def _torsional_motion(idx: int) -> numpy.ndarray:
+        """Determine the torsional motion of an atom."""
+        # If it is one of the axis atoms, it doesn't move
+        if idx in axis_idxs:
+            return numpy.zeros(3)
+        # If it is in the first group, subtract and cross
+        (xyz,) = coordinates(geo, idxs=(idx,))
+        if idx in group1:
+            dxyz = numpy.subtract(xyz, xyz1)
+            return numpy.cross(dxyz, axis1)
+        if idx in group2:
+            dxyz = numpy.subtract(xyz, xyz2)
+            return numpy.cross(dxyz, axis2)
+
+        raise ValueError(f"{idx} not in {axis_idxs} {group1} {group2}")
+
+    return _torsional_motion

automol/geom/__init__.py

@@ -47,6 +47,7 @@
 from .base._0core import atom_indices
 from .base._0core import dummy_atom_indices
 from .base._0core import masses
+from .base._0core import mass_weight_vector
 from .base._0core import total_mass
 from .base._0core import center_of_mass
 from .base._0core import mass_centered
@@ -197,6 +198,7 @@
     'atom_indices',
     'dummy_atom_indices',
     'masses',
+    'mass_weight_vector',
     'total_mass',
     'center_of_mass',
     'mass_centered',

automol/geom/base/_0core.py

@@ -662,7 +662,7 @@ def rotational_analysis(
 
 
 def translational_normal_modes(geo, mass_weight: bool = True) -> numpy.ndarray:
-    """Calculate the rotational normal modes.
+    """Calculate the translational normal modes.
 
     :param geo: The geometry
     :param mass_weight: Return mass-weighted normal modes?
@@ -692,6 +692,7 @@ def rotational_normal_modes(geo, mass_weight: bool = True) -> numpy.ndarray:
         rot_coo = numpy.concatenate([numpy.cross(xyz, rot_axis) for xyz in xyzs])
         rot_coos.append(rot_coo)
     rot_coos = numpy.transpose(rot_coos)
+
     if mass_weight:
         rot_coos *= mass_weight_vector(geo)[:, numpy.newaxis]
     return rot_coos

automol/geom/base/__init__.py

@@ -39,6 +39,7 @@
 from ._0core import atom_indices
 from ._0core import dummy_atom_indices
 from ._0core import masses
+from ._0core import mass_weight_vector
 from ._0core import total_mass
 from ._0core import center_of_mass
 from ._0core import mass_centered
@@ -130,6 +131,7 @@
     'atom_indices',
     'dummy_atom_indices',
     'masses',
+    'mass_weight_vector',
     'total_mass',
     'center_of_mass',
     'mass_centered',

automol/graph/base/_06heur.py

@@ -176,10 +176,9 @@ def rotational_bond_keys(
     :param gra: the graph
     :param lin_keys: keys to linear atoms in the graph
     :type lin_keys: list[int]
-    :param with_h_rotors: Include H rotors?
-    :type with_h_rotors: bool
-    :param with_ch_rotors: Include CH rotors?
-    :type with_ch_rotors: bool
+    :param with_h_rotors: Include rotors neighbored by only hydrogens on one side?
+    :param with_ch_rotors: Include rotors with C neighbored by only hydrogens?
+    :param with_rings_rotors: Include rotors in rings?
     :returns: The rotational bond keys
     :rtype: frozenset[frozenset[{int, int}]]
     """
@@ -212,12 +211,9 @@ def rotational_segment_keys(
     :param gra: the graph
     :param lin_keys: keys to linear atoms in the graph
     :type lin_keys: list[int]
-    :param with_h_rotors: Include H rotors?
-    :type with_h_rotors: bool
-    :param with_ch_rotors: Include CH rotors?
-    :type with_ch_rotors: bool
-    :param with_rings_rotors: Include atoms in a ring?
-    :type with_rings_rotors: bool
+    :param with_h_rotors: Include rotors neighbored by only hydrogens on one side?
+    :param with_ch_rotors: Include rotors with C neighbored by only hydrogens?
+    :param with_rings_rotors: Include rotors in rings?
     :returns: The rotational bond keys
     :rtype: frozenset[frozenset[{int, int}]]
     """

automol/tests/data/h2o2_freqs_proj.txt

@@ -0,0 +1,6 @@
+  1027.1
+  1350.1
+  1487.1
+  3851.8
+  3854.1
+

automol/tests/test_geom.py

@@ -786,14 +786,14 @@ def test__repulsion_energy():
 
 
 @pytest.mark.parametrize(
-    "smi,fml",
+    "smi,fml,tors",
     [
-        ("[OH]", "oh"),
-        ("OO", "h2o2"),
-        ("C1=CCCC1.[H]", "c4h7_h2"),
+        ("[OH]", "oh", False),
+        ("OO", "h2o2", True),
+        ("C1=CCCC1.[H]", "c4h7_h2", False),
     ],
 )
-def test__vibrational_analysis(smi, fml):
+def test__vibrational_analysis(smi, fml, tors):
     """Test automol.geom.vibrational_analysis."""
     print(f"Testing frequency projection for {fml} (SMILES {smi})")
     geo = geom.from_xyz_string((DATA_PATH / f"{fml}_geom.xyz").read_text())
@@ -803,6 +803,13 @@ def test__vibrational_analysis(smi, fml):
     print(freqs)
     assert numpy.allclose(freqs, ref_freqs, atol=1e-1), f"{freqs} !=\n{ref_freqs}"
 
+    # If requested, test projecting out of torsions
+    if tors:
+        ref_freqs = numpy.loadtxt(DATA_PATH / f"{fml}_freqs_proj.txt")
+        freqs, _ = geom.vibrational_analysis(geo, hess, tors=False)
+        print(freqs)
+        assert numpy.allclose(freqs, ref_freqs, atol=1e-1), f"{freqs} !=\n{ref_freqs}"
+
 
 if __name__ == "__main__":
     # __align()
@@ -821,4 +828,4 @@ def test__vibrational_analysis(smi, fml):
     # test__repulsion_energy()
     # test__dist_analysis()
     # test__argunique_coulomb_spectrum()
-    test__vibrational_analysis("OO", "h2o2")
+    test__vibrational_analysis("OO", "h2o2", True)