New: Complete vibrational analysis for non-linear polyatomics

Next I need to make sure it handles diatomics as well...

automol/geom/base/_0core.py

@@ -665,8 +665,8 @@ def rotational_analysis(
     return eig_vals, eig_vecs, geo_aligned
 
 
-def translational_normal_modes(geo, mass_weight: bool=False):
-    """Calculate the rotational normal modes of a geometry.
+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?
@@ -679,8 +679,8 @@ def translational_normal_modes(geo, mass_weight: bool=False):
     return trans_coos
 
 
-def rotational_normal_modes(geo, mass_weight: bool=False):
-    """Calculate the rotational normal modes of a geometry.
+def rotational_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?
@@ -694,9 +694,8 @@ def rotational_normal_modes(geo, mass_weight: bool=False):
     return rot_coos
 
 
-def external_normal_modes(geo, mass_weight: bool=False):
-    """Calculate the external (translational and rotational) normal modes for a
-    geometry.
+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?
@@ -707,6 +706,25 @@ def external_normal_modes(geo, mass_weight: bool=False):
     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):
     """Calculate the rotational constants.
     (these not sorted in to A,B,C).

automol/geom/base/_vib.py

@@ -5,7 +5,7 @@
 import numpy
 from qcelemental import constants as qcc
 
-from ._0core import coordinates, masses, rotational_analysis
+from ._0core import masses, projection_onto_internal_modes
 
 MatrixLike = Sequence[Sequence[float]] | numpy.ndarray
 
@@ -20,23 +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
+    # 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. compute eigenvalues and eigenvectors of the mass-weighted Hessian matrix
-    eig_vals, eig_vecs = numpy.linalg.eigh(hess_mw)
+    # 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
 
-    # 3. Project out translations and rotations
-    _, rot_axes, geo_aligned = rotational_analysis(geo)
-    eck_xyzs = numpy.array(coordinates(geo_aligned))
-    print("Eck 1", numpy.linalg.norm(eck_xyzs))
-    print(eck_xyzs)
-    print("Axes 1", numpy.linalg.norm(rot_axes))
-    print(rot_axes)
+    # 2. compute eigenvalues and eigenvectors of the mass-weighted Hessian matrix
+    eig_vals, eig_vecs = numpy.linalg.eigh(hess_int)
 
     # 3. un-mass-weight the normal coordinates
-    norm_coos = eig_vecs / mw_vec[:, numpy.newaxis]
+    #       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

@@ -787,18 +787,9 @@ def test__vibrational_analysis(data_directory_path):
     """Test automol.geom.vibrational_analysis."""
     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)
-
-    # from automol.geom.base import _vib2
-    # freqs2, _ = _vib2.normal_modes(geo, hess)
-    # print(freqs2)
-
-    geom.rotational_normal_modes(geo, mass_weight=True)
+    assert numpy.allclose(freqs, ref_freqs, atol=1e-1), f"{freqs} !=\n{ref_freqs}"
 
 
 if __name__ == "__main__":