# The analytic Hessian, vibrational frequencies & thermochemistry

[The gradient chapter](analytic-derivatives.md) showed that a first energy derivative needs no orbital
response, only the converged density and $W$. The second derivative loses that simplification — it needs
genuine coupled-perturbed orbital response, computed via [the CPHF/CPKS machinery](linear-response-cphf.md)
— and once it is in hand, a short, standard recipe converts it into the two things chemists actually
measure at a stationary point: vibrational frequencies and thermochemistry. This chapter derives the full
skeleton-plus-fold-back Hessian assembly, the mass-weighted translation/rotation projection that isolates
genuine vibration, and the exact ideal-gas partition-function formulas qc-rs evaluates — grounded in
`.design/76analytic-hessian-frequencies.md` and `python/qc/thermo.py` (a direct port of PySCF's
`hessian/thermo.py`). [The guide chapter](../20-guide/hessian-frequencies-thermo.md) carries the compact
usage digest; this chapter carries the full derivation.

## Why the second derivative reintroduces orbital response

Differentiating the stationarity condition $\partial E/\partial C=0$ **a second time** is exactly where
orbital response reappears — how the orbitals relax as *one* nucleus moves is itself a function of where
every *other* nucleus sits, and that dependence was invisible to the first derivative precisely because
Hellmann-Feynman made it vanish there. The full second derivative splits cleanly into two pieces:

$$
\frac{\partial^2E}{\partial R_A\partial R_B} = \underbrace{E^{[2]}_{AB}}_{\text{skeleton (fixed }C\text{)}}
+ \underbrace{\sum_{ai}\frac{\partial F_{ai}}{\partial R_A}\,U^B_{ai} + \text{c.c.}}_{\text{orbital response
— needs }U^B},
$$

where $U^B=\partial\mathbf C/\partial R_B$ is found by solving the coupled-perturbed equations
$(\mathbf A+\mathbf B)U^B=-\mathbf b^B$ — literally the same $A+B$ operator from [the linear-response
chapter](linear-response-cphf.md), reused here as a linear solve (one right-hand side per nuclear Cartesian
coordinate, $3n_{\text{atom}}$ solves total) rather than an eigenvalue problem.

## The skeleton term in full

The skeleton is the "easy half" — every term contracted with the *unperturbed* density and energy-weighted
density, exactly as in the gradient but one derivative order higher, using **three** distinct
second-derivative two-electron integral classes rather than the gradient's single `int2e_ip1`:

$$
E^{J,[2]}_{AB} - c_xE^{K,[2]}_{AB} = \tfrac12\sum_{\mu\nu\lambda\sigma}\Bigl(D_{\mu\nu}D_{\lambda\sigma} -
\tfrac{c_x}{2}D_{\mu\lambda}D_{\nu\sigma}\Bigr)\,\partial^2_{AB}(\mu\nu|\lambda\sigma),
$$

where the second derivative $\partial^2_{AB}$ falls into three block patterns depending on which pair of
atoms it differentiates:

| Block | Integral | What it differentiates |
|---|---|---|
| Diagonal (`A,A`) | `int2e_ipip1` | both derivatives on the same bra center |
| Off-diagonal (`A,B`) | `int2e_ip1ip2` | one derivative on the bra center, one on the ket center |
| Within-pair (`A,B`) | `int2e_ipvip1` | both derivatives on the bra pair, but on different bra shells |

alongside the one-electron/Pulay skeleton (second-derivative kinetic/nuclear-attraction integrals
contracted with $D$, second-derivative overlap contracted with $W$ — structurally the gradient's Pulay
term one order up) and, for KS references, the XC second-derivative contribution built from the same
grid machinery [derived in the DFT-grid chapter](dft-xc-quadrature.md).

## The CPHF right-hand side and the fold-back

The right-hand side driving each nucleus's coupled-perturbed solve reuses the *gradient's own*
first-derivative two-electron integral (`int2e_ip1`, [already built for the gradient](analytic-derivatives.md))
— nothing new needs to be computed for this piece, only recombined:

$$
F^{A,x}_{\mu\nu} = h^{A,x}_{\mu\nu} + \sum_{\lambda\sigma}D_{\lambda\sigma}\Bigl[(\mu\nu|\lambda\sigma)^{A,x}
- \tfrac{c_x}{2}(\mu\lambda|\nu\sigma)^{A,x}\Bigr],
$$

with $(\cdot)^{A,x}$ meaning the bra differentiated on atom $A$'s shells, and $h^{A,x}$ the first-derivative
core Hamiltonian from the gradient chapter. A moving AO basis means the overlap derivative
$S^{A,x}\ne0$, which fixes the occupied-occupied block of the response by the orthonormality constraint
itself — $U^{A,x}_{ij}=-\tfrac12S^{A,x}_{ij}$ with $S^{A,x}_{ij}=(C^{\mathsf T}S^{A,x}C)_{ij}$ — before the
virtual-occupied block is even solved. The occupied-occupied block's already-known response then folds
into the virtual-occupied CPHF right-hand side:

$$
\sum_{bj}(A+B)_{ai,bj}\,U^{A,x}_{bj} = -b^{A,x}_{ai}, \qquad
b^{A,x}_{ai} = F^{A,x}_{ai} - S^{A,x}_{ai}\varepsilon_i - \sum_{kl}(A+B)_{ai,kl}\Bigl(-\tfrac12S^{A,x}_{kl}\Bigr),
$$

which is exactly `linear_solve(apply_aplusb, rhs=-b, ...)` — [the same CPHF solver](linear-response-cphf.md)
stability analysis and QC-SCF already share, just with this geometry-specific right-hand side. Once $U^{A,x}$
converges, the orbital-energy response follows as a byproduct,
$\varepsilon^{A,x}_{ij}=b^{A,x}_{ij}\big|_{oo}+[(A+B)U^{A,x}]_{ij}+U^{A,x}_{ij}(\varepsilon_i-\varepsilon_j)$,
needed for the fold-back:

$$
E^{\text{resp}}_{AB} = 4\sum_{\mu\nu}F^{A,x}_{\mu\nu}\bigl(U^{B,y}C_{\text{occ}}^{\mathsf T}\bigr)_{\mu\nu}
- 4\sum_{\mu\nu}W^{B,y}_{\mu\nu}S^{A,x}_{\mu\nu} - 2\sum_{ij}\bigl(C^{\mathsf T}S^{A,x}C\bigr)_{ij}\,
\varepsilon^{B,y}_{ij}.
$$

**ROHF/ROKS reuses the exact same energy-weighted-density projector construction** [the gradient
chapter](analytic-derivatives.md) derived for $W$ ($W=\sum_\sigma P_\sigma F_\sigma P_\sigma$ from the spin
Fock matrices, never the effective-Fock eigenvalues) — the design note flags this as "the same subtlety
[the gradient chapter] flagged, now also driving the CPHF right-hand side," i.e. getting $W$ wrong here
breaks not only the skeleton term but the very right-hand side the coupled-perturbed solve depends on. One
error in $W$, two places it can silently corrupt.

:::{prf:algorithm} Assembling the molecular Hessian
:label: alg-hessian-full

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

1. Accumulate the skeleton blocks (one-electron + Pulay, the three two-electron second-derivative
   classes, and XC if present) contracted with the fixed $\mathbf D$ and $\mathbf W$.
2. Build $F^{A,x}$ and $S^{A,x}$ for every nucleus $A$ and Cartesian direction $x$ (reusing the gradient's
   own `int2e_ip1`/`int1e_ip*` integrals).
3. For each of the $3n_{\text{atom}}$ perturbations: fix the occupied-occupied block
   $U^{A,x}_{ij}=-\tfrac12S^{A,x}_{ij}$, fold it into the CPHF right-hand side $b^{A,x}$, and solve
   $(A+B)U^{A,x}_{\text{vo}}=-b^{A,x}_{\text{vo}}$ via the shared `linear_solve`.
4. Compute the orbital-energy response $\varepsilon^{A,x}$ as a byproduct of the converged $U^{A,x}$.
5. Fold every pair $(A,B)$ back into the Hessian via $E^{\text{resp}}_{AB}$ above.
6. Add the classical nuclear-repulsion Hessian $\partial^2V_{nn}/\partial R_A\partial R_B$ (closed form).
7. **Self-test — the acoustic sum rule**: $\sum_B\mathbf H_{AB}\approx\mathbf0$, the Hessian analogue of
   the gradient's translational-invariance check (a rigid translation of the whole molecule changes no
   second derivative either).
:::

## From the Hessian to normal modes: the projection in full

A Cartesian Hessian mixes translation, rotation, and genuine vibration in a single matrix — recovering the
vibrations alone needs an explicit projection, not just a diagonalization. qc-rs's `harmonic_analysis`
(`python/qc/thermo.py`, a direct port of PySCF's `hessian/thermo.py`) constructs this projection from first
principles rather than a canned formula:

1. **Mass-weight**: $\tilde H_{Ai,Bj}=H_{Ai,Bj}/\sqrt{m_Am_B}$ (masses in atomic units), symmetrized.
2. **Build the rigid translation/rotation basis directly in mass-weighted coordinates.** Three
   translation directions are simply $\sqrt{m_A}\,\hat{\mathbf e}_k$ for each Cartesian axis $k$; three
   rotation directions are $\sqrt{m_A}\,(\hat{\mathbf e}_k\times\mathbf r_A^{\text{rel}})$ about each axis,
   with $\mathbf r_A^{\text{rel}}$ the atom's displacement from the center of mass. A QR decomposition of
   these six columns orthonormalizes them and **automatically drops a null column for a linear
   molecule** (rotation about the molecular axis moves no mass, so that column is exactly zero and QR
   detects it via a vanishing diagonal element) — this is *why* qc-rs never needs a special-cased
   "is this molecule linear?" branch to get $3N-5$ vs. $3N-6$ right; the projection handles it
   automatically from the geometry itself.
3. **Project onto the orthogonal complement**: $P_{\text{vib}}=I-QQ^{\mathsf T}$ with $Q$ the
   translation/rotation basis; diagonalizing $P_{\text{vib}}$ itself and keeping eigenvectors with
   eigenvalue $>0.5$ gives an orthonormal basis for the vibrational subspace alone (dimension $3N-6$ or
   $3N-5$).
4. **Diagonalize the mass-weighted Hessian restricted to this vibrational subspace.** Each eigenvalue $k$
   is a force constant; $\tilde\nu=\sqrt{|k|}\cdot(\text{unit factor})$ gives the frequency in
   wavenumbers, with the *sign* of $k$ preserved as the sign of the reported frequency — a **negative**
   force constant gives a **negative** (i.e. imaginary) reported frequency, the unambiguous signature of a
   saddle point rather than a true minimum.
5. **Un-mass-weight** the surviving eigenvectors back to genuine Cartesian atomic displacements (dividing
   by $\sqrt{m_A}$ again) to get the normal modes actually visualized, and compute each mode's **reduced
   mass** as $\mu_k=1/\sum_{A}|\mathbf d_{A,k}|^2$ from the un-normalized displacement.

## IR intensities: reusing the normal modes directly

Given the converged normal modes and the **atomic polar tensor** (APT, $\partial\mu_i/\partial R_{A,x}$ —
the dipole-moment derivative with respect to nuclear displacement, available for RHF/UHF/RKS/UKS and ROHF
references), the double-harmonic-approximation IR intensity per mode is a direct projection of the APT
onto each mass-weighted normal mode:

$$
\frac{\partial\mu_i}{\partial Q_k} = \sum_{A,x}\frac{\partial\mu_i}{\partial R_{A,x}}\,
\frac{L_{A,x,k}}{\sqrt{m_A}}, \qquad I_k\,[\text{km/mol}] = 974.8801\times4.8032^2\sum_i\Bigl(
\frac{\partial\mu_i}{\partial Q_k}\Bigr)^2,
$$

with $L$ the orthonormal mass-weighted normal-mode matrix from the projection above and the constant
converting atomic units (electron charge per bohr, mass-weighted in amu$^{-1/2}$) into the conventional
km/mol intensity unit via the Debye/Å conversion factor. No new machinery is needed beyond what the
frequency analysis already built — the APT is itself a derivative-integral object closely related to the
gradient's own dipole and derivative-Fock integrals, and the intensity is a pure linear-algebra contraction
once both the APT and the normal modes exist.

## Ideal-gas thermochemistry: the exact partition functions

The harmonic frequencies feed standard ideal-gas statistical mechanics. qc-rs evaluates each partition
function in closed form rather than looking anything up in a table, using CODATA physical constants
throughout ($k_B$, $h$, $N_A$, the Hartree-to-Joule conversion, …):

**Vibrational** (harmonic oscillator, one factor per real positive mode, with vibrational temperature
$\Theta_k=hc\tilde\nu_k/k_B$ and $x_k=\Theta_k/T$):

$$
\text{ZPE} = \tfrac12\sum_k k_B\Theta_k, \qquad
E_{\text{vib}} = \sum_k k_B\Theta_k\Bigl(\tfrac12+\frac{1}{e^{x_k}-1}\Bigr), \qquad
S_{\text{vib}} = k_B\sum_k\Bigl(\frac{x_k}{e^{x_k}-1}-\ln(1-e^{-x_k})\Bigr).
$$

Only real, positive frequencies enter here — an imaginary mode (from an incompletely optimized geometry)
is excluded from the vibrational sum entirely rather than producing a nonsensical contribution, which is
also why `n_imaginary>0` should be treated as "these thermochemistry numbers are not meaningful yet," not
merely as a warning to note in passing.

**Translational** (an ideal gas of total mass $M$, at temperature $T$ and pressure $P$; this is exactly the
Sackur-Tetrode construction):

$$
q_{\text{trans}} = \Bigl(\frac{2\pi Mk_BT}{h^2}\Bigr)^{3/2}\frac{k_BT}{P}, \qquad
E_{\text{trans}} = \tfrac32k_BT, \qquad S_{\text{trans}} = k_B\bigl(\ln q_{\text{trans}}+\tfrac52\bigr).
$$

**Rotational**, from the principal moments of inertia $I_1,I_2,I_3$ (eigenvalues of the inertia tensor
built from the same nuclear masses/geometry — [the geometric-analysis chapter](../20-guide/properties/geometric-analysis.md)
derives the same inertia tensor for rotational-constant reporting) and rotational temperature
$\Theta_{r}=\hbar^2/(2Ik_B)$ per nonzero moment — a **linear** molecule has exactly one vanishing moment
(no energy cost to spin about the molecular axis), collapsing the general nonlinear formula to a simpler
one:

$$
\text{linear:}\quad q_{\text{rot}}=\frac{T}{\sigma\Theta_r}, \quad E_{\text{rot}}=k_BT, \quad
S_{\text{rot}}=k_B\bigl(\ln q_{\text{rot}}+1\bigr);
$$
$$
\text{nonlinear:}\quad q_{\text{rot}}=\frac{\sqrt\pi}{\sigma}\sqrt{\frac{T^3}{\Theta_{r,1}\Theta_{r,2}
\Theta_{r,3}}}, \quad E_{\text{rot}}=\tfrac32k_BT, \quad S_{\text{rot}}=k_B\bigl(\ln
q_{\text{rot}}+\tfrac32\bigr),
$$

with $\sigma$ the rotational symmetry number (qc-rs currently fixes $\sigma=1$ — a point-group-aware
symmetry-number refinement is a known follow-up, not yet implemented; for an asymmetric molecule like
water this default is exact, but for a symmetric one it would need correcting). Whether a molecule counts
as linear is decided directly from the inertia eigenvalues (one vanishingly small compared to the others),
the same automatic detection the normal-mode projection above uses — no separate "is it linear" flag is
threaded through from user input.

**Assembly** (ideal-gas relations $H=U+k_BT$, $G=H-TS$, an electronic spin-multiplicity entropy term
$S_{\text{elec}}=k_B\ln(\text{mult})$, and $C_p=C_v+k_B$ from the ideal-gas equation of state):

$$
U = E_{\text{trans}}+E_{\text{rot}}+E_{\text{vib}}+\text{ZPE}, \qquad H=U+k_BT, \qquad
S = S_{\text{trans}}+S_{\text{rot}}+S_{\text{vib}}+S_{\text{elec}}, \qquad G=H-TS.
$$

Verified example — water/STO-3G, the full harmonic + thermochemistry + IR-intensity pipeline:

```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)
np.round(np.sort(fr.frequencies), 2)   # [2041.42, 4494.04, 4796.72]  cm^-1
fr.n_imaginary                          # 0 -- a genuine minimum
fr.intensities                          # [263.64, 861.86, 459.39]  km/mol (IR intensities)
fr.thermo["zpe"], fr.thermo["g_tot"]     # (0.025817, 0.033348)  Hartree at 298.15 K, 101325 Pa
```

Three real positive frequencies, zero imaginary modes, and a set of physically reasonable IR intensities
together confirm this is a genuine minimum with a fully self-consistent harmonic analysis — and the
acoustic sum rule ($\sum_B H_{AB}\approx0$, [verified directly in the guide chapter](../20-guide/hessian-frequencies-thermo.md))
is the Hessian-level counterpart of the gradient's own translational-invariance check.

:::{exercise}
:label: ex-hessian-thermo-theory

1. The translation/rotation basis construction detects a linear molecule by a QR decomposition producing
   a vanishing diagonal element for the axial-rotation column, rather than a special-cased check on atom
   count or geometry. Explain in one sentence why $\sqrt{m_A}(\hat{\mathbf e}_\parallel\times\mathbf
   r_A^{\text{rel}})$ is exactly zero for every atom when $\hat{\mathbf e}_\parallel$ is the molecular
   axis of a linear molecule, and why that is precisely what a linear molecule "having one fewer
   rotational degree of freedom" means physically.
2. The vibrational entropy/energy sums explicitly exclude negative (imaginary) frequencies rather than,
   say, taking $|k|$ and including them anyway. Why would including an imaginary mode in the vibrational
   partition function produce numbers that are not merely inaccurate but physically meaningless?
3. The design note flags the ROHF/ROKS energy-weighted density $W$ as "one error, two places" for the
   Hessian specifically (beyond the gradient). Identify the second place, and explain briefly why getting
   $W$ wrong there corrupts more than just one energy term.
:::

:::{solution} ex-hessian-thermo-theory
:class: dropdown

1. Rotating about the molecular axis of a linear molecule leaves every atom exactly where it was (each
   atom lies *on* that axis, so $\mathbf r_A^{\text{rel}}$ is parallel to $\hat{\mathbf e}_\parallel$, and
   the cross product of two parallel vectors is identically zero) — there is no atomic displacement at
   all associated with that "rotation," for any atom, at any mass. A rigid-body motion that displaces
   nothing is not a genuine degree of freedom, which is exactly what "one fewer rotational degree of
   freedom for a linear molecule" means: the QR decomposition's vanishing diagonal element is simply
   detecting, numerically, that this basis vector was already the zero vector.
2. The harmonic-oscillator partition function and its derived entropy/energy formulas are built from
   $\Theta_k=hc\tilde\nu_k/k_B$ appearing inside an $e^{-\Theta_k/T}$-type Boltzmann factor — these
   formulas assume a genuine, positive-curvature harmonic well with a well-defined ground state and
   equally spaced excited levels. A negative force constant does not correspond to a bound harmonic
   oscillator at all — it describes motion down an unstable direction (the saddle-point reaction
   coordinate), for which there is no well-defined vibrational ground state or level spacing to plug into
   a partition function. Taking $|k|$ and computing anyway would produce a number that looks like an
   entropy or energy contribution but corresponds to no physical vibrational mode whatsoever — worse than
   omitting it, because it would silently contaminate $G$ and $S$ with a fabricated contribution instead
   of visibly signaling (via `n_imaginary>0`) that the geometry is not a genuine minimum.
3. The first place is the skeleton Hessian term itself (exactly as in the gradient, $W$ enters the Pulay
   overlap-derivative contraction). The second, Hessian-specific place is the **CPHF right-hand side**
   $b^{A,x}$ and the fold-back term $E^{\text{resp}}_{AB}$, both of which explicitly contain $W^{B,y}$
   contracted against $S^{A,x}$. Getting $W$ wrong therefore does not just bias one skeleton contribution
   — it corrupts the very perturbation that drives every one of the $3n_{\text{atom}}$ coupled-perturbed
   solves, so the error propagates into the orbital response $U^{A,x}$ itself and hence into every
   fold-back term built from it, not merely into a single additive piece of the final Hessian.
:::

With the SCF energy, its gradient, its Hessian, and now a complete vibrational/thermochemical analysis
derived from first principles, [the ECP chapter](ecp-theory.md) and [the solvation
chapter](solvation-theory.md) turn to the remaining Hamiltonian modifications this manual's Part II has not
yet covered — effective core potentials and implicit solvation — both of which the analytic-gradient and
(where implemented) analytic-Hessian machinery this chapter derived already extends to cover.