Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

fix projected frequencies #368

Open
wants to merge 2 commits into
base: v1.4.5
Choose a base branch
from
Open

fix projected frequencies #368

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
640e135
fix projection
t-young31 Nov 24, 2024
4e6a7c1
self review
t-young31 Nov 24, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
86 changes: 52 additions & 34 deletions autode/hessians.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,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 = [
Expand Down Expand Up @@ -151,20 +153,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)
Expand All @@ -174,14 +165,10 @@ 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()

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

return t1, t2, t3, np.array(t4), np.array(t5), np.array(t6)

Expand Down Expand Up @@ -213,19 +200,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))
shoubhikraj marked this conversation as resolved.
Show resolved Hide resolved
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

My suggestion editing from this line:

 tr_bas = np.array(tr_vecs).transpose()
U_s, s_v, _ = np.linalg.svd(tr_bas, full_matrices=True)

if self.atoms.are_linear() and s_v[5]/s_v[0] > 1.e-4:
    logger.warning("Molecule detected as linear, but lowest singular "
                             "value from trans. rot. basis is {s_v[5]:.4e}")

return U_s

If the rank of matrix is k, then first k columns of SVD will be the orthonormal basis for the column space. Since we started from only the translation and rotation vectors, the first 5 or 6 columns of U will be translation and rotation basis - and the rest will be orthogonal vectors corresponding to vibration (since U always produces an orthonormal basis, as far as i understand!)


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]
shoubhikraj marked this conversation as resolved.
Show resolved Hide resolved

Expand Down Expand Up @@ -400,13 +389,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
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this whole part is not needed, and the original code was fine(?). Projecting the Hessian should mean that we are safe to assume the first 5 or 6 columns to be zero. If they are not zero - then that is probably the measure of how much numerical error there is in the Hessian. (btw, I have been trying to understand the math - and so far I have been able to mathematically prove that the first 5 / 6 columns of the projected Hessian definitely have to be zero. For the first 5 / 6 rows it also makes intuitive sense - but I haven't been able to find the mathematical proof yet)

norms = np.linalg.norm(H, axis=0)
shoubhikraj marked this conversation as resolved.
Show resolved Hide resolved
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
shoubhikraj marked this conversation as resolved.
Show resolved Hide resolved
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same with this part - I think it is not required.

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":
Expand Down Expand Up @@ -591,7 +609,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()

Expand Down
4 changes: 1 addition & 3 deletions autode/wrappers/G09.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
1 change: 1 addition & 0 deletions doc/changelog.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
45 changes: 42 additions & 3 deletions tests/test_hessian.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
12 changes: 8 additions & 4 deletions tests/test_wrappers/test_gaussian.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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():
Expand Down
Loading