# Gradients & geometry optimization

So far every calculation used a **fixed** geometry — the coordinates you typed. But the geometry you typed is
rarely the *equilibrium* one. This chapter is about the **forces** on the nuclei (the energy gradient) and
using them to **relax a molecule to its equilibrium structure**. It is what turns "an energy at a guessed
geometry" into "the geometry nature actually adopts."

## Theory: the gradient is the force

Within the [Born–Oppenheimer](../10-foundations/many-electron-and-bo.md) picture the electronic energy is a
function of the nuclear positions, $E(\mathbf R)$ — the **potential energy surface (PES)**. The **force** on
atom $A$ is minus the gradient of that surface,

$$
\mathbf F_A = -\frac{\partial E}{\partial \mathbf R_A},
$$

and an **equilibrium structure** is a stationary point where every force vanishes, $\partial E/\partial
\mathbf R = 0$ (a minimum on the PES). qc-rs computes the gradient **analytically** — directly, from
derivative integrals, not by finite differences — which is fast and precise. A built-in self-test confirms
the forces obey **translational invariance**: $\sum_A \mathbf F_A \approx 0$ (a rigid shift of the whole
molecule cannot change the energy).

### Why the gradient needs no orbital response — but a Pulay term

Naively, $\partial E/\partial R_A$ looks like it should need $\partial\mathbf C/\partial R_A$ — how the
orbitals *themselves* respond as a nucleus moves. The **Hellmann–Feynman theorem** says it does not: because
the converged SCF energy is stationary with respect to orbital rotations ($\mathbf g=0$, the
[SCF convergence theory](../10-foundations/scf-convergence-theory.md) condition), every orbital-response term
cancels out of the *first* derivative, and only the operators' **explicit** $R$-dependence survives:

$$
\frac{\partial E}{\partial R_A}
= \sum_{\mu\nu} D_{\mu\nu}\,\frac{\partial h_{\mu\nu}}{\partial R_A}
+ \tfrac12\sum_{\mu\nu\lambda\sigma} D_{\mu\nu}D_{\lambda\sigma}\,\frac{\partial(\mu\nu|\lambda\sigma)}{\partial R_A}
\;-\; \sum_{\mu\nu} W_{\mu\nu}\,\frac{\partial S_{\mu\nu}}{\partial R_A}
\;+\; \frac{\partial V_{nn}}{\partial R_A},
$$

with the **energy-weighted (Lagrangian) density** $W_{\mu\nu}=2\sum_i \varepsilon_i\,C_{\mu i}C_{\nu i}$ (RHF;
UHF sums per spin, ROHF uses the projector form $W=\sum_\sigma P_\sigma F_\sigma P_\sigma$) standing in for
the ordinary density in the **last** one-electron term — the **Pulay force**, $-\operatorname{Tr}(WS^{[1]})$.

It exists for a purely geometric reason: qc-rs's Gaussian basis functions are **attached to the nuclei**, so
as a nucleus moves its AOs move with it — the overlap matrix itself depends on $R_A$ even though the orbital
*coefficients* need no further differentiation. Skip the Pulay term and the gradient is wrong by exactly this
basis-following contribution, one of the most common gradient bugs in the field. Two distinct kinds of
$R$-dependence appear in the operator terms above: an **AO-centered** part (the differentiated Gaussian
itself moves) and a **Hellmann–Feynman operator-centered** part (the nuclear-attraction operator
$-Z_B/|\mathbf r-\mathbf R_B|$ depends explicitly on *every* nucleus $B$, contributing a force at $B$ even
when the differentiated shell sits elsewhere). qc-rs's `qc-grad` crate assembles exactly these terms — kinetic,
nuclear attraction (both pieces), two-electron repulsion, the Pulay overlap term, $V_{nn}$, and ECP/PCM/D3-D4
when active — mirroring the [energy decomposition](scf.md) you already know from `energy_components`. The
full derivation, including the two-electron and DFT-XC gradient contractions, is in
[Analytic derivatives](../10-foundations/analytic-derivatives.md).

## The gradient at a geometry

After an SCF, the forces are one accessor away:

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

done = qc.chk.new(atom=water, ao="cc-pvdz", unit="angstrom").scf(ref="r").run()
grad = np.asarray(done.scf.gradient)          # shape (natom, 3), atomic units (Ha/bohr)

grad.shape                 # (3, 3)
np.abs(grad).max()         # 0.014317   the largest force component
grad.sum(axis=0)           # ~[0, 0, 0]  translational-invariance self-test
```

`done.scf.gradient` (equivalently `qc.grad(done)`) returns the `[natom, 3]` array of $\partial E/\partial
\mathbf R$ in **hartree/bohr**. A **non-zero** gradient (here max 0.0143) means the molecule is *not* at
equilibrium — the forces point downhill toward a better structure. The gradient covers **RHF/UHF/ROHF and
KS-DFT** (including hybrids), and automatically includes any **ECP, PCM, and DFT-D3/D4** contributions you
turned on.

## Geometry optimization

To *follow* those forces to the minimum, add an **`.opt()`** step. It drives the analytic gradient with the
**geomeTRIC** internal-coordinate optimizer, taking steps until the forces (and the displacement) fall below
threshold:

```python
opt = qc.chk.new(atom=water, ao="cc-pvdz", unit="angstrom").scf(ref="r", xc="b3lyp").opt().run()

opt.opt.converged     # True
opt.opt.energy        # -76.420628   energy at the optimized geometry
opt.opt.e_traj        # [-76.420349, ..., -76.420628]   energy per optimization step
opt.opt.gmax_traj     # max-gradient per step   (converges toward 0)
opt.opt.grms_traj     # rms-gradient per step
opt.coordinates()     # the optimized coordinates (bohr)
```

`.opt()` chains onto the SCF: each optimization step is a full SCF at a slightly moved geometry, and geomeTRIC
proposes the next move from the gradient. The final checkpoint holds the **optimized** structure and its
electronic state, ready for properties or a frequency analysis. The optimizer prints its own step-by-step
progress (energy, gradient, trust radius) as it runs.

```python
# read the relaxed geometry back
c = np.asarray(opt.coordinates()) * 0.52917721092   # bohr -> angstrom
# for this B3LYP/cc-pVDZ water: O–H ≈ 0.9687 Å, H–O–H ≈ 102.7°
```

:::{tip} `opt()` options
`opt(coordsys=..., maxiter=...)` exposes the common geomeTRIC controls — `coordsys` picks the coordinate
system (`"tric"` translation-rotation internal coordinates is the robust default) and `maxiter` caps the
number of steps. The full set is in the SCF/optimization reference.
:::

## Reading a converged optimization

The `*_traj` accessors are the story of the optimization — use them to check it behaved:

```python
import numpy as np
e = np.asarray(opt.opt.e_traj)
print("steps          :", len(e))                    # 4
print("energy lowered :", round(e[0] - e[-1], 6), "Ha")   # 0.000279
print("final max-grad :", float(np.asarray(opt.opt.gmax_traj)[-1]))  # ~1e-5, below threshold
```

A healthy optimization shows the energy **decreasing monotonically** and the max-gradient **shrinking toward
zero**. geomeTRIC's convergence criteria (Gaussian-style) require the energy change, the RMS/max gradient,
and the RMS/max displacement to *all* fall below their tolerances — which is why `converged=True` is a
stronger statement than "the energy stopped changing."

## Worked example: how much does relaxing help?

```python
import qc
water = "O 0 0 0.117; H 0 0.757 -0.469; H 0 -0.757 -0.469"   # a slightly-off input geometry

single = qc.chk.new(atom=water, ao="cc-pvdz", unit="angstrom").scf(ref="r", xc="b3lyp").run()
opt    = qc.chk.new(atom=water, ao="cc-pvdz", unit="angstrom").scf(ref="r", xc="b3lyp").opt().run()

print(f"single point : {single.scf.energy:.6f}")   # -76.420349
print(f"optimized    : {opt.opt.energy:.6f}")       # -76.420628
print(f"relaxation   : {(single.scf.energy - opt.opt.energy)*627.509:.3f} kcal/mol")  # 0.175
```

The optimized energy is lower (as it must be — relaxing can only descend the PES), here by ~0.18 kcal/mol
because the input was already close. For a poor starting geometry the gain is far larger, and the *structure*
(bond lengths, angles) is often what you actually wanted.

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

1. After an SCF at your input geometry, `np.abs(done.scf.gradient).max()` is `2e-6`. Is the molecule at
   equilibrium? What would `.opt()` do from here?
2. Why is `sum(gradient, axis=0) ≈ 0` a useful *self-consistency check* on a gradient implementation, no
   matter the molecule?
3. An optimization returns `converged=False` after `maxiter` steps, but `e_traj` is still dropping steadily.
   What is the likely cause and the fix?
:::

:::{solution} ex-grad
:class: dropdown

1. Essentially yes — a max force of `2e-6` Ha/bohr is below typical optimization thresholds, so it is at (a)
   stationary point. `.opt()` would take one or two tiny steps and report `converged=True` almost
   immediately.
2. A rigid translation of the whole molecule cannot change the energy, so the net force must be zero:
   $\sum_A \mathbf F_A = 0$ exactly. A computed gradient that violates this reveals a bug or an inconsistency
   (missing term, wrong sign) — it is a free, molecule-independent test.
3. The optimizer simply ran out of steps before reaching the (still-descending) minimum. Raise `maxiter`
   (`opt(maxiter=...)`), or restart the optimization from the last geometry, which the checkpoint already
   holds.
:::

With energies, correlation, forces, and optimized structures in hand, the remaining guide chapters add the
**environment** ([solvation & dispersion](solvation-dispersion.md)), **visualization**, and **logging**, then
the large [molecular-properties suite](properties/index.md).