Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

New: Complete vibrational analysis for non-linear polyatomics #563

Closed
wants to merge 2 commits into from
Closed
Show file tree
Hide file tree
Changes from 1 commit
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
New: Complete vibrational analysis for non-linear polyatomics
Next I need to make sure it handles diatomics as well...
  • Loading branch information
avcopan committed Aug 30, 2024
commit d34883de3ac84edd4590ab0b115968bd73f34190
32 changes: 25 additions & 7 deletions automol/geom/base/_0core.py
Original file line number Diff line number Diff line change
Expand Up @@ -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?
Expand All @@ -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?
Expand All @@ -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?
Expand All @@ -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).
Expand Down
23 changes: 11 additions & 12 deletions automol/geom/base/_vib.py
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -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

Expand Down
13 changes: 2 additions & 11 deletions automol/tests/test_geom.py
Original file line number Diff line number Diff line change
Expand Up @@ -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__":
Expand Down