The analytic Hessian, vibrational frequencies & thermochemistry

The gradient chapter 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 — 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 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, 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.

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) — 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 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 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.

Algorithm 5 (Assembling the molecular Hessian)

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 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:

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) is the Hessian-level counterpart of the gradient’s own translational-invariance check.

Exercise 12

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 to Exercise 12

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 and the solvation chapter 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.