# Hessians, vibrational frequencies & thermochemistry
The [gradient](gradients-geomopt.md) 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](../10-foundations/scf-convergence-theory.md)), 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](../10-foundations/scf-convergence-theory.md) 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](scf.md)), so the Hessian assembles the "skeleton + fold-back" recipe on top of it:
```{mermaid}
flowchart LR
S["Skeleton term
2nd-deriv integrals × D, W"] --> SUM["+"]
G["1st-derivative Fock F^B,x
(already built for the gradient)"] --> RHS["CPHF right-hand side b^B"]
RHS --> CPHF["Solve (A+B) U^B = −b^B
(reuses the SCF response engine)"]
CPHF --> FOLD["Fold-back term
Σ ∂F/∂R_A · U^B"]
FOLD --> SUM
SUM --> H["Molecular Hessian ∂²E/∂R_A∂R_B"]
```
:::{prf:algorithm} Assembling the molecular Hessian
:label: alg-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:
```python
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:
```python
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](gradients-geomopt.md)) first for a real analysis.
## Worked example
```python
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}
:label: ex-hessian
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?
:::
:::{solution} ex-hessian
:class: dropdown
1. One imaginary frequency means the Hessian has a **negative eigenvalue** — the structure is a **saddle
point** (e.g. a transition state), not a minimum, consistent with an incomplete optimization. They should
re-run `.opt()` to full convergence (checking `opt.converged` and the gradient trajectory) before trusting
any frequency as a real vibration.
2. A linear molecule has only **2** independent rotational degrees of freedom (rotation about its own axis
carries no energy, since the moment of inertia about that axis is zero), so only 5 of the 6 rigid-body
directions are removed, leaving $3N-5$ vibrations. A nonlinear molecule has all 3 rotations, removing 6
and leaving $3N-6$. Water, being *bent* rather than linear, has the full 3 rotational degrees of freedom
to remove; the mode a linear triatomic would have "extra" relative to water is exactly the bending motion
that in water is already counted among its 3 vibrations, not a missing one.
3. Hellmann–Feynman makes the *first* derivative stationary-in-the-orbitals, so orbital response cancels;
differentiating that stationarity condition a second time is precisely what reintroduces the orbitals'
own response to a *second* nuclear displacement, which only a CPHF solve can supply.
:::
With energies, forces, and now vibrational/thermochemical analysis in hand, the remaining guide chapters turn
to the **environment** — [solvation & dispersion](solvation-dispersion.md) — and then visualization, logging,
and the full properties suite.