Merge pull request #566 from avcopan/dev

New: Vibrational analysis with projection of translations and rotations

automol/geom/__init__.py

@@ -50,11 +50,15 @@
 from .base._0core import total_mass
 from .base._0core import center_of_mass
 from .base._0core import mass_centered
+from .base._0core import aligned_to_principal_axes
 from .base._0core import reduced_mass
 from .base._0core import inertia_tensor
 from .base._0core import principal_axes
 from .base._0core import moments_of_inertia
 from .base._0core import rotational_constants
+from .base._0core import rotational_analysis
+from .base._0core import translational_normal_modes
+from .base._0core import rotational_normal_modes
 # # geometric measurements
 from .base._0core import distance
 from .base._0core import central_angle
@@ -196,11 +200,15 @@
     'total_mass',
     'center_of_mass',
     'mass_centered',
+    'aligned_to_principal_axes',
     'reduced_mass',
     'inertia_tensor',
     'principal_axes',
     'moments_of_inertia',
     'rotational_constants',
+    'rotational_analysis',
+    'translational_normal_modes',
+    'rotational_normal_modes',
     # # geometric measurements
     'distance',
     'central_angle',

automol/geom/base/_0core.py

@@ -173,7 +173,7 @@ def string(geo, angstrom: bool = True, mode: Optional[ArrayLike] = None):
 
     # If requested, include the mode in the string
     if mode is not None:
-        assert len(lines) == natms, f"Invalid mode for geometry:\n{mode}\n{geo}"
+        mode = numpy.reshape(mode, (natms, 3))
         mode_strs = [f"{m[0]:10.6f} {m[1]:10.6f} {m[2]:10.6f}" for m in mode]
         lines = ["  ".join([s, x]) for s, x in zip(lines, mode_strs)]
 
@@ -515,6 +515,15 @@ def masses(geo, amu=True):
     return amas
 
 
+def mass_weight_vector(geo) -> numpy.ndarray:
+    """Get the mass-weighting vector for a geometry.
+
+    :param geo: The geometry
+    :return: The mass-weighting vector
+    """
+    return numpy.sqrt(numpy.repeat(masses(geo), 3))
+
+
 def total_mass(geo, amu=True):
     """Calculate the total mass of a molecular geometry.
 
@@ -598,49 +607,129 @@ def inertia_tensor(geo, amu=True):
     return ine
 
 
-def principal_axes(geo, amu=True):
-    """Determine the principal axes of rotation for a molecular geometry.
+def moments_of_inertia(geo, amu: bool = True, drop_null: bool = False):
+    """Calculate the moments of inertia along the xyz axes
+    (these not sorted in to A,B,C).
 
     :param geo: molecular geometry
     :type geo: automol geometry data structure
     :param amu: parameter to control electron mass -> amu conversion
     :type amu: bool
     :rtype: tuple(tuple(float))
     """
-    _, paxs = rotational_analysis(geo, amu=amu)
-    return paxs
+    moms, _ = rotational_analysis(geo, amu=amu, drop_null=drop_null)
+    return moms
 
 
-def moments_of_inertia(geo, amu=True):
-    """Calculate the moments of inertia along the xyz axes
-    (these not sorted in to A,B,C).
+def principal_axes(geo, drop_null: bool = False):
+    """Determine the principal axes of rotation for a molecular geometry.
 
     :param geo: molecular geometry
-    :type geo: automol geometry data structure
     :param amu: parameter to control electron mass -> amu conversion
-    :type amu: bool
     :rtype: tuple(tuple(float))
     """
-    moms, _ = rotational_analysis(geo, amu=amu)
-    return moms
+    _, rot_axes = rotational_analysis(geo, drop_null=drop_null)
+    return rot_axes
+
+
+def aligned_to_principal_axes(geo):
+    """Project the geometry onto principal axes of rotation, to standardize its
+    orientation.
+
+    :param geo: The geometry
+    :return: The geometry with standard orientation
+    """
+    return transform_by_matrix(mass_centered(geo), principal_axes(geo), trans=False)
 
 
 def rotational_analysis(
-    geo, amu: bool = True
+    geo, drop_null: bool = False, amu: bool = True
 ) -> tuple[tuple[float, ...], numpy.ndarray]:
     """Do a rotational analysis, generating the moments of inertia and principal axes.
 
     :param geo: A molecular geometry
+    :param drop_null: Drop null rotations from the analysis? (1 if linear, 3 if atom)
     :param amu: Use atomic mass units instead of electron mass unit?, defaults to True
     :return: The eigenvalues and eigenvetors of the inertia tensor
     """
     ine = inertia_tensor(geo, amu=amu)
     eig_vals, eig_vecs = numpy.linalg.eigh(ine)
+    # Sort the eigenvectors and values
     eig_vals = tuple(eig_vals)
-    return eig_vals, eig_vecs
+    # Drop null rotations, if requested
+    ndrop = 0 if not drop_null else 3 if is_atom(geo) else 1 if is_linear(geo) else 0
+    assert numpy.allclose(eig_vals[:ndrop], 0.0), f"ndrop={ndrop} eig_vals={eig_vals}"
+    return eig_vals[ndrop:], eig_vecs[:, ndrop:]
+
+
+def translational_normal_modes(geo, mass_weight: bool = True) -> numpy.ndarray:
+    """Calculate the rotational normal modes.
+
+    :param geo: The geometry
+    :param mass_weight: Return mass-weighted normal modes?
+    :return: The translational normal modes, as a 3N x 3 matrix
+    """
+    natms = count(geo)
+    trans_coos = numpy.tile(numpy.eye(3), (natms, 1))
+    if mass_weight:
+        trans_coos *= mass_weight_vector(geo)[:, numpy.newaxis]
+    return trans_coos
+
+
+def rotational_normal_modes(geo, mass_weight: bool = True) -> numpy.ndarray:
+    """Calculate the rotational normal modes.
+
+    For each atom, the direction of rotation around a given axis is given as the cross
+    product of its position with the axis.
+
+    :param geo: The geometry
+    :param mass_weight: Return mass-weighted normal modes?
+    :return: The rotational normal modes, as a 3N x (3 or 2) matrix
+    """
+    _, rot_axes = rotational_analysis(geo, drop_null=True)
+    xyzs = coordinates(mass_centered(geo))
+    rot_coos = []
+    for rot_axis in rot_axes.T:
+        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
+
+
+def external_normal_modes(geo, mass_weight: bool = True) -> numpy.ndarray:
+    """Calculate the external (translational and rotational) normal modes.
+
+    :param geo: The geometry
+    :param mass_weight: Return mass-weighted normal modes?
+    :return: The external normal modes, as a 3N x (6 or 5) matrix
+    """
+    trans_coos = translational_normal_modes(geo, mass_weight=mass_weight)
+    rot_coos = rotational_normal_modes(geo, mass_weight=mass_weight)
+    return numpy.hstack([trans_coos, rot_coos])
+
+
+def projection_onto_internal_modes(geo) -> numpy.ndarray:
+    """Get the matrix for projecting onto mass-weighted internal 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 mass_weight: Return a mass-weighted projection?
+    :return: The projection onto internal normal modes, as a 3N x (3N - 6 or 5) matrix
+    """
+    ext_coos = external_normal_modes(geo, mass_weight=True)
+    ext_dim = numpy.shape(ext_coos)[-1]
+    coo_basis, *_ = numpy.linalg.svd(ext_coos, full_matrices=True)
+    int_proj = coo_basis[:, ext_dim:]
+    return int_proj
 
 
-def rotational_constants(geo, amu=True):
+def rotational_constants(geo, amu=True, drop_null: bool = False):
     """Calculate the rotational constants.
     (these not sorted in to A,B,C).
 
@@ -650,7 +739,7 @@ def rotational_constants(geo, amu=True):
     :type amu: bool
     :rtype: tuple(float)
     """
-    moms = moments_of_inertia(geo, amu=amu)
+    moms = moments_of_inertia(geo, amu=amu, drop_null=drop_null)
     cons = numpy.divide(1.0, moms) / 4.0 / numpy.pi / phycon.SOL
     cons = tuple(cons)
     return cons
@@ -1149,20 +1238,23 @@ def transform(geo, func, idxs=None):
     return from_data(symbs, xyzs)
 
 
-def transform_by_matrix(geo, mat):
+def transform_by_matrix(geo, mat, trans: bool = True):
     """Transform the coordinates of a molecular geometry by multiplying
-    it by some input transfomration matrix.
+    it by some input transformation matrix.
 
     :param geo: molecular geometry
     :type geo: automol molecular geometry data structure
     :param mat: transformation matrix
     :type mat: tuple(tuple(float))
+    :param trans: Whether to transpose the matrix for multiplication as r . M^T
     :rtype: automol moleculer geometry data structure
     """
+    if trans:
+        mat = numpy.transpose(mat)
 
     symbs = symbols(geo)
     xyzs = coordinates(geo)
-    xyzs = numpy.dot(xyzs, numpy.transpose(mat))
+    xyzs = numpy.dot(xyzs, mat)
 
     return from_data(symbs, xyzs)
 

automol/geom/base/__init__.py

@@ -42,11 +42,15 @@
 from ._0core import total_mass
 from ._0core import center_of_mass
 from ._0core import mass_centered
+from ._0core import aligned_to_principal_axes
 from ._0core import reduced_mass
 from ._0core import inertia_tensor
 from ._0core import principal_axes
 from ._0core import moments_of_inertia
 from ._0core import rotational_constants
+from ._0core import rotational_analysis
+from ._0core import translational_normal_modes
+from ._0core import rotational_normal_modes
 # # geometric measurements
 from ._0core import distance
 from ._0core import central_angle
@@ -130,11 +134,15 @@
     'total_mass',
     'center_of_mass',
     'mass_centered',
+    'aligned_to_principal_axes',
     'reduced_mass',
     'inertia_tensor',
     'principal_axes',
     'moments_of_inertia',
     'rotational_constants',
+    'rotational_analysis',
+    'translational_normal_modes',
+    'rotational_normal_modes',
     # # geometric measurements
     'distance',
     'central_angle',

automol/geom/base/_vib.py

@@ -5,9 +5,7 @@
 import numpy
 from qcelemental import constants as qcc
 
-from phydat import ptab
-
-from ._0core import symbols
+from ._0core import masses, projection_onto_internal_modes
 
 MatrixLike = Sequence[Sequence[float]] | numpy.ndarray
 
@@ -22,16 +20,22 @@ def vibrational_analysis(
     :param freq: Return frequencies (cm^-1), instead of force constants (a.u.)?
     :return: The vibrational frequencies (or force constants) and normal modes
     """
-    # 1. build mass-weighted Hessian matrix
-    symbs = symbols(geo)
-    mw_vec = numpy.repeat(list(map(ptab.to_mass, symbs)), 3) ** -0.5
-    hess_mw = mw_vec[:, numpy.newaxis] * hess * mw_vec[numpy.newaxis, :]
+    # 1. Mass-weight the Hessian matrix
+    mw_vec = numpy.sqrt(numpy.repeat(masses(geo), 3))
+    hess_mw = hess / mw_vec[:, numpy.newaxis] / mw_vec[numpy.newaxis, :]
+
+    # 2. Project onto the space of internal motions
+    #       K = Qmwt Hmw Qmw = (PI)t Hmw (PI) = It Hint I
+    int_proj = projection_onto_internal_modes(geo)
+    hess_int = int_proj.T @ hess_mw @ int_proj
 
     # 2. compute eigenvalues and eigenvectors of the mass-weighted Hessian matrix
-    eig_vals, eig_vecs = numpy.linalg.eigh(hess_mw)
+    eig_vals, eig_vecs = numpy.linalg.eigh(hess_int)
 
     # 3. un-mass-weight the normal coordinates
-    norm_coos = mw_vec[:, numpy.newaxis] * eig_vecs
+    #       Qmw = PI (see above)    Q = Qmw / mw_vec
+    norm_coos_mw = numpy.dot(int_proj, eig_vecs) / mw_vec[:, numpy.newaxis]
+    norm_coos = norm_coos_mw / mw_vec[:, numpy.newaxis]
     if not freq:
         return eig_vals, norm_coos
 

automol/tests/test_geom.py

@@ -2,6 +2,7 @@
 """
 
 # import pytest
+import io
 from pathlib import Path
 
 import numpy
@@ -785,14 +786,32 @@ def test__repulsion_energy():
 
 def test__vibrational_analysis(data_directory_path):
     """Test automol.geom.vibrational_analysis."""
+    # Linear diatomic
+    geo = (("O", (0.0, 0.0, 0.204769)), ("H", (0.0, 0.0, -1.638153)))
+    hess_io = io.StringIO(
+        """
+            -0.00008   0.        0.        0.00008   0.        0.
+            0.        0.000021  0.        0.       -0.000021  0.
+            0.        0.        0.511701  0.        0.       -0.511701
+            0.00008   0.        0.       -0.00008   0.        0.
+            0.       -0.000021  0.        0.        0.000021  0.
+            0.        0.       -0.511701  0.        0.        0.511701
+        """
+    )
+    hess = numpy.loadtxt(hess_io)
+    ref_freqs = (3776.5,)
+    freqs, _ = geom.vibrational_analysis(geo, hess)
+    print(freqs)
+    assert len(freqs) == len(ref_freqs), f"{freqs} != {ref_freqs}"
+    assert numpy.allclose(freqs, ref_freqs, atol=1e-1), f"{freqs} != {ref_freqs}"
+
+    # Non-linear polyatomic (TS)
     geo = geom.from_xyz_string((data_directory_path / "c4h7_h2_geom.xyz").read_text())
     hess = numpy.loadtxt(data_directory_path / "c4h7_h2_hess.txt")
-    freqs = numpy.loadtxt(data_directory_path / "c4h7_h2_freqs.txt")
-    print(geo)
-    print(hess)
-    print(freqs)
+    ref_freqs = numpy.loadtxt(data_directory_path / "c4h7_h2_freqs.txt")
     freqs, _ = geom.vibrational_analysis(geo, hess)
     print(freqs)
+    assert numpy.allclose(freqs, ref_freqs, atol=1e-1), f"{freqs} !=\n{ref_freqs}"
 
 
 if __name__ == "__main__":