Vibrational frequencies (via ASE)

vibe-qc doesn’t ship an analytic Hessian yet, but ASE’s ase.vibrations.Vibrations builds one by finite differences of the forces — and vibe-qc’s analytic forces go through the same VibeQC ASE calculator that drives geometry optimisation. That gives you the three basic products:

  • normal-mode frequencies (in cm⁻¹)

  • zero-point vibrational energy

  • per-mode atomic displacement vectors (for animation / IR prediction)

for any molecule vibe-qc can treat (RHF, UHF, RKS, UKS).

Workflow

Three steps:

  1. Optimise the geometry. Hessian at a non-stationary point is meaningless — you’ll get imaginary translations / rotations mixed into the normal modes.

  2. Run the finite-difference Hessian. ASE displaces each atom in each Cartesian direction, recomputes forces, builds the numerical Hessian, and mass-weights it.

  3. Diagonalise to get frequencies and normal modes.

Example: water

from ase import Atoms
from ase.optimize import BFGSLineSearch
from ase.vibrations import Vibrations
from vibeqc.ase import VibeQC

atoms = Atoms(
    symbols=["O", "H", "H"],
    positions=[
        [ 0.0,  0.00,  0.00],
        [ 0.0,  0.76,  0.59],
        [ 0.0, -0.76,  0.59],
    ],
)
atoms.calc = VibeQC(basis="sto-3g")

# 1. Relax — Hessian is only meaningful at a stationary point
BFGSLineSearch(atoms, logfile=None).run(fmax=0.01)

# 2. Finite-difference Hessian via ASE. ``name=`` is a cache prefix;
#    rerunning with the same name reuses previous displacement results.
vib = Vibrations(atoms, name="water_vib")
vib.run()

# 3. Summary table
vib.summary()

Typical output (HF/STO-3G, after BFGS to fmax = 0.01 eV/Å):

---------------------
  #    meV     cm^-1
---------------------
  0    0.0      0.0      <- translation
  1    0.0      0.2      <- translation
  2    4.3     34.3      <- rotation
  3    4.6     37.4      <- rotation
  4    0.1i     0.8i     <- near-zero modes (numerical noise)
  5    6.0i    48.8i
  6  268.9  2169.2       HOH bend
  7  513.2  4139.6       symmetric O-H stretch
  8  544.4  4390.6       antisymmetric O-H stretch
---------------------
Zero-point energy: 0.668 eV

The three real internal modes match the chemistry: a bend near 2000 cm⁻¹ and two stretches near 4000 cm⁻¹. HF/STO-3G over-estimates by 15–35% — the usual scaling factor for this level of theory is ~0.89. Multiplying the printed frequencies by 0.89 gives 1930, 3684, 3907 cm⁻¹ — within 2% of the experimental water frequencies (1595, 3657, 3756 cm⁻¹). Bigger basis + DFT gets you closer without the scaling factor.

Pulling the numpy arrays

vib.get_frequencies() returns complex frequencies in cm⁻¹ (imaginary ones signal saddle points or numerical noise below the threshold). vib.get_energies() gives meV.

import numpy as np

freqs = vib.get_frequencies()            # complex-valued ndarray
zpe = vib.get_zero_point_energy()        # eV

print("Real frequencies (cm^-1):",
      np.real(freqs[np.abs(np.imag(freqs)) < 1e-6]))
print(f"ZPE = {zpe * 96.485:.3f} kJ/mol")

Animating a mode

To write an XYZ-frame animation of, say, mode 8 (antisymmetric stretch):

vib.write_mode(n=8, kT=None, nimages=30)
# -> water_vib.8.traj — open with `ase gui water_vib.8.traj`

You can also render each displaced geometry to XYZ and feed it into your favourite orbital viewer (Jmol, Avogadro, VMD) for a publication-quality animation.

Caveats

  • Finite differences are slow. A molecule with N atoms needs 6N single-point force evaluations (±δ along 3N directions) on top of the optimisation. Water is 18 force calls; a 20-atom molecule is 120. That’s usable for small organics, painful for anything bigger.

  • Step size matters. ASE’s default displacement is 0.01 Å, which is usually fine. If a mode looks noisy, reduce it (Vibrations(..., delta=0.005)) and re-run.

  • Results depend on SCF tolerance. Tighten conv_tol_grad on the RHF/RKS options for cleaner low-frequency modes; the default 1e-6 starts to hit the floor around 50 cm⁻¹.

  • No IR intensities yet. Those need dipole-moment derivatives; the next wishlist item after an analytic Hessian ships.

Next

Once you have vibrational frequencies, you can compute thermodynamic functions (ZPE, internal energy, entropy, enthalpy, Gibbs free energy at temperature) via ASE’s IdealGasThermo — see thermodynamics.