Hessians, vibrational frequencies & thermochemistry

Hessians, vibrational frequencies & thermochemistry#

The gradient told you whether you are at equilibrium. The Hessian — the matrix of second derivatives of the energy — tells you what kind of stationary point it is, and unlocks the two things chemists actually measure at a minimum: vibrational frequencies (infrared/Raman spectra) and thermochemistry (the free energy that predicts equilibria and rates). This is a real, tested qc-rs feature (m.scf.hessian, qc.thermo.frequencies) that the rest of the manual has not yet covered.

Theory: why the second derivative needs orbital response#

Differentiate the energy once more and something changes. The gradient — the first derivative — needed no orbital-response term because Hellmann–Feynman collapses it away: the SCF is stationary in the orbitals (\(\mathbf g=0\), SCF convergence theory), so \(\partial E/\partial R_A\) only ever sees the operators’ explicit \(R\)-dependence. Differentiating that stationarity condition a second time is where the orbitals’ own response reappears — because how the orbitals relax as one nucleus moves is itself a function of where every other nucleus is:

\[ \frac{\partial^2 E}{\partial R_A\,\partial R_B} = \underbrace{E^{[2]}_{AB}}_{\text{skeleton: 2nd-derivative integrals}\;\times\;D,\,W} \;+\; \underbrace{\sum_{ai} \frac{\partial F_{ai}}{\partial R_A}\,U^{B}_{ai} \;+\; \text{c.c.}}_{\text{orbital response — needs }U^{B}}. \]

The skeleton term is the easy half — second-derivative integrals (\(\partial^2 T\), \(\partial^2 V_{ne}\), \(\partial^2(\mu\nu|\lambda\sigma)\), \(\partial^2 S\) via \(W\)) contracted with the unperturbed density, exactly like the gradient but one derivative order up. The orbital-response term is new physics: \(U^B = \partial\mathbf C/\partial R_B\), how the orbitals themselves relax as nucleus \(B\) moves, found by solving the coupled-perturbed Hartree–Fock/Kohn–Sham (CPHF/CPKS) equations

\[ (\mathbf A+\mathbf B)\,U^B = -\,\mathbf b^B, \]

the same orbital-Hessian operator \(\mathbf A+\mathbf B\) from the SCF convergence theory augmented-Hessian machinery — reused here as a linear solve rather than an eigenvalue problem, with a geometry-derivative right-hand side \(\mathbf b^B\) built from the first-derivative Fock \(F^{B,x}\). This is the one genuinely new ingredient beyond the gradient: qc-rs’s CPHF/CPKS response engine already exists (it also powers stability analysis), so the Hessian assembles the “skeleton + fold-back” recipe on top of it:

        flowchart LR
    S["Skeleton term<br/>2nd-deriv integrals × D, W"] --> SUM["+"]
    G["1st-derivative Fock F^B,x<br/>(already built for the gradient)"] --> RHS["CPHF right-hand side b^B"]
    RHS --> CPHF["Solve (A+B) U^B = −b^B<br/>(reuses the SCF response engine)"]
    CPHF --> FOLD["Fold-back term<br/>Σ ∂F/∂R_A · U^B"]
    FOLD --> SUM
    SUM --> H["Molecular Hessian ∂²E/∂R_A∂R_B"]

Algorithm 2 (Assembling the molecular Hessian)

Input: converged SCF (\(\mathbf C\), \(\varepsilon\), \(\mathbf D\)), the CPHF response engine. Output: the Hessian \(\mathbf H \in \mathbb R^{3n_{\text{atom}}\times 3n_{\text{atom}}}\).

  1. Accumulate the skeleton blocks — one-electron, two-electron (4-center or RI), and (for KS) XC second-derivative terms — each contracted with the fixed density \(\mathbf D\) and energy-weighted density \(\mathbf W\).

  2. Build the geometry-derivative Fock \(F^{B,x}\) and overlap \(S^{B,x}\) for every nucleus \(B\) (reusing the same first-derivative integrals the gradient already needs).

  3. For each perturbation \(B\), solve \((\mathbf A+\mathbf B)U^B = -\mathbf b^B\) for the orbital response \(U^B\) (the CPHF/CPKS solve — one linear solve per nucleus, \(3n_{\text{atom}}\) total).

  4. Fold \(U^B\) back: add \(\sum_{ai}(\partial F_{ai}/\partial R_A)\,U^B_{ai} + \text{c.c.}\) to the skeleton.

  5. Add the classical nuclear-repulsion Hessian \(\partial^2 V_{nn}/\partial R_A\partial R_B\) (closed form).

  6. Self-test: the acoustic sum rule \(\sum_B \mathbf H_{AB} \approx \mathbf 0\) must hold (a rigid translation of the whole molecule changes no second derivative either) — the Hessian analogue of the gradient’s translational-invariance check.

From the Hessian to a spectrum: normal modes#

A Hessian in Cartesian coordinates mixes translation, rotation, and genuine vibration. Untangling them is a short, standard recipe:

  1. Mass-weight the Hessian, \(\tilde H_{Ai,Bj} = H_{Ai,Bj}/\sqrt{m_A m_B}\) (masses in atomic units).

  2. Project out rigid translation and rotation — 6 directions for a nonlinear molecule, 5 for a linear one (rotation about the molecular axis carries no energy) — using the mass-weighted translation/rotation projector.

  3. Diagonalize the projected \(\tilde H\): each eigenvalue \(k\) is a normal-mode force constant, and the corresponding eigenvector is the normal mode — the pattern of atomic displacement for that vibration.

  4. Convert to a frequency: \(\tilde\nu\,[\text{cm}^{-1}] = \sqrt{k}\cdot\text{(unit factor)}\). A negative \(k\) gives an imaginary frequency — the signature of a saddle point, not a true minimum (one imaginary mode along the reaction coordinate is exactly what a transition state looks like).

A nonlinear molecule with \(N\) atoms therefore has \(3N-6\) vibrational modes; a linear one has \(3N-5\) (one fewer rotation to remove). Water (\(N=3\), nonlinear) has 3; H₂ (\(N=2\), linear) has 1 — both verified below.

Thermochemistry: from frequencies to free energy#

The harmonic frequencies feed standard ideal-gas statistical mechanics (translational + rotational + vibrational partition functions) to produce the quantities a chemist actually wants at a given temperature and pressure:

  • Zero-point energy \(\text{ZPE} = \tfrac12\sum_k h c\,\tilde\nu_k\) — the vibrational ground-state energy every harmonic mode contributes even at \(T=0\).

  • Enthalpy \(H\), entropy \(S\), Gibbs free energy \(G=H-TS\), and the heat capacities \(C_v,C_p\) — each built from the translational, rotational, and vibrational partition functions in the usual way, added to the electronic energy.

\(G\) is what actually predicts chemical equilibria and, via transition-state theory, reaction rates — the harmonic frequencies are the bridge from “one SCF energy at one geometry” to that.

Usage#

The Hessian is a result on the SCF step, just like the gradient:

import qc, numpy as np
water = "O 0 0 0.117; H 0 0.757 -0.469; H 0 -0.757 -0.469"
m = qc.chk.new(atom=water, ao="sto-3g", unit="angstrom").scf(ref="r").run()

H = np.asarray(m.scf.hessian)      # shape (3*natom, 3*natom) = (9, 9) for water
H.sum(axis=1)                       # ~0 in every row: the acoustic sum-rule self-test

qc.thermo.frequencies(mychk) takes it the rest of the way — mass-weighting, projection, diagonalization, and the thermochemistry — in one call:

fr = qc.thermo.frequencies(m)

fr.frequencies     # array of harmonic wavenumbers, cm⁻¹
fr.n_imaginary     # count of imaginary (negative-k) modes: 0 at a true minimum
fr.norm_mode       # the normal-mode displacement vectors
fr.thermo          # dict: temperature, pressure, zpe, e_tot, h_tot, g_tot, s_tot, cv_tot, cp_tot

A converged geometry is a prerequisite in spirit, not just in practice: frequencies computed away from a true minimum (\(\mathbf g\ne 0\)) are not meaningful harmonic frequencies — always run .opt() (gradients & geometry optimization) first for a real analysis.

Worked example#

import qc, numpy as np
water = "O 0 0 0.117; H 0 0.757 -0.469; H 0 -0.757 -0.469"
m = qc.chk.new(atom=water, ao="sto-3g", unit="angstrom").scf(ref="r").run()

fr = qc.thermo.frequencies(m)
print("frequencies (cm⁻¹):", np.round(np.sort(fr.frequencies), 2))
print("imaginary modes    :", fr.n_imaginary)
print("ZPE (Ha)           :", round(fr.thermo["zpe"], 6))
print("G at 298.15 K (Ha) :", round(fr.thermo["g_tot"], 6))

# frequencies (cm⁻¹): [2041.42 4494.04 4796.72]
# imaginary modes    : 0
# ZPE (Ha)           : 0.025817
# G at 298.15 K (Ha) : 0.033348

# a linear molecule: 3N-5 = 1 mode, not 3N-6
h2 = qc.chk.new(atom="H 0 0 0; H 0 0 0.74", ao="sto-3g", unit="angstrom").scf(ref="r").run()
print("H2 modes:", len(qc.thermo.frequencies(h2).frequencies))   # 1

Three real, positive frequencies and zero imaginary modes confirm water’s geometry is a genuine minimum, not a saddle point — and the count itself (3 for nonlinear water, 1 for linear H₂) is a direct check that the translation/rotation projection worked.

Exercise 7

  1. A colleague computes frequencies at a geometry they only partially optimized. n_imaginary comes back as

    1. What does that mean, and what should they check before trusting the frequencies?

  2. Why does H₂ (linear, \(N=2\)) have exactly \(3N-5=1\) vibrational mode while a bent triatomic like water (\(N=3\), nonlinear) has \(3N-6=3\) — where did the “missing” rotational mode for water go, relative to a linear molecule with the same atom count?

  3. The Hessian needs a CPHF solve that the gradient does not. In one sentence, why does the first derivative avoid it while the second cannot?

With energies, forces, and now vibrational/thermochemical analysis in hand, the remaining guide chapters turn to the environment — solvation & dispersion — and then visualization, logging, and the full properties suite.