Vibrational analysis

Harmonic vibrational analysis diagonalizes the mass-weighted Hessian to calculate frequencies and normal modes.

Water frequencies

Three atoms, three vibrations. Relaxed and analysed with GFN2-xTB and again with PBE/def2-SVP:

Mode	GFN2-xTB	PBE/def2-SVP	Experiment
bend	1538.6 (−56.2)	1608.4 (+13.6)	1594.75
symmetric stretch	3643.4 (−13.6)	3689.7 (+32.7)	3657.05
antisymmetric stretch	3651.7 (−104.3)	3788.8 (+32.9)	3755.93

All values in cm⁻¹, deviation from experiment in brackets. The experimental numbers are the gas-phase fundamentals of H₂¹⁶O from Benedict, Gailar and Plyler, J. Chem. Phys. 24, 1139 (1956).

Read that comparison carefully, because it is not a clean like-for-like. A fundamental is the observed 0 → 1 transition and already contains anharmonicity, which pulls an O-H stretch down by roughly 150 cm⁻¹ from its harmonic value. Both calculations here are harmonic. So a harmonic number that lands on the fundamental is partly cancelling error, and the honest reading is the pattern rather than any single line.

The pattern is this. GFN2-xTB is low across the board, by 56 cm⁻¹ on the bend and 104 cm⁻¹ on the antisymmetric stretch. PBE is high by +14 to +33 cm⁻¹, which is what a harmonic treatment of a functional that slightly overbinds looks like once anharmonicity is left out.

The more useful failure is in the splitting. Experiment puts the two O-H stretches 99 cm⁻¹ apart. PBE gives 99 cm⁻¹. GFN2-xTB gives 8.2 cm⁻¹, so its two stretches are effectively degenerate. If you are using the tight-binding method to pre-screen conformers or to get a zero-point energy, that hardly matters. If you are trying to assign a spectrum, it makes the method useless for the O-H region, and no amount of tighter convergence will fix it, because it is the parametrization and not the numerics.

Zero-point energy from the same runs: 0.5529 eV with GFN2-xTB, 0.5642 eV with PBE. That is half the sum of the mode energies, so it inherits the same errors.

Normal modes

Ordering by frequency is not the same as knowing which mode is which. The assignment below is read off the eigenvectors: whether the two O-H bonds lengthen together or in opposition. Displacements are exaggerated so they are visible; the real zero-point amplitude of an O-H stretch is around 0.07 Å.

Bend, 1539 cm⁻¹. The bonds barely change length; the angle opens and closes.

Symmetric stretch, 3643 cm⁻¹. Both bonds lengthen and shorten together.

Antisymmetric stretch, 3652 cm⁻¹. One bond lengthens while the other shortens.

Frequency calculation

from atomli.build import molecule
from atomli.calculators.xtb import XTB
from atomli.optimize import BFGS
from atomli.vibrations import Vibrations

atoms = molecule("H2O")
atoms.calc = XTB(method="gfn2")
BFGS(atoms).run(fmax=1e-4)          # relax first, always

vib = Vibrations(atoms)
vib.run()
vib.summary()

Energies come back in electronvolt and frequencies in reciprocal centimetre. Imaginary modes carry their magnitude in the imaginary part of the returned complex array, which is what summary marks with a trailing i.

On this machine that whole sequence costs 0.67 ms for the single point, 0.89 ms for the relaxation and 0.83 ms for the Hessian, each measured from a freshly built molecule and calculator. At three atoms that construction is most of the number and the SCF itself is a rounding error, which is the point: a tight-binding frequency job on a small molecule is not something you need to plan around.

Geometry relaxation

The relaxation in that snippet is not a formality. A Hessian is well defined at any geometry, but only at a stationary point does its spectrum mean frequencies. Everywhere else the six rigid-body modes, which should sit at zero because translating or rotating a molecule costs nothing, pick up the residual gradient and come back as large numbers that look exactly like real vibrations.

Here are those six modes for the same molecule at three geometries, largest magnitude shown:

Geometry	Residual force	Largest rigid-body mode
relaxed to fmax 1e-4	9.0e-5 eV/Å	56.3 cm⁻¹
ASE database geometry, unrelaxed	0.75 eV/Å	440.1 cm⁻¹
one O-H pulled out to 1.20 Å	3.89 eV/Å	1078.9 cm⁻¹

The database geometry is a perfectly reasonable guess, and it is still 0.75 eV/Å from the GFN2 minimum. That is enough to push a rigid-body mode to 440 cm⁻¹, which is in the range where a real low-frequency torsion would live. Nothing in the output says it is spurious.

The stretched case is worse in a way that matters more. Its three highest modes come out at 1307, 1602, 3502 cm⁻¹, and the rigid-body contamination now reaches 1079 cm⁻¹, so the two sets overlap and there is no longer any way to tell from the numbers alone which is which. A frequency list from an unrelaxed structure is not approximately right. It is a different quantity.

Check it rather than assume it:

import numpy as np
assert np.abs(atoms.get_forces()).max() < 1e-4

Imaginary frequencies

A negative Hessian eigenvalue gives an imaginary frequency, and the mode it belongs to is a direction in which the energy goes down. At a converged minimum there should be none. When there is one, either the geometry is not stationary, as above, or the structure is a genuine saddle point.

Linear water is the clean example. At exactly 180° with equal bond lengths the bending gradient vanishes by symmetry, so relaxing the bond length alone lands on a true stationary point:

Quantity	Value
O-H length at the stationary point	0.9247 Å
residual force	1.1e-14 eV/Å
energy above the bent minimum	1.358 eV
imaginary modes	1982.2i, 1982.2i cm⁻¹

Two imaginary modes, not one, and they are degenerate to 3e-11 cm⁻¹. That is correct and worth understanding: a linear molecule bends in two independent planes, and both directions run downhill from here. So this is a second-order saddle, not a transition state. A transition state has exactly one imaginary mode. If you are hunting for one and find two, you are on the wrong stationary point.

The residual force of 1e-14 eV/Å is the reason the imaginary frequencies here are trustworthy. Compare it with the unrelaxed cases above: there, imaginary or inflated modes were an artefact; here they are the physics. The only thing separating the two situations is whether the gradient vanished.

Mode count

Vibrations returns 3N numbers for N atoms, one per Cartesian degree of freedom. Six of them are not vibrations: three translations of the whole molecule and three rotations of it. That leaves 3N − 6 real vibrations. A linear molecule has only two distinguishable rotations, because spinning it about its own axis moves nothing, so it keeps one more vibration and the count is 3N − 5.

Molecule	Shape	3N results	Translations + rotations	Vibrations
H2O	bent	9	6	3
CO2	linear	9	5	4

Both have three atoms and both return nine numbers; they differ in how many of those nine are real vibrations. Relaxed CO₂ with GFN2-xTB gives 601, 601, 1425, 2594 cm⁻¹ for its four, and the first two are the degenerate bending pair, the same two-planes argument as linear water except that here the molecule is a minimum and they come out real.

Atomli does not project the rigid-body modes out, so all 3N appear in the returned array and you identify them by magnitude. At the water minimum the six rigid-body modes are all below 56 cm⁻¹ while the smallest true vibration is 1539 cm⁻¹, a gap of more than a factor of twenty, so the split is unambiguous. Some of those six come back marginally imaginary. That is finite-difference noise on a mode whose true frequency is zero, and it is expected.

Hessian methods

Molecular quantum-chemistry calculations carry an analytical second derivative. Everything else has to displace atoms and difference the forces. Vibrations picks the route from what the attached calculator says it can produce, so every calculator works:

method	What runs
"auto" (default)	Analytical when the calculator has a Hessian for this system, finite differences otherwise
"analytical"	Analytical, or an error naming why it is unavailable
"finite_difference"	Displacements, even where an analytical Hessian exists

After run(), vib.method_used reports which one ran:

vib = Vibrations(atoms).run()
vib.method_used            # "analytical"

An analytical Hessian is available for molecular QC calculations other than SKALA. Force fields, machine-learning potentials, xTB, MLIP, SKALA, and any periodic cell go through finite differences, which needs no special support from the calculator at all.

Both routes are the same physics, and on relaxed water they agree to 1.06 cm⁻¹ across all three vibrations:

analytical = Vibrations(atoms).run()
numeric = Vibrations(atoms, method="finite_difference").run()

They do not cost the same, and not in the direction you would guess. Measured on the same geometries, same machine:

Molecule	Atoms	Analytical	Finite difference
H2O	3	4.0 s	0.6 s
CH4	5	16.4 s	1.6 s

The displacement route wins by roughly 7× on water and 10× on methane, and the gap widens with size rather than closing. Finite differences on N atoms is 6N single points that each reuse the converged wavefunction of the reference geometry, which is cheap; the coupled-perturbed solve behind the analytical Hessian is not yet competitive with that in this engine. Take the default unless you specifically need the analytical route, and if a Hessian is the bottleneck in your workflow, measure both before assuming.

Finite-difference controls

delta is the displacement in angstrom and nfree is the number of displacements per degree of freedom, 2 or 4. They only affect the finite-difference route.

vib = Vibrations(atoms, delta=0.005, nfree=4, method="finite_difference")

indices restricts the analysis to a subset of atoms, which is how a partial Hessian for an adsorbate on a frozen slab is built. Restricting the set selects finite differences, because a partial Hessian is a displacement quantity.

vib = Vibrations(atoms, indices=[12, 13, 14])
vib.run()
vib.get_frequencies().shape    # (9,) -- three atoms x three directions

A partial Hessian has no rigid-body modes to discard, because the frozen atoms hold the fragment in place. Do not subtract six from that count.

Modes and results

vib.get_energies()       # complex, eV, ascending
vib.get_frequencies()    # complex, cm^-1, ascending
vib.get_zero_point_energy()
vib.get_mode(8)          # (n_atoms, 3) Cartesian displacement
vib.get_hessian_2d()     # (3 n_active, 3 n_active), eV/angstrom^2

get_mode is what the figures above are built from: take the displacement array, scale it, and add it to the equilibrium positions at a few phases.

import numpy as np

mode = vib.get_mode(8)
scale = 0.15 / np.linalg.norm(mode, axis=1).max()   # peak displacement, in angstrom

frames = [atoms.copy() for _ in range(24)]
for index, frame in enumerate(frames):
    phase = np.sin(2 * np.pi * index / len(frames))
    frame.set_positions(atoms.get_positions() + scale * phase * mode)

Normalize by the row norms, not by np.abs(mode).max(). The largest single Cartesian component understates the atom’s actual displacement by up to a factor of √3, which is enough to turn a stretch into a structure that visibly falls apart.

write_mode does the same thing straight to disk, one file per mode as vib.<n>.xyz, or every non-zero mode when no index is given:

vib.write_mode(8)

Hessian API

Calculators that have an analytical Hessian expose it directly, shaped (n_atoms, 3, n_atoms, 3) in electronvolt per angstrom squared:

hessian = atoms.calc.get_hessian(atoms)

Calculators without one raise ASE’s PropertyNotImplementedError. Reach for Vibrations instead when you want frequencies regardless of the backend.