# 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](../20-guide/gradients-geomopt.md) 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](linear-response-cphf.md)).
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](linear-response-cphf.md) 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](post-hf-correlation.md)), 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](density-fitting-ri.md)). 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](density-fitting-ri.md), 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](solvation-theory.md) 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:

```python
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](dft-xc-quadrature.md)) 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](linear-response-cphf.md).

:::{exercise}
:label: ex-analytic-grad-theory

1. 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.
2. 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").
3. 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} ex-analytic-grad-theory
:class: dropdown

1. 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.
2. 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.
3. $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](analytic-hessian-thermo.md) — the second nuclear derivative *does* need orbital
response (the coupled-perturbed equations [derived next](linear-response-cphf.md)), precisely because a
second derivative of a stationary quantity is no longer itself stationary.