From 2fbca10e7387682ea96277d21cff72777ce73ad8 Mon Sep 17 00:00:00 2001 From: avcopan Date: Fri, 30 Aug 2024 15:58:34 -0500 Subject: [PATCH] New: Projection of torsions / internal rotations from Hessian The resulting frequencies look reasonable to me, at a glance. --- automol/geom/_2vib.py | 202 +++++++++++++++++++++---- automol/geom/__init__.py | 2 + automol/geom/base/_0core.py | 3 +- automol/geom/base/__init__.py | 2 + automol/graph/base/_06heur.py | 16 +- automol/tests/data/h2o2_freqs_proj.txt | 6 + automol/tests/test_geom.py | 19 ++- 7 files changed, 202 insertions(+), 48 deletions(-) create mode 100644 automol/tests/data/h2o2_freqs_proj.txt diff --git a/automol/geom/_2vib.py b/automol/geom/_2vib.py index 24355fc1..f27c9fc9 100644 --- a/automol/geom/_2vib.py +++ b/automol/geom/_2vib.py @@ -1,44 +1,36 @@ """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. @@ -46,6 +38,11 @@ def vibrational_analysis( :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 diff --git a/automol/geom/__init__.py b/automol/geom/__init__.py index 0b875d14..bccc0049 100644 --- a/automol/geom/__init__.py +++ b/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', diff --git a/automol/geom/base/_0core.py b/automol/geom/base/_0core.py index bbdb76fc..82b65680 100644 --- a/automol/geom/base/_0core.py +++ b/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 diff --git a/automol/geom/base/__init__.py b/automol/geom/base/__init__.py index 38fccdae..be488884 100644 --- a/automol/geom/base/__init__.py +++ b/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', diff --git a/automol/graph/base/_06heur.py b/automol/graph/base/_06heur.py index 2bffed83..5199fb49 100644 --- a/automol/graph/base/_06heur.py +++ b/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}]] """ diff --git a/automol/tests/data/h2o2_freqs_proj.txt b/automol/tests/data/h2o2_freqs_proj.txt new file mode 100644 index 00000000..9186f7a5 --- /dev/null +++ b/automol/tests/data/h2o2_freqs_proj.txt @@ -0,0 +1,6 @@ + 1027.1 + 1350.1 + 1487.1 + 3851.8 + 3854.1 + diff --git a/automol/tests/test_geom.py b/automol/tests/test_geom.py index 044ce8f4..31f0e5be 100644 --- a/automol/tests/test_geom.py +++ b/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)