From 640e13582bc7417b0fca527e2eaf7e93b36accca Mon Sep 17 00:00:00 2001 From: t-young31 Date: Sun, 24 Nov 2024 15:57:42 +0000 Subject: [PATCH 1/2] fix projection --- autode/hessians.py | 91 +++++++++++++++++----------- autode/wrappers/G09.py | 4 +- tests/test_hessian.py | 45 +++++++++++++- tests/test_wrappers/test_gaussian.py | 12 ++-- 4 files changed, 108 insertions(+), 44 deletions(-) diff --git a/autode/hessians.py b/autode/hessians.py index c8ef05be6..405c0c94d 100644 --- a/autode/hessians.py +++ b/autode/hessians.py @@ -15,6 +15,7 @@ Union, TYPE_CHECKING, ) + from autode.wrappers.keywords import Functional, GradientKeywords from autode.log import logger from autode.config import Config @@ -38,6 +39,8 @@ from autode.wrappers.keywords import GradientKeywords, Keywords from autode.values import Distance, Gradient +THRESHOLD_ROTATION_FREQ = Frequency(50, "cm-1") + class Hessian(ValueArray): implemented_units = [ @@ -151,20 +154,9 @@ def _tr_vecs( """ n_atoms = len(self.atoms) - if n_atoms > 2: - # Get an orthonormal basis shifted from the principal rotation axis - _rot_M = np.array( - [ - [1.0, 0.0, 0.0], - [0.0, 0.09983341664682815, -0.9950041652780258], - [0.0, 0.9950041652780258, 0.09983341664682815], - ] - ) - - _, (e_x, e_y, e_z) = np.linalg.eigh(_rot_M.dot(self.atoms.moi)) - else: - # Get a random orthonormal basis in 3D - (e_x, e_y, e_z), _ = np.linalg.qr(np.random.rand(3, 3)) + e_x = np.array([1.0, 0.0, 0.0]) + e_y = np.array([0.0, 1.0, 0.0]) + e_z = np.array([0.0, 0.0, 1.0]) t1 = np.tile(e_x, reps=n_atoms) t2 = np.tile(e_y, reps=n_atoms) @@ -174,16 +166,16 @@ def _tr_vecs( t4, t5, t6 = [], [], [] for atom in self.atoms: - t4 += np.cross(e_x, atom.coord - com).tolist() - t5 += np.cross(e_y, atom.coord - com).tolist() - t6 += np.cross(e_z, atom.coord - com).tolist() + r = atom.coord - com + t4 += np.cross(e_x, r).tolist() + t5 += np.cross(e_y, r).tolist() + t6 += np.cross(e_z, r).tolist() - if any(np.isclose(np.linalg.norm(t_i), 0.0) for t_i in (t4, t5, t6)): - # Found linear dependency in rotation vectors, attempt to remove - # by initialising different random orthogonal vectors - return self._tr_vecs() + t4 = np.array(t4) + t5 = np.array(t5) + t6 = np.array(t6) - return t1, t2, t3, np.array(t4), np.array(t5), np.array(t6) + return t1, t2, t3, t4, t5, t6 @cached_property def _proj_matrix(self) -> np.ndarray: @@ -213,19 +205,21 @@ def _proj_matrix(self) -> np.ndarray: # Construct M^1/2, which as it's diagonal, is just the roots of the # diagonal elements - masses = np.repeat( - [atom.mass for atom in self.atoms], repeats=3, axis=np.newaxis - ) - m_half = np.diag(np.sqrt(masses)) + masses = np.repeat([atom.mass for atom in self.atoms], repeats=3) + m_half = np.sqrt(masses) for t_i in (t1, t2, t3, t4, t5, t6): - t_i[:] = np.dot(m_half, np.array(t_i)) - t_i /= np.linalg.norm(t_i) + t_i *= m_half # Generate a transform matrix D with the first columns as translation/ # rotation vectors with the remainder as random orthogonal columns - M = np.random.rand(3 * len(self.atoms), 3 * len(self.atoms)) - 0.5 - M[:, :6] = np.column_stack((t1, t2, t3, t4, t5, t6)) + M = np.eye(3 * len(self.atoms)) + + i = 0 + for t_i in (t1, t2, t3, t4, t5, t6): + if not np.isclose(np.linalg.norm(t_i), 0): + M[:, i] = t_i + i += 1 return np.linalg.qr(M)[0] @@ -400,13 +394,42 @@ def frequencies_proj(self) -> List[Frequency]: "Could not calculate projected frequencies, must " "have atoms set" ) - n_tr = self.n_tr # Number of translational+rotational modes - lambdas = np.linalg.eigvalsh(self._proj_mass_weighted[n_tr:, n_tr:]) - trans_rot_freqs = [Frequency(0.0) for _ in range(n_tr)] + H = self._proj_mass_weighted + norms = np.linalg.norm(H, axis=0) + max_norm = np.max(norms) + n_zeroed_modes = 0 # Number of modes that have been well projected out of the hessian + for norm in norms: + if norm / max_norm > 0.1 or n_zeroed_modes == n_tr: + break + else: + n_zeroed_modes += 1 + + if n_zeroed_modes != n_tr: + logger.warn( + f"Number of well zeroed eigenvectors of the hessian " + f"was [{n_zeroed_modes}] should be [{n_tr}]" + ) + + lambdas = np.linalg.eigvalsh(H[n_zeroed_modes:, n_zeroed_modes:]) + trans_rot_freqs = [Frequency(0.0) for _ in range(n_zeroed_modes)] vib_freqs = self._eigenvalues_to_freqs(lambdas) + n_rotational_modes_in_vib_freqs = n_tr - n_zeroed_modes + for i in np.argsort(np.abs(np.array(vib_freqs))): + freq = vib_freqs[i] + if ( + n_rotational_modes_in_vib_freqs > 0 + and freq.real < THRESHOLD_ROTATION_FREQ + ): + logger.warning( + "Found a vibrational mode that should be a rotation with " + f"frequency [{freq}] cm-1. Forcing to zero" + ) + vib_freqs[i] = Frequency(0.0) + n_rotational_modes_in_vib_freqs -= 1 + return trans_rot_freqs + vib_freqs def copy(self, *args, **kwargs) -> "Hessian": @@ -591,7 +614,7 @@ def _gradient(self, species) -> "Gradient": return species.gradient.flatten() @property - def _init_gradient(self) -> "Gradient": + def _init_gradient(self) -> "np.ndarray": """Gradient at the initial geometry of the species""" return np.array(self._species.gradient).flatten() diff --git a/autode/wrappers/G09.py b/autode/wrappers/G09.py index b2cfcbf05..b18e85e9d 100644 --- a/autode/wrappers/G09.py +++ b/autode/wrappers/G09.py @@ -546,9 +546,7 @@ def coordinates_from(self, calc: "CalculationExecutor") -> Coordinates: coords: List[List[float]] = [] for i, line in enumerate(calc.output.file_lines): - if "Input orientation" in line or ( - "Standard orientation" in line and len(coords) == 0 - ): + if "Input orientation" in line or "Standard orientation" in line: coords.clear() xyz_lines = calc.output.file_lines[ i + 5 : i + 5 + calc.molecule.n_atoms diff --git a/tests/test_hessian.py b/tests/test_hessian.py index 2a3d63396..0c54b550d 100644 --- a/tests/test_hessian.py +++ b/tests/test_hessian.py @@ -3,6 +3,7 @@ import pickle import numpy as np import autode as ade +from scipy.stats import special_ortho_group from autode.utils import work_in_tmp_dir, ProcessPool from . import testutils import multiprocessing as mp @@ -588,15 +589,13 @@ def test_nwchem_hessian_co2(): keywords=ade.HessianKeywords(), ) calc.set_output_filename("CO2_hess_nwchem.out") - print(co2.hessian) - print(co2.hessian._mass_weighted) assert_correct_co2_frequencies( hessian=co2.hessian, expected=(659.76, 1406.83, 2495.73) ) @testutils.work_in_zipped_dir(os.path.join(here, "data", "hessians.zip")) -def test_imag_mode(): +def test_sn2_imag_mode(): """ Ensure the imaginary mode for an SN2 reaction is close to that obtained from ORCA by checking forwards (f) and backwards (b) displaced geometries @@ -917,3 +916,43 @@ def test_hessian_pickle_and_unpickle(): assert reloaded_hessian.shape == (3 * mol.n_atoms, 3 * mol.n_atoms) assert reloaded_hessian.atoms == mol.atoms + + +def test_hessian_proj_freqs_acetylene(): + # fmt: off + raw_hessian = np.array([ + [5.5863, -0.0, -0.0, -4.1043, 0.0, 0.0, -1.5072, 0.0, 0.0, 0.0252, 0.0, 0.0], + [-0.0, 0.1913, 0.0, 0.0, -0.115, 0.0, 0.0, -0.0932, 0.0, 0.0, 0.017, 0.0], + [-0.0, 0.0, 0.1913, 0.0, 0.0, -0.115, 0.0, 0.0, -0.0932, 0.0, 0.0, 0.017], + [-4.1043, 0.0, 0.0, 5.5805, -0.0, -0.0, 0.0253, 0.0, 0.0, -1.5014, 0.0, 0.0], + [0.0, -0.115, 0.0, -0.0, 0.192, 0.0, 0.0, 0.017, 0.0, 0.0, -0.094, 0.0], + [0.0, 0.0, -0.115, -0.0, 0.0, 0.192, 0.0, 0.0, 0.017, 0.0, 0.0, -0.094], + [-1.5072, 0.0, 0.0, 0.0253, 0.0, 0.0, 1.4822, -0.0, -0.0, -0.0003, 0.0, 0.0], + [0.0, -0.0932, 0.0, 0.0, 0.017, 0.0, -0.0, 0.0579, 0.0, 0.0, 0.0183, 0.0], + [0.0, 0.0, -0.0932, 0.0, 0.0, 0.017, -0.0, 0.0, 0.0579, 0.0, 0.0, 0.0183], + [0.0252, 0.0, 0.0, -1.5014, 0.0, 0.0, -0.0003, 0.0, 0.0, 1.4765, -0.0, -0.0], + [0.0, 0.017, 0.0, 0.0, -0.094, 0.0, 0.0, 0.0183, 0.0, -0.0, 0.0587, 0.0], + [0.0, 0.0, 0.017, 0.0, 0.0, -0.094, 0.0, 0.0, 0.0183, -0.0, 0.0, 0.0587] + ]) + # fmt: on + atoms = Atoms( + [ + Atom("C", 0.0000, 0.0000, -0.6043), + Atom("C", 0.0000, 0.0000, 0.6042), + Atom("H", 0.0000, 0.0000, -1.6791), + Atom("H", 0.0000, 0.0000, 1.6797), + ] + ) + + for _ in range(100): + h = Hessian(atoms=atoms, input_array=raw_hessian) + freqs = h.frequencies_proj + assert len([v for v in freqs if abs(v.to("cm-1")) < 10.0]) == 5 + assert ( + len([v for v in freqs if np.isclose(v.to("cm-1"), 694.4, atol=1)]) + == 2 + ) + assert len([v for v in freqs if 10 < v.to("cm-1") < 600]) == 0 + + rotation_matrix = special_ortho_group.rvs(3) + atoms.coordinates = np.dot(rotation_matrix, atoms.coordinates.T).T diff --git a/tests/test_wrappers/test_gaussian.py b/tests/test_wrappers/test_gaussian.py index a116ad6d7..963941fa9 100644 --- a/tests/test_wrappers/test_gaussian.py +++ b/tests/test_wrappers/test_gaussian.py @@ -204,13 +204,17 @@ def test_gauss_optts_calc(): assert bond_added test_mol.calc_thermo(calc=calc, ss="1atm", lfm_method="igm") - assert calc.terminated_normally assert calc.optimiser.converged - assert len(test_mol.imaginary_frequencies) == 1 + assert test_mol.imaginary_frequencies is not None + + # NOTE: autodE does not remove eigenmodes that are not + # well zeroed from projected hessian. This is in contrast to + # other electronic structure codes + assert len(test_mol.imaginary_frequencies) == 2 - assert -40.324 < test_mol.free_energy < -40.322 - assert -40.301 < test_mol.enthalpy < -40.298 + assert -40.324 < test_mol.free_energy < -40.320 + assert -40.300 < test_mol.enthalpy < -40.297 def test_bad_gauss_output(): From 4e6a7c16e621f9580284143e9296ac5d9e76bdf7 Mon Sep 17 00:00:00 2001 From: t-young31 Date: Sun, 24 Nov 2024 16:02:21 +0000 Subject: [PATCH 2/2] self review --- autode/hessians.py | 7 +------ doc/changelog.rst | 1 + 2 files changed, 2 insertions(+), 6 deletions(-) diff --git a/autode/hessians.py b/autode/hessians.py index 405c0c94d..3684f5a1c 100644 --- a/autode/hessians.py +++ b/autode/hessians.py @@ -15,7 +15,6 @@ Union, TYPE_CHECKING, ) - from autode.wrappers.keywords import Functional, GradientKeywords from autode.log import logger from autode.config import Config @@ -171,11 +170,7 @@ def _tr_vecs( t5 += np.cross(e_y, r).tolist() t6 += np.cross(e_z, r).tolist() - t4 = np.array(t4) - t5 = np.array(t5) - t6 = np.array(t6) - - return t1, t2, t3, t4, t5, t6 + return t1, t2, t3, np.array(t4), np.array(t5), np.array(t6) @cached_property def _proj_matrix(self) -> np.ndarray: diff --git a/doc/changelog.rst b/doc/changelog.rst index 38d9b3e3c..e6dcb7389 100644 --- a/doc/changelog.rst +++ b/doc/changelog.rst @@ -12,6 +12,7 @@ Functionality improvements Bug Fixes ********* - Fixes coordinate extraction in some G16 output files +- Fixes the projected frequency calculation to make it rotation invariant by not discarding non-zeroed modes from the projected Hessian 1.4.4