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

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.