Analytic derivatives: the nuclear gradient

Analytic derivatives: the nuclear gradient#

Every geometry optimization, vibrational analysis, and ab initio molecular dynamics step needs the force on every nucleus, \(-\partial E/\partial X_A\). Computing it by finite difference costs \(6N_{\text{atom}}\) extra SCF energies (perturbing each Cartesian coordinate forward and backward) and is numerically noisy; computing it analytically costs roughly one extra pass over the same derivative integral classes the energy build already used, at machine precision. This chapter derives the full analytic RHF/UHF/ROHF gradient — why no coupled-perturbed solve is needed, the exact energy-weighted density for every SCF reference including the ROHF correctness subtlety, and the two-electron gradient’s quartet/RI structure — grounded in .design/74analytic-nuclear-gradients.md. The gradients/geometry-optimization chapter carries the compact usage digest; this chapter carries the derivation in full.

The total derivative and the key simplification#

For a variational SCF energy \(E(\{C\},\{\varepsilon\};X)\), the total derivative with respect to a nuclear Cartesian coordinate \(X_A\) has, in principle, three kinds of terms: how the Hamiltonian’s operators depend on \(X_A\) at fixed AO basis and density, how the AO basis functions themselves move (since each GTO is centered on a nucleus), and how the converged MO coefficients would respond to the geometry perturbation. The last of these — orbital response — is what makes gradients expensive in general (coupled-perturbed Hartree-Fock/Kohn-Sham, derived in full in the next chapter). But for the energy gradient specifically, it vanishes identically, by a direct consequence of the variational principle:

\[ \frac{\partial E}{\partial X_A}\bigg|_{\text{via }C} = \sum_{\mu\nu}\frac{\partial E}{\partial D_{\mu\nu}}\cdot\underbrace{\frac{\partial D_{\mu\nu}}{\partial C}\frac{\partial C}{\partial X_A}}_{\text{orbital response}} = 0, \]

because \(\partial E/\partial C=0\) at SCF convergence — this is the Hellmann-Feynman theorem applied to a variational energy functional. The catch is that HF/KS orbitals are not merely eigenfunctions of a fixed operator; they satisfy an orthonormality constraint \(C^{\mathsf T}SC=I\) that itself depends on the geometry (through the AO overlap \(S\)), and differentiating that constraint is what reintroduces a geometry-dependent term — the Pulay force — even though the orbital-rotation response drops out entirely. Collecting every surviving term (RHF form; UHF/ROHF/KS are per-spin/effective generalizations of the same skeleton):

\[ \frac{\partial E}{\partial X_A} = \underbrace{\sum_{\mu\nu}D_{\mu\nu}\frac{\partial h_{\mu\nu}}{\partial X_A}}_{\text{core-Hamiltonian}} + \underbrace{\tfrac12\sum_{\mu\nu\lambda\sigma}D_{\mu\nu}D_{\lambda\sigma} \frac{\partial(\mu\nu|\lambda\sigma)}{\partial X_A}\bigg|_{JK}}_{\text{two-electron }(J-c_xK)} - \underbrace{\sum_{\mu\nu}W_{\mu\nu}\frac{\partial S_{\mu\nu}}{\partial X_A}}_{\text{Pulay}} + \frac{\partial E_{\text{xc}}}{\partial X_A} + \frac{\partial E_{\text{ECP}}}{\partial X_A} + \frac{\partial E_{\text{PCM}}}{\partial X_A} + \frac{\partial V_{nn}}{\partial X_A}. \]

The practical consequence: a converged HF or pure/hybrid KS gradient needs only the converged density \(D\) and one further object, the energy-weighted density \(W\) — no iterative coupled-perturbed solve at all. (The Z-vector/CPHF machinery is needed only for gradients of non-variational quantities — MP2 and other post-SCF correlation energies — which is genuinely out of scope for the plain SCF gradient this chapter covers, and is instead where linear-response theory becomes essential.)

The energy-weighted (Lagrangian) density \(W\)#

\(W\) is exactly the object the Pulay term needs to absorb the orbital-orthonormality constraint’s geometry dependence — it plays the role, in the overlap term, that \(D\) plays in the Hamiltonian term.

RHF: \(W = C_{\text{occ}}^{\mathsf T}\operatorname{diag}(2\varepsilon_{\text{occ}})\,C_{\text{occ}}\) — structurally identical to building \(D\) from occupied MOs, but weighting each orbital by its orbital energy instead of by its occupation number.

UHF: \(W = \sum_\sigma C_{\text{occ}}^{\sigma\mathsf T}\operatorname{diag}(\varepsilon^\sigma_{\text{occ}})\, C_{\text{occ}}^\sigma\) — the same construction, independently per spin channel.

ROHF/ROKS — the correctness landmine. ROHF diagonalizes one effective Roothaan Fock operator shared by both spins (the same fact that made ROHF-MP2 ambiguous), so the eigenvalues stored as orbital_energies are eigenvalues of that effective operator, not genuine Lagrangian multipliers for the true spin-resolved constraint — naively reusing the RHF/UHF formula with these eigenvalues gives a plausible-looking but wrong gradient (correct only in the trivial case of zero open shells). The correct construction instead rebuilds \(W\) directly from the converged spin Fock matrices \(F_\alpha,F_\beta\) (already available from the converged ROHF density) and per-spin AO projectors:

\[ W = W_\alpha + W_\beta, \qquad W_\sigma = P_\sigma\,F_\sigma\,P_\sigma, \qquad P_\alpha = C_{\text{occ}}^{\alpha\mathsf T}C_{\text{occ}}^\alpha,\quad P_\beta = C_{\text{occ}}^{\beta\mathsf T}C_{\text{occ}}^\beta, \]

with \(P_\alpha\) projecting onto the \(\alpha\)-occupied space (singly and doubly occupied orbitals) and \(P_\beta\) onto the \(\beta\)-occupied space (doubly occupied only) — this is the canonical construction (matching PySCF’s grad/rohf.py:make_rdm1e, not a simpler MO-basis block rule an earlier design draft had guessed and then had to correct). This is, by qc-rs’s own design note, “the single subtlest piece of correctness in the whole gradient” — validated by finite-difference comparison before being trusted, since an error here produces a gradient that looks reasonable (right order of magnitude, right qualitative direction) without being numerically correct.

The one-electron (core-Hamiltonian) and Pulay gradient#

\(h=T+V_{ne}\), so \(\partial h/\partial X = \partial T/\partial X + \partial V_{ne}/\partial X\). The kinetic term is purely an AO-basis-center derivative (the kinetic operator itself carries no nuclear coordinate). The nuclear attraction term splits into two physically distinct pieces: the AO-basis-center part (how the bra/ket Gaussian moves) and the genuinely Hellmann-Feynman part — the operator \(1/|\mathbf r-\mathbf R_B|\) itself depends on nucleus \(B\)’s position, so differentiating it with respect to \(R_B\) gives a term centered on the nucleus, scaled by that nucleus’s effective charge \(Z_{\text{eff},B}\):

\[ \sum_{\mu\nu}D_{\mu\nu}\frac{\partial h_{\mu\nu}}{\partial X_A} = \underbrace{\sum_{\mu\nu\in A}D_{\mu\nu} \left(\frac{\partial T_{\mu\nu}}{\partial X_A}+\frac{\partial V_{ne,\mu\nu}^{\text{basis}}}{\partial X_A}\right)}_{\text{AO-center, scattered to shells on atom }A} + \underbrace{Z_{\text{eff},A}\sum_{\mu\nu}D_{\mu\nu}\frac{\partial}{\partial R_A}\frac{1}{|\mathbf r-\mathbf R_A|}\bigg|_{\mu\nu}}_{\text{Hellmann-Feynman, one term per nucleus}}. \]

The Pulay term is structurally simple once \(W\) is in hand — it is just \(-\sum_{\mu\nu}W_{\mu\nu}\, \partial S_{\mu\nu}/\partial X_A\), an AO-center overlap derivative contracted against \(W\) exactly the way the Hamiltonian terms contract against \(D\). Every AO-center contribution here is naturally organized as a loop over shell pairs, scattering each pair’s contribution to whichever atom owns the differentiated shell; the Hellmann-Feynman piece is instead a loop over nuclei, since the operator itself sits on the nucleus rather than on a basis function. Adding the classical nuclear-repulsion gradient \(\partial V_{nn}/\partial X_A\) (an elementary Coulomb’s-law derivative between point charges) completes the one-electron side.

The two-electron gradient: why it has no Hellmann-Feynman term#

In chemist’s notation, \((\mu\nu|\lambda\sigma)=\iint\chi_\mu(\mathbf r_1)\chi_\nu(\mathbf r_1)\,r_{12}^{-1}\,\chi_\lambda(\mathbf r_2)\chi_\sigma(\mathbf r_2)\,d\mathbf r_1\,d\mathbf r_2\). The operator \(r_{12}^{-1}\) carries no nuclear coordinate at all — unlike \(V_{ne}\)’s \(1/|\mathbf r-\mathbf R_B|\), there is no nucleus for this operator to sit on. Consequently the two-electron gradient depends on geometry only through the four AO centers of \(\mu,\nu,\lambda,\sigma\), and has no Hellmann-Feynman term whatsoever — every contribution is a basis-center (Pulay-type) derivative. This is the structural reason the two-electron gradient needs only basis-center derivative integrals and never an integral centered on a bare nucleus.

For the exact 4-center path, the derivative is taken on the first bra shell (libcint’s int2e_ip1), and because a GTO is \(\chi(\mathbf r-\mathbf A)\), \(\partial\chi/\partial\mathbf A=-\nabla\chi\) — an electron-coordinate gradient with a sign flip. The quartet’s full derivative needs all four centers differentiated in turn, but permutational and translational symmetry recovers most of them for free: rotating which shell of the quartet plays the “first” (differentiated) role covers the remaining centers, and the whole-molecule identity \(\sum_A\mathbf F_A=0\) (translational invariance — the total force on a free molecule vanishes) lets the very last center’s contribution be inferred from the others rather than computed directly. The energy-side Coulomb/exchange pattern carries over unchanged into the gradient contraction — \(J\) contracts the total density against both index pairs, \(K\) contracts the spin density with one index crossed — with the same \(\tfrac12\), \(a_x\) (exact-exchange fraction), and range-separation weight scalars from the energy-side Fock build applied at assembly.

The RI two-electron gradient: an extra, mandatory piece#

The RI-factorized two-electron energy is \((\mu\nu|\lambda\sigma)_{\text{RI}}=\sum_{PQ}(\mu\nu|P)\, V^{-1}_{PQ}\,(Q|\lambda\sigma)\) (derived fully in the RI chapter). Differentiating this factorized form requires differentiating both the 3-center integral \((\mu\nu|P)\) and the 2-center metric \(V_{PQ}=(P|Q)\) — and because \(V^{-1}\) itself depends on geometry, its derivative must be accounted for through the chain rule. Defining the fitted charge \(\gamma_P=\sum_QV^{-1}_{PQ}c_Q\) with \(c_Q=\sum_{\mu\nu}(Q|\mu\nu)D_{\mu\nu}\) (exactly the resident whitened \(B\)-factor’s \(q_P\) from the RI chapter, no new object needed), the Coulomb energy is \(E_J=\tfrac12\gamma^{\mathsf T}V\gamma\) (using \(c=V\gamma\)), and differentiating it through the \(V^{-1}\)-chain-rule identity \(\partial(K^{-1}R)=K^{-1}(\partial R - \partial K\cdot q)\) (the same identity the PCM chapter uses for its own response derivative) gives a clean two-term result with no leftover explicit \(\partial V^{-1}/\partial X\):

\[ \frac{\partial E_J}{\partial X} = \sum_P\gamma_P\sum_{\mu\nu}\frac{\partial(P|\mu\nu)}{\partial X}D_{\mu\nu} \;-\; \tfrac12\sum_{PQ}\gamma_P\gamma_Q\,\frac{\partial(P|Q)}{\partial X}. \]

The analogous exchange-gradient expression needs one further fitted object, the doubly-MO-transformed metric factor \(G_{PQ}=\sum_{ij}\tilde d_{P,i,j}\tilde d_{Q,j,i}\) (deliberately not symmetric in \(i,j\) — one index carries the occupation weighting, the other must not), but the structure is the same: an AO-center term from \(\partial(P|\mu\nu)/\partial X\) and a metric term from \(\partial(P|Q)/\partial X\).

Three distinct derivative-integral classes and their atom-scatter targets:

Integral	Differentiates	Scatters to
`int3c2e_ip1`	\((P\|\mu\nu)\), AO centers	atom(\(\mu\)), atom(\(\nu\))
`int3c2e_ip2`	\((P\|\mu\nu)\), auxiliary center	atom(\(P\))
`int2c2e_ip1`	\((P\|Q)\), auxiliary centers	atom(\(P\)), atom(\(Q\))

The aux-center scatter is the one genuinely new structural piece the RI gradient has that the 4-center gradient does not: the auxiliary fitting functions are themselves atom-centered Gaussians, so they carry real, physical forces, and — critically — this auxiliary-basis response is not optional for qc-rs the way some codes treat it as a togglable approximation. Since qc-rs’s default JK-fitting basis (default_jk_aux) is atom-centered, omitting the aux-center terms would drop a real part of the physical force and fail a finite-difference check outright; only a hypothetical geometry-independent auxiliary basis could ever skip them.

A structural difference from the 4-center case is worth flagging directly: for 4-center, \(\sum_A\mathbf F_A=0\) is a check on internal consistency alone (permutation/transpose bookkeeping), since every contribution is already balanced by construction. For RI, the same sum-to-zero identity is a genuine detector of a missing auxiliary-response term — dropping the \((P|\mu\nu)\) aux-center derivative or the metric derivative breaks translational invariance measurably, not just subtly. And because RI is itself an approximation, the RI gradient should agree with the exact 4-center gradient only up to the RI fitting error (the same few-\(\times10^{-5}\)-scale offset seen in RI vs 4-center energies) — not to floating-point precision.

Verified example — water/cc-pVDZ, comparing the exact 4-center gradient against the RI gradient, and confirming translational invariance for both:

import qc, numpy as np
water = "O 0 0 0.117; H 0 0.757 -0.469; H 0 -0.757 -0.469"

m4c = qc.chk.new(atom=water, ao="cc-pvdz", unit="angstrom").scf(ref="r").run()
g4c = np.array(m4c.scf.gradient)
g4c.sum(axis=0)   # [-1.4e-29, 2.2e-15, -8.0e-15] -- translational invariance to ~1e-15

mri = qc.chk.new(atom=water, ao="cc-pvdz", unit="angstrom").ints(eri="ri-ram").scf(ref="r").run()
gri = np.array(mri.scf.gradient)
gri.sum(axis=0)   # [1.5e-27, 1.5e-15, -4.4e-16] -- also zero to numerical noise

np.abs(g4c - gri).max()   # 2.05e-05 -- matches the RI fitting-error scale seen in energies,
                          # not a floating-point-level agreement

Both gradients are translationally invariant to numerical noise (\(\sim10^{-15}\)), and the 4c/RI difference (\(\sim2\times10^{-5}\)) sits at exactly the RI approximation-error scale rather than the floating-point noise floor — precisely the two distinct diagnostic behaviors the design note predicts.

What analytic gradients do not need (and what does)#

This chapter’s whole derivation hinges on the energy being a stationary functional of the orbitals — that is what collapsed the orbital-response term to zero via Hellmann-Feynman. Two important consequences follow, both worth stating explicitly since they mark the boundary of this chapter’s scope:

DFT XC, ECP, and PCM contributions each add their own gradient term (\(\partial E_{\text{xc}}/\partial X\), etc.) but none of them reintroduce orbital response — they are all still variational contributions to the same stationary \(E\), so the same Hellmann-Feynman collapse applies to each. The XC term needs the grid weight’s own geometry derivative (the Becke partition weight derivative \(\partial W_A/\partial R_C\) derived in the DFT-grid chapter) plus the AO-follows- atom derivative \(\partial\chi/\partial X=-\nabla\chi\); the ECP term needs \(\partial\langle\chi|U_{\text{ECP}} |\chi\rangle/\partial X\) with no Pulay piece of its own (ECP is a genuine two-center operator matrix, each center differentiated directly, not an orthonormality constraint); PCM needs cavity/tessera-geometry and response derivatives on top of the same tessera-potential integral derivative.
Any non-variational quantity’s gradient is a different problem entirely. MP2 (and any other post-SCF correlation energy) is not stationary with respect to the reference orbitals — the reference determinant was optimized for the SCF energy, not the correlation energy — so its gradient genuinely needs the orbital response back, via the Z-vector/coupled-perturbed equations. That is explicitly out of scope for the SCF gradient this chapter derives, and is instead the subject of the next chapter.

Exercise 10

The verified example shows the 4c/RI gradient difference (\(\sim2\times10^{-5}\)) is roughly three orders of magnitude larger than either gradient’s own deviation from translational invariance (\(\sim10^{-15}\)). Explain why these two numbers measure completely different things, and why neither should be expected to shrink the other.
ROHF’s orbital_energies are eigenvalues of the effective Roothaan Fock operator, and using them directly in the RHF/UHF-style \(W\) formula gives a plausible-looking but wrong gradient. Explain in one sentence why this error would be easy to miss without a finite-difference check (i.e. why the wrong gradient still “looks reasonable”).
The 2-electron gradient has no Hellmann-Feynman term, but the 1-electron nuclear-attraction gradient does. Both operators formally act between two AO centers. What is the one structural difference between \(r_{12}^{-1}\) and \(1/|\mathbf r-\mathbf R_B|\) that explains this?

Solution to Exercise 10

The \(\sim10^{-15}\) translational-invariance residual measures pure floating-point rounding in a sum that is exactly zero by an underlying physical symmetry (a free molecule feels no net force) — it would not shrink with a better auxiliary basis or a different SCF method, only with higher-precision arithmetic. The \(\sim2\times10^{-5}\) 4c/RI difference measures a genuine physical/mathematical approximation — the same resolution-of-the-identity truncation that produces the RI energy fitting error — and it would shrink systematically if the auxiliary basis were enlarged. Neither number constrains the other: one is a floor set by machine precision, the other a ceiling set by basis quality.
The wrong gradient is still built from the same converged density \(D\), the same AO integrals, and a \(W\)-like object with the right shape/dimensions and roughly the right order of magnitude — it is not a crash, not a NaN, and not wildly off in scale, so nothing about running the calculation flags it as suspicious. Only comparing against a numerical (finite-difference) reference gradient — which requires no assumption about which Lagrangian multipliers are “correct” — would reveal that the direction/ magnitude is subtly off specifically in the open-shell block, since the error vanishes exactly when n_open=0 and is easy to overlook on a system with only a small open-shell contribution.
\(r_{12}^{-1}\) is the interaction between two electrons — it depends only on the electron coordinates \(\mathbf r_1,\mathbf r_2\), never on a nuclear position, so there is no nuclear coordinate for a Hellmann-Feynman-style operator-center derivative to act on; every geometry dependence enters only through where the AO basis functions themselves are centered. \(1/|\mathbf r-\mathbf R_B|\), by contrast, is the interaction between an electron and a nucleus — the operator itself contains the nuclear position \(\mathbf R_B\) explicitly, so differentiating the operator with respect to that same \(\mathbf R_B\) is a well-defined, nonzero, genuinely Hellmann-Feynman term in addition to the basis-center piece.

The energy-weighted density \(W\) and the Pulay-force structure derived here reappear, generalized, in the analytic Hessian chapter — the second nuclear derivative does need orbital response (the coupled-perturbed equations derived next), precisely because a second derivative of a stationary quantity is no longer itself stationary.