Gradients & geometry optimization

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 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 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 you already know from energy_components. The full derivation, including the two-electron and DFT-XC gradient contractions, is in Analytic derivatives.

The gradient at a geometry#

After an SCF, the forces are one accessor away:

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:

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.

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

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?#

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 6

  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?

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