Born–Oppenheimer Surfaces

A numerical Born–Oppenheimer surface is a map from nuclear geometry to a consistently identified electronic energy plus nuclear repulsion. Smooth-looking samples are not enough: state changes, incompletely converged electronic solutions, or basis-motion errors can create false barriers and incorrect forces. The molecular Born–Oppenheimer treatment explains when separating electronic and nuclear motion is justified.

Helpful background. Excited-state targeting supplies state tracking; basis-set convergence separates representation errors.

Energy surfaces and geometry conventions

At nuclear coordinates $\mathbf R$, solve the electronic problem and define

$$U_a(\mathbf R)=E_a^{\rm el}(\mathbf R)+E_{\rm nuc}(\mathbf R), \qquad E_{\rm nuc}=\sum_{A<B}\frac{Z_AZ_B}{R_{AB}}$$

in atomic units. Record whether an output energy already includes $E_{\rm nuc}$; adding it twice changes both minima and dissociation limits. Specify coordinate units, total charge, spin sector, and electronic approximation at every geometry.

Choose scan points to resolve curvature and avoided crossings, not only to look evenly spaced. Reuse the preceding solution as an initial guess, but also test independent starts where multiple self-consistent branches exist. Forward and backward scans can reveal hysteresis or branch switching. The lowest energy found at each point need not form one smooth metastable branch.

Derivatives in a moving basis

For a normalized generalized eigenpair $H(\mathbf R)c=E S(\mathbf R)c$ with $c^\dagger Sc=1$, differentiation gives

$$\partial_A E =c^\dagger(\partial_A H-E\,\partial_A S)c.$$

Terms containing $\partial_Ac$ cancel by stationarity. The overlap derivative is essential when basis functions move with nuclei. In many-electron orbital methods, analogous Pulay and response terms must be included according to the energy functional; the bare Hellmann–Feynman expectation is not universally the complete force in a finite basis.

Check analytic forces against centered finite differences of tightly converged energies over several step sizes. If electronic energy noise is $\delta E$, a first derivative can inherit noise of order $\delta E/h$ and a second derivative of order $\delta E/h^2$. Making the displacement arbitrarily small eventually worsens the estimate.

State identity and interpolation

At neighboring geometries, track overlaps, symmetry, spin, and appropriate density or orbital descriptors. Near an avoided crossing, the energy-sorted adiabatic state changes character smoothly. A diabatic model instead tries to retain character while allowing off-diagonal coupling; these are different representations.

For two diabatic energies $U_1(R),U_2(R)$ with real coupling $g$,

$$U_\pm(R)=\frac{U_1+U_2}{2} \pm\sqrt{\left(\frac{U_1-U_2}{2}\right)^2+g^2}.$$

At a crossing of the uncoupled curves, the adiabatic gap is $2|g|$. A scan spacing larger than the crossing width can miss this structure, causing an interpolant to invent a cusp or incorrect curvature. A vanishing gap also signals that nuclear dynamics may require multiple surfaces.

Fit an interpolation whose smoothness matches the intended use. Vibrational frequencies require second derivatives, so energy agreement alone is insufficient. Reserve points not used in the fit and compare both energies and forces there. Check dissociation limits and avoid extrapolating a local polynomial far outside its sampled region.

A harmonic surface check

For $U(R)=U_0+\tfrac12k(R-R_0)^2$, centered energy differences give the force exactly, and the centered second difference returns $k$ exactly in exact arithmetic. The vibrational frequency is $\omega=\sqrt{k/\mu}$. Adding a constant offset changes neither derivative nor frequency. These identities test geometry units, finite differences, and the separation of reference energy from curvature.

Vibrational numerics develops the nuclear eigenproblem once the surface is sufficiently accurate.

Exercises

If the energy noise is reduced by a factor of $100$, how much smaller can a finite-difference Hessian displacement be at unchanged noise amplification?

Solution

The Hessian noise scales as $\delta E/h^2$. Reducing $\delta E$ by $100$ allows $h$ to decrease by a factor of $10$ at the same noise level. Truncation error and correlated solver errors still need independent checks.

References

• Born, M., and K. Huang, Dynamical Theory of Crystal Lattices, Oxford University Press, 1954 — electronic-channel and nuclear-motion framework.
• Pulay, P., “Ab initio calculation of force constants and equilibrium geometries in polyatomic molecules. I. Theory,” Molecular Physics 17, 197–204 (1969), doi:10.1080/00268976900100941 — force derivatives in geometry-dependent bases.