From f118b3b8e336fa25ff52477b570a9a5be5eb5193 Mon Sep 17 00:00:00 2001 From: avcopan Date: Fri, 30 Aug 2024 12:26:35 -0500 Subject: [PATCH] New: Vibrational analysis works for diatomics as well Also, defines rotational modes to match the current frame of the molecule, rather than the Eckart frame. If all 3 rotations are being projected out, this isn't an issue, but if we ant to choose only 2, the frames need to match. --- automol/geom/base/_0core.py | 52 ++++++++++++++++++++----------------- automol/tests/test_geom.py | 22 ++++++++++++++++ 2 files changed, 50 insertions(+), 24 deletions(-) diff --git a/automol/geom/base/_0core.py b/automol/geom/base/_0core.py index d26be69b..4003d0ac 100644 --- a/automol/geom/base/_0core.py +++ b/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)] @@ -607,7 +607,7 @@ def inertia_tensor(geo, amu=True): return ine -def moments_of_inertia(geo, amu=True): +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). @@ -617,21 +617,19 @@ def moments_of_inertia(geo, amu=True): :type amu: bool :rtype: tuple(tuple(float)) """ - moms, *_ = rotational_analysis(geo, amu=amu) + moms, _ = rotational_analysis(geo, amu=amu, drop_null=drop_null) return moms -def principal_axes(geo): +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)) """ - _, paxs, *_ = rotational_analysis(geo) - return paxs + _, rot_axes = rotational_analysis(geo, drop_null=drop_null) + return rot_axes def aligned_to_principal_axes(geo): @@ -641,28 +639,27 @@ def aligned_to_principal_axes(geo): :param geo: The geometry :return: The geometry with standard orientation """ - *_, geo_aligned = rotational_analysis(geo) - return geo_aligned + return transform_by_matrix(mass_centered(geo), principal_axes(geo), trans=False) def rotational_analysis( - geo, amu: bool = True -) -> tuple[tuple[float, ...], numpy.ndarray, Any]: + 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, along with a - geometry aligned to the principal axes + :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 - idxs = numpy.argsort(eig_vals) - eig_vals = tuple(eig_vals[idxs]) - eig_vecs = eig_vecs[:, idxs] - geo_aligned = transform_by_matrix(mass_centered(geo), eig_vecs, trans=False) - return eig_vals, eig_vecs, geo_aligned + eig_vals = tuple(eig_vals) + # 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: @@ -682,13 +679,20 @@ def translational_normal_modes(geo, mass_weight: bool = True) -> numpy.ndarray: 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, geo_aligned = rotational_analysis(geo) - xyzs = coordinates(geo_aligned) - rot_coos = numpy.vstack([numpy.cross(xyz, rot_axes) for xyz in xyzs]) + _, 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 @@ -725,7 +729,7 @@ def projection_onto_internal_modes(geo) -> numpy.ndarray: 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). @@ -735,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 diff --git a/automol/tests/test_geom.py b/automol/tests/test_geom.py index c454ac1c..9a0f482b 100644 --- a/automol/tests/test_geom.py +++ b/automol/tests/test_geom.py @@ -2,6 +2,7 @@ """ # import pytest +import io from pathlib import Path import numpy @@ -785,10 +786,31 @@ 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") 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}"