From 2a382589886e912fd5d30290bd666107d400f4e6 Mon Sep 17 00:00:00 2001 From: avcopan Date: Fri, 30 Aug 2024 13:23:24 -0500 Subject: [PATCH 1/5] Adds H2O2 as a test case for frequency projection --- automol/tests/data/h2o2_freqs.txt | 6 ++++++ automol/tests/data/h2o2_geom.xyz | 6 ++++++ automol/tests/data/h2o2_hess.txt | 12 ++++++++++++ automol/tests/test_geom.py | 8 ++++++++ 4 files changed, 32 insertions(+) create mode 100644 automol/tests/data/h2o2_freqs.txt create mode 100644 automol/tests/data/h2o2_geom.xyz create mode 100644 automol/tests/data/h2o2_hess.txt diff --git a/automol/tests/data/h2o2_freqs.txt b/automol/tests/data/h2o2_freqs.txt new file mode 100644 index 00000000..853fc2df --- /dev/null +++ b/automol/tests/data/h2o2_freqs.txt @@ -0,0 +1,6 @@ + 355.7 + 1029.3 + 1350.1 + 1491.2 + 3851.8 + 3854.2 diff --git a/automol/tests/data/h2o2_geom.xyz b/automol/tests/data/h2o2_geom.xyz new file mode 100644 index 00000000..827d21cb --- /dev/null +++ b/automol/tests/data/h2o2_geom.xyz @@ -0,0 +1,6 @@ + 4 + +O 0.698619 0.118502 -0.054402 +O -0.698619 -0.118502 -0.054402 +H 1.019206 -0.649819 0.435212 +H -1.019206 0.649819 0.435212 \ No newline at end of file diff --git a/automol/tests/data/h2o2_hess.txt b/automol/tests/data/h2o2_hess.txt new file mode 100644 index 00000000..8e701e0c --- /dev/null +++ b/automol/tests/data/h2o2_hess.txt @@ -0,0 +1,12 @@ +4.038979999999999793e-01 -9.803489999999999427e-02 6.458379999999999677e-02 -2.920800000000000063e-01 -2.046950000000000158e-02 2.360330000000000067e-02 -7.534230000000000094e-02 1.387009999999999910e-01 -8.501179999999999848e-02 -3.647580000000000267e-02 -2.019680000000000086e-02 -3.175299999999999824e-03 +-9.803489999999999427e-02 4.128060000000000063e-01 -2.317750000000000088e-01 -2.046950000000000158e-02 -6.724800000000000222e-02 2.662240000000000179e-03 8.404190000000000271e-02 -3.576199999999999934e-01 2.233670000000000100e-01 3.446239999999999720e-02 1.206239999999999916e-02 5.745120000000000261e-03 +6.458379999999999677e-02 -2.317750000000000088e-01 1.604299999999999893e-01 -2.360330000000000067e-02 -2.662240000000000179e-03 -2.052930000000000033e-02 -5.824010000000000298e-02 2.249269999999999881e-01 -1.418089999999999906e-01 1.725959999999999994e-02 9.509999999999999357e-03 1.907929999999999956e-03 +-2.920800000000000063e-01 -2.046950000000000158e-02 -2.360330000000000067e-02 4.038979999999999793e-01 -9.803489999999999427e-02 -6.458379999999999677e-02 -3.647580000000000267e-02 -2.019680000000000086e-02 3.175299999999999824e-03 -7.534230000000000094e-02 1.387009999999999910e-01 8.501179999999999848e-02 +-2.046950000000000158e-02 -6.724800000000000222e-02 -2.662240000000000179e-03 -9.803489999999999427e-02 4.128060000000000063e-01 2.317750000000000088e-01 3.446239999999999720e-02 1.206239999999999916e-02 -5.745120000000000261e-03 8.404190000000000271e-02 -3.576199999999999934e-01 -2.233670000000000100e-01 +2.360330000000000067e-02 2.662240000000000179e-03 -2.052930000000000033e-02 -6.458379999999999677e-02 2.317750000000000088e-01 1.604299999999999893e-01 -1.725959999999999994e-02 -9.509999999999999357e-03 1.907929999999999956e-03 5.824010000000000298e-02 -2.249269999999999881e-01 -1.418089999999999906e-01 +-7.534230000000000094e-02 8.404190000000000271e-02 -5.824010000000000298e-02 -3.647580000000000267e-02 3.446239999999999720e-02 -1.725959999999999994e-02 1.189530000000000032e-01 -1.187219999999999942e-01 7.866810000000000469e-02 -7.134820000000000051e-03 2.178969999999999909e-04 -3.168380000000000155e-03 +1.387009999999999910e-01 -3.576199999999999934e-01 2.249269999999999881e-01 -2.019680000000000086e-02 1.206239999999999916e-02 -9.509999999999999357e-03 -1.187219999999999942e-01 3.466210000000000124e-01 -2.165199999999999902e-01 2.178969999999999909e-04 -1.062639999999999895e-03 1.102649999999999941e-03 +-8.501179999999999848e-02 2.233670000000000100e-01 -1.418089999999999906e-01 3.175299999999999824e-03 -5.745120000000000261e-03 1.907929999999999956e-03 7.866810000000000469e-02 -2.165199999999999902e-01 1.388570000000000082e-01 3.168380000000000155e-03 -1.102649999999999941e-03 1.043880000000000068e-03 +-3.647580000000000267e-02 3.446239999999999720e-02 1.725959999999999994e-02 -7.534230000000000094e-02 8.404190000000000271e-02 5.824010000000000298e-02 -7.134820000000000051e-03 2.178969999999999909e-04 3.168380000000000155e-03 1.189530000000000032e-01 -1.187219999999999942e-01 -7.866810000000000469e-02 +-2.019680000000000086e-02 1.206239999999999916e-02 9.509999999999999357e-03 1.387009999999999910e-01 -3.576199999999999934e-01 -2.249269999999999881e-01 2.178969999999999909e-04 -1.062639999999999895e-03 -1.102649999999999941e-03 -1.187219999999999942e-01 3.466210000000000124e-01 2.165199999999999902e-01 +-3.175299999999999824e-03 5.745120000000000261e-03 1.907929999999999956e-03 8.501179999999999848e-02 -2.233670000000000100e-01 -1.418089999999999906e-01 -3.168380000000000155e-03 1.102649999999999941e-03 1.043880000000000068e-03 -7.866810000000000469e-02 2.165199999999999902e-01 1.388570000000000082e-01 diff --git a/automol/tests/test_geom.py b/automol/tests/test_geom.py index 9a0f482b..ca7bc80b 100644 --- a/automol/tests/test_geom.py +++ b/automol/tests/test_geom.py @@ -805,6 +805,14 @@ def test__vibrational_analysis(data_directory_path): 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 (H2O2) + geo = geom.from_xyz_string((data_directory_path / "h2o2_geom.xyz").read_text()) + hess = numpy.loadtxt(data_directory_path / "h2o2_hess.txt") + ref_freqs = numpy.loadtxt(data_directory_path / "h2o2_freqs.txt") + freqs, _ = geom.vibrational_analysis(geo, hess) + print(freqs) + assert numpy.allclose(freqs, ref_freqs, atol=1e-1), f"{freqs} !=\n{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") From 241988b754b8967dd6b87f0aba2b7cb93c088551 Mon Sep 17 00:00:00 2001 From: avcopan Date: Fri, 30 Aug 2024 13:24:44 -0500 Subject: [PATCH 2/5] Refactor: Add numbering to higher-level geom module files --- automol/geom/{_molsym.py => _0molsym.py} | 0 automol/geom/{_conv.py => _1conv.py} | 8 ++-- automol/geom/__init__.py | 60 ++++++++++++------------ automol/geom/_align.py | 2 +- automol/geom/_extra.py | 2 +- automol/geom/_ring.py | 2 +- automol/geom/base/_0core.py | 1 - 7 files changed, 37 insertions(+), 38 deletions(-) rename automol/geom/{_molsym.py => _0molsym.py} (100%) rename automol/geom/{_conv.py => _1conv.py} (99%) diff --git a/automol/geom/_molsym.py b/automol/geom/_0molsym.py similarity index 100% rename from automol/geom/_molsym.py rename to automol/geom/_0molsym.py diff --git a/automol/geom/_conv.py b/automol/geom/_1conv.py similarity index 99% rename from automol/geom/_conv.py rename to automol/geom/_1conv.py index a0cf663d..25101d98 100644 --- a/automol/geom/_conv.py +++ b/automol/geom/_1conv.py @@ -14,11 +14,11 @@ from .. import vmat from ..extern import molfile, py3dmol_, rdkit_ -from ..geom import _molsym from ..graph import base as graph_base from ..inchi import base as inchi_base from ..util import ZmatConv, dict_, heuristic, vector, zmat_conv from ..zmat import base as zmat_base +from . import _0molsym from .base import ( central_angle, coordinates, @@ -848,9 +848,9 @@ def external_symmetry_factor(geo, chiral_center=True): if is_atom(geo): ext_sym_fac = 1.0 else: - pg_obj = _molsym.point_group_from_geometry(geo) - ext_sym_fac = _molsym.point_group_symmetry_number(pg_obj) - if _molsym.point_group_is_chiral(pg_obj) and chiral_center: + pg_obj = _0molsym.point_group_from_geometry(geo) + ext_sym_fac = _0molsym.point_group_symmetry_number(pg_obj) + if _0molsym.point_group_is_chiral(pg_obj) and chiral_center: ext_sym_fac *= 0.5 return ext_sym_fac diff --git a/automol/geom/__init__.py b/automol/geom/__init__.py index 6c844a8b..e546a373 100644 --- a/automol/geom/__init__.py +++ b/automol/geom/__init__.py @@ -105,41 +105,41 @@ from .base._vib import vibrational_analysis # L4 # MolSym interface -from ._molsym import point_group_from_geometry +from ._0molsym import point_group_from_geometry # conversion functions: # # conversions -from ._conv import graph -from ._conv import graph_without_stereo -from ._conv import connectivity_graph_deprecated -from ._conv import zmatrix -from ._conv import zmatrix_with_conversion_info -from ._conv import update_zmatrix -from ._conv import amchi -from ._conv import amchi_with_numbers -from ._conv import inchi -from ._conv import inchi_with_numbers -from ._conv import chi -from ._conv import chi_with_sort -from ._conv import smiles -from ._conv import rdkit_molecule -from ._conv import py3dmol_view -from ._conv import display +from ._1conv import graph +from ._1conv import graph_without_stereo +from ._1conv import connectivity_graph_deprecated +from ._1conv import zmatrix +from ._1conv import zmatrix_with_conversion_info +from ._1conv import update_zmatrix +from ._1conv import amchi +from ._1conv import amchi_with_numbers +from ._1conv import inchi +from ._1conv import inchi_with_numbers +from ._1conv import chi +from ._1conv import chi_with_sort +from ._1conv import smiles +from ._1conv import rdkit_molecule +from ._1conv import py3dmol_view +from ._1conv import display # # derived properties -from ._conv import is_connected -from ._conv import linear_atoms -from ._conv import closest_unbonded_atom_distances -from ._conv import could_be_forming_bond -from ._conv import ts_reacting_electron_direction -from ._conv import external_symmetry_factor +from ._1conv import is_connected +from ._1conv import linear_atoms +from ._1conv import closest_unbonded_atom_distances +from ._1conv import could_be_forming_bond +from ._1conv import ts_reacting_electron_direction +from ._1conv import external_symmetry_factor # # derived operations -from ._conv import apply_zmatrix_conversion -from ._conv import undo_zmatrix_conversion -from ._conv import set_distance -from ._conv import set_central_angle -from ._conv import set_dihedral_angle +from ._1conv import apply_zmatrix_conversion +from ._1conv import undo_zmatrix_conversion +from ._1conv import set_distance +from ._1conv import set_central_angle +from ._1conv import set_dihedral_angle # # interfaces -from ._conv import ase_atoms -from ._conv import from_ase_atoms +from ._1conv import ase_atoms +from ._1conv import from_ase_atoms # extra functions: from ._extra import are_torsions_same from ._extra import is_unique diff --git a/automol/geom/_align.py b/automol/geom/_align.py index ca9f7257..be9898a9 100644 --- a/automol/geom/_align.py +++ b/automol/geom/_align.py @@ -6,7 +6,7 @@ import numpy from ..graph.base import _01networkx -from ._conv import graph as graph_conv +from ._1conv import graph as graph_conv from .base import distance_matrix, is_atom, reorder diff --git a/automol/geom/_extra.py b/automol/geom/_extra.py index cf42b8c9..cba9ce7e 100644 --- a/automol/geom/_extra.py +++ b/automol/geom/_extra.py @@ -5,7 +5,7 @@ import numpy from ..graph import base as graph_base -from ._conv import graph, inchi +from ._1conv import graph, inchi from .base import ( almost_equal_coulomb_spectrum, almost_equal_dist_matrix, diff --git a/automol/geom/_ring.py b/automol/geom/_ring.py index eaeaca27..de010b3b 100644 --- a/automol/geom/_ring.py +++ b/automol/geom/_ring.py @@ -8,7 +8,7 @@ from phydat import phycon from ..graph import base as graph_base -from ._conv import ( +from ._1conv import ( graph, ) from .base import ( diff --git a/automol/geom/base/_0core.py b/automol/geom/base/_0core.py index d247f594..153be5f1 100644 --- a/automol/geom/base/_0core.py +++ b/automol/geom/base/_0core.py @@ -131,7 +131,6 @@ def set_coordinates(geo, xyz_dct, angstrom=False): :type geo: automol geometry data structure :param xyz_dct: new xyz values for a set of indices :type xyz_dct: dict[int: tuple(float)] - :rtype: automol molecular graph data structure """ symbs = symbols(geo) From 9ceaa25dd8de928e6567e7b8c4d105a616fabaab Mon Sep 17 00:00:00 2001 From: avcopan Date: Fri, 30 Aug 2024 13:38:40 -0500 Subject: [PATCH 3/5] Refactors vibrational analysis test --- automol/tests/conftest.py | 10 ------- automol/tests/data/oh_freqs.txt | 1 + automol/tests/data/oh_geom.xyz | 4 +++ automol/tests/data/oh_hess.txt | 6 ++++ automol/tests/test_geom.py | 53 +++++++++++---------------------- 5 files changed, 28 insertions(+), 46 deletions(-) delete mode 100644 automol/tests/conftest.py create mode 100644 automol/tests/data/oh_freqs.txt create mode 100644 automol/tests/data/oh_geom.xyz create mode 100644 automol/tests/data/oh_hess.txt diff --git a/automol/tests/conftest.py b/automol/tests/conftest.py deleted file mode 100644 index ab9d8081..00000000 --- a/automol/tests/conftest.py +++ /dev/null @@ -1,10 +0,0 @@ -"""Pytest fixtures.""" - -from pathlib import Path - -import pytest - - -@pytest.fixture -def data_directory_path(): - return Path(__file__).parent / "data" diff --git a/automol/tests/data/oh_freqs.txt b/automol/tests/data/oh_freqs.txt new file mode 100644 index 00000000..0425a66f --- /dev/null +++ b/automol/tests/data/oh_freqs.txt @@ -0,0 +1 @@ + 3776.5 diff --git a/automol/tests/data/oh_geom.xyz b/automol/tests/data/oh_geom.xyz new file mode 100644 index 00000000..a1fdef65 --- /dev/null +++ b/automol/tests/data/oh_geom.xyz @@ -0,0 +1,4 @@ + 2 + +O 0.000000 0.000000 0.108359 +H 0.000000 0.000000 -0.866873 \ No newline at end of file diff --git a/automol/tests/data/oh_hess.txt b/automol/tests/data/oh_hess.txt new file mode 100644 index 00000000..c1c71473 --- /dev/null +++ b/automol/tests/data/oh_hess.txt @@ -0,0 +1,6 @@ +-7.978209999999999631e-05 0.000000000000000000e+00 0.000000000000000000e+00 7.978209999999999631e-05 0.000000000000000000e+00 0.000000000000000000e+00 +0.000000000000000000e+00 2.062429999999999835e-05 0.000000000000000000e+00 0.000000000000000000e+00 -2.062429999999999835e-05 0.000000000000000000e+00 +0.000000000000000000e+00 0.000000000000000000e+00 5.117009999999999614e-01 0.000000000000000000e+00 0.000000000000000000e+00 -5.117009999999999614e-01 +7.978209999999999631e-05 0.000000000000000000e+00 0.000000000000000000e+00 -7.978209999999999631e-05 0.000000000000000000e+00 0.000000000000000000e+00 +0.000000000000000000e+00 -2.062429999999999835e-05 0.000000000000000000e+00 0.000000000000000000e+00 2.062429999999999835e-05 0.000000000000000000e+00 +0.000000000000000000e+00 0.000000000000000000e+00 -5.117009999999999614e-01 0.000000000000000000e+00 0.000000000000000000e+00 5.117009999999999614e-01 diff --git a/automol/tests/test_geom.py b/automol/tests/test_geom.py index ca7bc80b..42432c7c 100644 --- a/automol/tests/test_geom.py +++ b/automol/tests/test_geom.py @@ -1,15 +1,16 @@ """ test automol.geom """ -# import pytest -import io from pathlib import Path import numpy +import pytest import automol from automol import geom +DATA_PATH = Path(__file__).parent / "data" + C2H2CLF_GEO = ( ("F", (2.994881276150, -1.414434615111, -0.807144415388)), ("C", (1.170155936996, 0.359360756989, -0.513323178859)), @@ -784,46 +785,26 @@ def test__repulsion_energy(): assert not automol.geom.has_low_relative_repulsion_energy(geo2, ref_geo, "lj_12_6") -def test__vibrational_analysis(data_directory_path): +@pytest.mark.parametrize( + "smi,geo_file,hess_file,freqs_file", + [ + ("[OH]", "oh_geom.xyz", "oh_hess.txt", "oh_freqs.txt"), + ("OO", "h2o2_geom.xyz", "h2o2_hess.txt", "h2o2_freqs.txt"), + ("C1=CCCC1.[H]", "c4h7_h2_geom.xyz", "c4h7_h2_hess.txt", "c4h7_h2_freqs.txt"), + ], +) +def test__vibrational_analysis(smi, geo_file, hess_file, freqs_file): """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 (H2O2) - geo = geom.from_xyz_string((data_directory_path / "h2o2_geom.xyz").read_text()) - hess = numpy.loadtxt(data_directory_path / "h2o2_hess.txt") - ref_freqs = numpy.loadtxt(data_directory_path / "h2o2_freqs.txt") - freqs, _ = geom.vibrational_analysis(geo, hess) - print(freqs) - assert numpy.allclose(freqs, ref_freqs, atol=1e-1), f"{freqs} !=\n{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") + print(f"Testing frequency projection for {smi}") + geo = geom.from_xyz_string((DATA_PATH / geo_file).read_text()) + hess = numpy.loadtxt(DATA_PATH / hess_file) + ref_freqs = numpy.loadtxt(DATA_PATH / freqs_file) 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__": - data_directory_path = Path(__file__).parent / "data" # __align() # test__change_zmatrix_row_values() # test__inchi_with_sort() @@ -840,4 +821,4 @@ def test__vibrational_analysis(data_directory_path): # test__repulsion_energy() # test__dist_analysis() # test__argunique_coulomb_spectrum() - test__vibrational_analysis(data_directory_path) + test__vibrational_analysis("OO", "h2o2_geom.xyz", "h2o2_hess.txt", "h2o2_freqs.txt") From 1c7a2de98f3046efd7916e11353935db2e6fc457 Mon Sep 17 00:00:00 2001 From: avcopan Date: Fri, 30 Aug 2024 13:58:58 -0500 Subject: [PATCH 4/5] Refactor projection to allow removal of other modes --- automol/geom/_2vib.py | 80 +++++++++++++++++++++++++++++++++++ automol/geom/__init__.py | 4 +- automol/geom/base/_0core.py | 31 -------------- automol/geom/base/__init__.py | 3 -- automol/geom/base/_vib.py | 51 ---------------------- automol/tests/test_geom.py | 20 ++++----- 6 files changed, 92 insertions(+), 97 deletions(-) create mode 100644 automol/geom/_2vib.py delete mode 100644 automol/geom/base/_vib.py diff --git a/automol/geom/_2vib.py b/automol/geom/_2vib.py new file mode 100644 index 00000000..24355fc1 --- /dev/null +++ b/automol/geom/_2vib.py @@ -0,0 +1,80 @@ +"""Vibrational analysis.""" + +from collections.abc import Sequence + +import numpy +from qcelemental import constants as qcc + +from .base import 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 +) -> tuple[tuple[float, ...], numpy.ndarray]: + """Get vibrational modes and frequencies from the Hessian matrix. + + :param geo: The molecular geometry + :param hess: The Hessian matrix + :param trans: Keep translations, instead of removing them? + :param rot: Keep rotations, instead of removing them? + :param wavenum: Return frequencies (cm^-1), instead of force constants (a.u.)? + :return: The vibrational frequencies (or force constants) and normal modes + """ + # 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 + proj = normal_mode_projection(geo, trans=trans, rot=rot) + hess_proj = proj.T @ hess_mw @ proj + + # 2. Compute eigenvalues and eigenvectors of the mass-weighted Hessian matrix + eig_vals, eig_vecs = numpy.linalg.eigh(hess_proj) + + # 3. Un-mass-weight the normal coordinates + # Qmw = PI (see above) Q = Qmw / mw_vec + norm_coos_mw = numpy.dot(proj, eig_vecs) / mw_vec[:, numpy.newaxis] + norm_coos = norm_coos_mw / mw_vec[:, numpy.newaxis] + if not wavenum: + return eig_vals, norm_coos + + # 4. Get wavenumbers from a.u. force constants + har2J = qcc.conversion_factor("hartree", "J") + amu2kg = qcc.conversion_factor("atomic_mass_unit", "kg") + bohr2m = qcc.conversion_factor("bohr", "meter") + sol = qcc.get("speed of light in vacuum") * 100 # in cm / s + to_inv_cm = numpy.sqrt(har2J / (amu2kg * bohr2m * bohr2m)) / (sol * 2 * numpy.pi) + freqs = numpy.sqrt(numpy.complex_(eig_vals)) * to_inv_cm + freqs = tuple(map(float, numpy.real(freqs) - numpy.imag(freqs))) + + return freqs, norm_coos diff --git a/automol/geom/__init__.py b/automol/geom/__init__.py index e546a373..0b875d14 100644 --- a/automol/geom/__init__.py +++ b/automol/geom/__init__.py @@ -101,8 +101,6 @@ from .base._2intmol import has_low_relative_repulsion_energy from .base._2intmol import total_repulsion_energy from .base._2intmol import repulsion_energy -# vibrational analysis -from .base._vib import vibrational_analysis # L4 # MolSym interface from ._0molsym import point_group_from_geometry @@ -140,6 +138,8 @@ # # interfaces from ._1conv import ase_atoms from ._1conv import from_ase_atoms +# vibrational analysis +from ._2vib import vibrational_analysis # extra functions: from ._extra import are_torsions_same from ._extra import is_unique diff --git a/automol/geom/base/_0core.py b/automol/geom/base/_0core.py index 153be5f1..bbdb76fc 100644 --- a/automol/geom/base/_0core.py +++ b/automol/geom/base/_0core.py @@ -697,37 +697,6 @@ def rotational_normal_modes(geo, mass_weight: bool = True) -> numpy.ndarray: 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, drop_null: bool = False): """Calculate the rotational constants. (these not sorted in to A,B,C). diff --git a/automol/geom/base/__init__.py b/automol/geom/base/__init__.py index bcbe68db..38fccdae 100644 --- a/automol/geom/base/__init__.py +++ b/automol/geom/base/__init__.py @@ -93,7 +93,6 @@ from ._2intmol import has_low_relative_repulsion_energy from ._2intmol import total_repulsion_energy from ._2intmol import repulsion_energy -from ._vib import vibrational_analysis __all__ = [ @@ -185,6 +184,4 @@ 'has_low_relative_repulsion_energy', 'total_repulsion_energy', 'repulsion_energy', - # vibrational analysis - 'vibrational_analysis', ] diff --git a/automol/geom/base/_vib.py b/automol/geom/base/_vib.py deleted file mode 100644 index 5252cfcb..00000000 --- a/automol/geom/base/_vib.py +++ /dev/null @@ -1,51 +0,0 @@ -"""Vibrational analysis.""" - -from collections.abc import Sequence - -import numpy -from qcelemental import constants as qcc - -from ._0core import masses, projection_onto_internal_modes - -MatrixLike = Sequence[Sequence[float]] | numpy.ndarray - - -def vibrational_analysis( - geo, hess: MatrixLike, freq: bool = True -) -> tuple[tuple[float, ...], numpy.ndarray]: - """Get vibrational modes and frequencies from the Hessian matrix. - - :param geo: The molecular geometry - :param hess: The Hessian matrix - :param freq: Return frequencies (cm^-1), instead of force constants (a.u.)? - :return: The vibrational frequencies (or force constants) and normal modes - """ - # 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_int) - - # 3. un-mass-weight the normal coordinates - # 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 - - # 4. get wavenumbers from a.u. force constants - har2J = qcc.conversion_factor("hartree", "J") - amu2kg = qcc.conversion_factor("atomic_mass_unit", "kg") - bohr2m = qcc.conversion_factor("bohr", "meter") - sol = qcc.get("speed of light in vacuum") * 100 # in cm / s - to_inv_cm = numpy.sqrt(har2J / (amu2kg * bohr2m * bohr2m)) / (sol * 2 * numpy.pi) - freqs = numpy.sqrt(numpy.complex_(eig_vals)) * to_inv_cm - freqs = tuple(map(float, numpy.real(freqs) - numpy.imag(freqs))) - - return freqs, norm_coos diff --git a/automol/tests/test_geom.py b/automol/tests/test_geom.py index 42432c7c..044ce8f4 100644 --- a/automol/tests/test_geom.py +++ b/automol/tests/test_geom.py @@ -786,19 +786,19 @@ def test__repulsion_energy(): @pytest.mark.parametrize( - "smi,geo_file,hess_file,freqs_file", + "smi,fml", [ - ("[OH]", "oh_geom.xyz", "oh_hess.txt", "oh_freqs.txt"), - ("OO", "h2o2_geom.xyz", "h2o2_hess.txt", "h2o2_freqs.txt"), - ("C1=CCCC1.[H]", "c4h7_h2_geom.xyz", "c4h7_h2_hess.txt", "c4h7_h2_freqs.txt"), + ("[OH]", "oh"), + ("OO", "h2o2"), + ("C1=CCCC1.[H]", "c4h7_h2"), ], ) -def test__vibrational_analysis(smi, geo_file, hess_file, freqs_file): +def test__vibrational_analysis(smi, fml): """Test automol.geom.vibrational_analysis.""" - print(f"Testing frequency projection for {smi}") - geo = geom.from_xyz_string((DATA_PATH / geo_file).read_text()) - hess = numpy.loadtxt(DATA_PATH / hess_file) - ref_freqs = numpy.loadtxt(DATA_PATH / freqs_file) + print(f"Testing frequency projection for {fml} (SMILES {smi})") + geo = geom.from_xyz_string((DATA_PATH / f"{fml}_geom.xyz").read_text()) + hess = numpy.loadtxt(DATA_PATH / f"{fml}_hess.txt") + ref_freqs = numpy.loadtxt(DATA_PATH / f"{fml}_freqs.txt") freqs, _ = geom.vibrational_analysis(geo, hess) print(freqs) assert numpy.allclose(freqs, ref_freqs, atol=1e-1), f"{freqs} !=\n{ref_freqs}" @@ -821,4 +821,4 @@ def test__vibrational_analysis(smi, geo_file, hess_file, freqs_file): # test__repulsion_energy() # test__dist_analysis() # test__argunique_coulomb_spectrum() - test__vibrational_analysis("OO", "h2o2_geom.xyz", "h2o2_hess.txt", "h2o2_freqs.txt") + test__vibrational_analysis("OO", "h2o2") From 2fbca10e7387682ea96277d21cff72777ce73ad8 Mon Sep 17 00:00:00 2001 From: avcopan Date: Fri, 30 Aug 2024 15:58:34 -0500 Subject: [PATCH 5/5] 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)