fix projection

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()
 

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

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

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():