The Molecular Electronic Hamiltonian and Born-Oppenheimer Approximation

Every commercial and academic electronic structure code — Gaussian, ORCA, Psi4, NWChem, VASP, Q-Chem — is built on the Born-Oppenheimer separation as its foundational assumption; the electronic energy $E_{el}(\mathbf{R})$ computed at a chosen level of theory (Hartree-Fock, DFT, coupled cluster) is precisely the quantity defined here, and geometry optimizations are literally a numerical search for minima of $U(\mathbf{R})$ on the resulting potential energy surface. Drug discovery pipelines depend on this at every stage: docking algorithms that estimate how a candidate molecule fits a protein active site implicitly assume each conformer sits at a local minimum of a Born-Oppenheimer surface, and free-energy perturbation calculations used to rank binding affinities are built from ensembles of nuclear configurations sampled on that same surface. In materials science, the entire enterprise of ab initio molecular dynamics (Car-Parrinello and Born-Oppenheimer molecular dynamics) propagates nuclei classically on an electronic energy surface recomputed at every timestep — directly implementing the two-tier structure derived here to simulate battery electrolytes, catalytic surfaces, and phase transitions in solids. Vibrational (infrared and Raman) spectroscopy, one of the most common analytical tools in any chemistry laboratory, is interpreted entirely through harmonic (or anharmonic) analysis of the curvature of $U(\mathbf{R})$ near its minimum — every reported vibrational frequency is a Born-Oppenheimer-surface property. The approximation's known failure modes are equally important industrially and scientifically. Photochemistry and photovoltaics — solar cell materials, photosynthetic pigments, photoswitches used in optogenetics — routinely involve conical intersections where two electronic potential energy surfaces cross, and there the clean separation into a single $U(\mathbf{R})$ breaks down; understanding non-adiabatic dynamics near these crossings is central to modern photochemical theory and to engineering better photostable dyes and light-harvesting complexes. Isotope effects exploited in mechanistic organic chemistry (kinetic isotope effects comparing C–H versus C–D bond cleavage rates) trace back to the mass-dependence built into the nuclear Schrödinger equation atop the (isotope-independent) Born-Oppenheimer electronic surface — the electronic surface $U(\mathbf{R})$ is identical for a molecule and its deuterated isotopologue, and all kinetic isotope effects arise purely from the different nuclear masses in the vibrational problem, a clean illustration of exactly where the Born-Oppenheimer boundary is drawn.

Ask the tutor: How would the derivation change if I wanted to include the diagonal Born-Oppenheimer correction directly in the electronic energy rather than treating it as a separate perturbative afterthought?

Consider the exact, non-relativistic molecular Hamiltonian for a collection of $N$ electrons and $M$ nuclei, written in atomic units ($\hbar = m_e = e = 4\pi\varepsilon_0 = 1$) so that energies come out in hartrees and lengths in bohrs. Every particle carries kinetic energy, and every pair of charged particles interacts through a Coulomb term. There are four fundamentally different kinds of interaction to track: nuclear kinetic energy, electronic kinetic energy, electron–electron repulsion, electron–nucleus attraction, and nucleus–nucleus repulsion. The full Hamiltonian is genuinely a many-body operator in $3N + 3M$ coordinates — for water alone that is 10 electrons and 3 nuclei, i.e., a function of 39 spatial coordinates. No exact solution exists in closed form for any molecule beyond one electron. The physical intuition that rescues the problem is a mass argument. A proton is 1836 times heavier than an electron; a carbon-12 nucleus is roughly 21,900 times heavier. Classically, a light particle bound to (or orbiting) a much heavier one moves far faster and adjusts to any change in the heavy particle's position almost instantly — think of a fly buzzing around a slowly rotating merry-go-round. Quantum mechanically the same intuition survives: the electronic wavefunction can "relax" to a new nuclear configuration on a timescale much shorter than the nuclei need to move appreciably. This licenses treating the nuclear positions as fixed, classical parameters when solving for the electronic wavefunction, and then afterward asking how the nuclei move on the electronic energy landscape those calculations produce. The picture that emerges is a two-tier structure. At the bottom tier, for every fixed nuclear geometry $\mathbf{R}$, one solves an electronic Schrödinger equation to get electronic wavefunctions $\psi_{el,n}(\mathbf{r};\mathbf{R})$ and electronic energies $E_{el,n}(\mathbf{R})$, where $n$ labels the electronic state (ground, first excited, etc.) and the semicolon signals that $\mathbf{R}$ enters only as a parameter, not a dynamical variable. At the top tier, the nuclei move under an effective potential $U_n(\mathbf{R}) = E_{el,n}(\mathbf{R}) + V_{NN}(\mathbf{R})$ — the potential energy surface — obeying their own Schrödinger equation for vibration, rotation, and translation. Every equilibrium bond length, every harmonic force constant, every reaction barrier height that a chemist quotes is a feature of this surface, generated entirely from the electronic-tier calculation.

The key mathematical object separating the two tiers is the total wavefunction ansatz. Rather than solving for one monstrous $\Psi(\mathbf{r},\mathbf{R})$ depending on all electron and nuclear coordinates simultaneously, Born-Oppenheimer factors it:

This factorization is an approximation, not an identity — it is exact only in the limit that the nuclear mass is infinite. The error introduced is controlled and, for most ground-state chemistry away from surface crossings, extremely small. The next figure makes explicit which Coulomb terms live in which tier of the Hamiltonian, since correctly sorting these terms is the crux of the whole derivation.

The exact non-relativistic molecular Hamiltonian, for $N$ electrons (indices $i,j$, mass $m_e=1$ in atomic units) and $M$ nuclei (indices $A,B$, masses $M_A$, charges $Z_A$), is

where:

• $\nabla_A^2, \nabla_i^2$ are Laplacians in the coordinates of nucleus $A$ and electron $i$ respectively (kinetic energy operators, $-\tfrac12\nabla^2$ in atomic units);
• $M_A$ is the mass of nucleus $A$ in units of the electron mass (so $M_A \gg 1$, typically $10^3$–$10^5$);
• $r_{iA}$ is the distance between electron $i$ and nucleus $A$; $r_{ij}$ is the distance between electrons $i$ and $j$; $R_{AB}$ is the distance between nuclei $A$ and $B$;
• $Z_A$ is the atomic number (nuclear charge) of nucleus $A$.

Reading: the five terms are, in order, nuclear kinetic energy $\hat{T}_N$, electronic kinetic energy $\hat{T}_e$, electron–nucleus attraction $\hat{V}_{eN}$ (negative, stabilizing), electron–electron repulsion $\hat{V}_{ee}$ (positive, destabilizing), and nucleus–nucleus repulsion $\hat{V}_{NN}$ (positive, destabilizing). This is simply Coulomb's law applied to every pair of charged particles plus a kinetic energy operator for each particle — nothing more exotic than that, but the coupling between all these terms through the shared coordinates is what makes it unsolvable exactly. Since $\hat{T}_N$ carries a factor $1/M_A$ that is smaller than the corresponding electronic factor by $10^3$–$10^5$, the Born-Oppenheimer approximation neglects it relative to the electronic terms when constructing the equation that determines the electronic wavefunction. What remains, for fixed nuclear positions $\mathbf{R}$, is the electronic Hamiltonian:

where:

• the nuclear positions $\mathbf{R} = (\mathbf{R}_1,\dots,\mathbf{R}_M)$ appear only as fixed numbers inside $r_{iA}$, not as dynamical variables to differentiate with respect to;
• $\hat{H}_{el}$ acts only on electronic coordinates $\mathbf{r}=(\mathbf{r}_1,\dots,\mathbf{r}_N)$.

Reading: this is a genuine, self-contained eigenvalue problem — electrons in the fixed external field of point charges, repelling one another. Solving

where:

• $n$ indexes electronic states (n = 0 for the ground state);
• $E_{el,n}(\mathbf{R})$ is the electronic energy eigenvalue, a function of the parameter $\mathbf{R}$ — changing $\mathbf{R}$ and resolving gives a new number;

gives the electronic wavefunctions and energies at one geometry. Repeating this calculation at many geometries $\mathbf{R}$ traces out the potential energy surface once the constant nuclear repulsion is added back:

Reading: $U_n(\mathbf{R})$ is exactly the potential energy surface that governs nuclear motion — the curve or surface with minima at stable geometries and saddle points at transition states, central to Potential Energy Surfaces. Nuclei then obey their own Schrödinger equation,

where:

• $v$ labels nuclear (vibrational/rotational) quantum states within electronic state $n$;
• $U_n(\mathbf{R})$ now plays the role of an ordinary potential energy function, exactly as in any textbook particle-in-a-potential problem.

Reading: this equation is precisely a "particle moving in a potential" problem, with the potential itself computed quantum mechanically from the electronic structure. It is the bridge from electronic structure theory to molecular spectroscopy and reaction dynamics.

Example 1 — Foundational: writing out $\hat H_{el}$ term by term for HeH$^+$. The helium hydride cation HeH$^+$ has two nuclei (He, $Z_A=2$; H, $Z_B=1$) and two electrons (it is isoelectronic with H$_2$, having lost one electron relative to neutral HeH). Fix the internuclear distance at $R = 1.463\ \text{bohr}$ (its approximate equilibrium bond length) and label electrons 1 and 2. The electronic Hamiltonian is

Every term is dictated directly by the general formula: two electronic kinetic energy operators (one per electron), four electron–nucleus attraction terms (two electrons $\times$ two nuclei, with $Z_A = 2$ for helium and $Z_B=1$ for hydrogen), and one electron–electron repulsion term $1/r_{12}$ — there is only one electron pair since $N=2$. Note there is no nuclear kinetic energy operator here at all: that is the entire content of the Born-Oppenheimer approximation as applied to this system. The nuclear repulsion, computed separately and added afterward as a constant shift at this fixed geometry, is

Solving $\hat H_{el}\psi_{el} = E_{el}\psi_{el}$ (numerically, since even two electrons in a molecule require approximate methods) at $R=1.463$ bohr gives an electronic energy of roughly $E_{el} \approx -4.36$ hartree (near the Hartree–Fock limit); adding $V_{NN}$ gives the total electronic-plus-nuclear-repulsion energy at that geometry, $U(R) \approx -2.99$ hartree, one point on the potential energy curve $U(R)$ for HeH$^+$. Repeating this whole calculation at a range of $R$ values and finding the minimum of $U(R)$ is exactly how the equilibrium bond length is predicted computationally — no experimental input required at all.

Example 2 — Applied: estimating the Born-Oppenheimer breakdown parameter for HF versus HI. The fractional error introduced by Born-Oppenheimer scales as roughly $\kappa = \sqrt{m_e/\mu}$, where $\mu$ is a reduced nuclear mass relevant to the vibration (from the order-of-magnitude derivation above). Compare hydrogen fluoride (H–F) to hydrogen iodide (H–I), both containing a light hydrogen but bonded to increasingly heavy partners. For HF, using atomic mass units and $m_e = 1/1836.15$ amu, the reduced mass of the two nuclei is

Converting to electron-mass units, $\mu_{HF} = 0.9572 \times 1836.15 = 1757.6\ m_e$, so

For HI, the reduced mass is

so in electron-mass units $\mu_{HI} = 1.00007\times 1836.15 = 1836.3\ m_e$, giving

The two ratios are nearly identical (both $\approx 0.023$–$0.024$) because in a hydride X–H the reduced mass is dominated by the light hydrogen atom regardless of how heavy the partner X is — the reduced mass saturates near $\mu \to m_H$ as $m_X \to \infty$. This shows explicitly why Born-Oppenheimer errors in ordinary hydride vibrational spectroscopy are essentially independent of the heavy partner and are governed almost entirely by the hydrogen mass, at the level of a few percent correction to already-small vibrational/electronic energy ratios — consistent with the fact that Born-Oppenheimer breakdown corrections to vibrational term values are typically parameterized empirically (Dunham-type Born-Oppenheimer breakdown corrections) at the sub-percent level of the vibrational spacing itself, i.e., of order $\kappa^2\sim 5\times10^{-4}$ relative to electronic energies.

Limiting case check: setting the nuclear mass to infinity, $M_A \to \infty$, sends every $1/M_A \to 0$, so $\hat T_N \to 0$ identically and $\kappa \to 0$: the Born-Oppenheimer approximation becomes exact in this limit, which is precisely the "clamped nuclei" idealization it is built on. This also confirms the formula behaves correctly at its defining limit rather than being an ad hoc fit.

Example 3 — Exam/competition level: estimating a diagonal Born-Oppenheimer correction (DBOC) using perturbation theory. A rigorous way to quantify how good Born-Oppenheimer is (beyond order-of-magnitude scaling) is to treat the discarded nuclear kinetic terms as a perturbation and evaluate the leading correction to the electronic energy. From Step 3 of the first derivation, the term retained within the Born-Oppenheimer electronic equation, if kept, is the expectation value of $\hat T_N$ evaluated in the electronic wavefunction itself — this is called the diagonal Born-Oppenheimer correction (DBOC):

Consider a highly simplified one-dimensional caricature of H$_2^+$: model the single electron's ground-state wavefunction dependence on the internuclear distance $R$ near equilibrium as a Gaussian in the reduced nuclear coordinate,

where $\alpha = 0.85\ \text{bohr}^{-2}$ describes how the electron localizes near one proton (a crude stand-in for the true molecular orbital). Differentiating twice with respect to the nuclear coordinate $R$ (playing the role of $\nabla_A$, with the electron coordinate $x$ held fixed during differentiation, exactly as required in the DBOC formula) gives

Taking the expectation value over the electronic Gaussian distribution, where $\langle (x-R/2)^2\rangle = 1/(4\alpha)$ for a Gaussian of this width (standard Gaussian second moment $1/(4\alpha)$),

Using $\alpha = 0.85\ \text{bohr}^{-2}$,

For a single proton of mass $M_p = 1836.15\ m_e$ (using the reduced nuclear mass for relative motion of the two protons, $M = M_p/2 = 918.08\ m_e$ for H$_2^+$'s symmetric stretch coordinate), the DBOC estimate is

Converting to more familiar units, $1.157\times10^{-4}\ \text{hartree} \times 27.211\ \text{eV/hartree} = 3.148\times10^{-3}\ \text{eV} \approx 3.15\ \text{meV}$, or equivalently about $25.4\ \text{cm}^{-1}$. This is tiny compared to a typical H–H bond dissociation energy of several eV, confirming numerically — with genuinely asymmetric, computed inputs rather than a hand-waved order of magnitude — that the diagonal correction neglected by Born-Oppenheimer is a few thousandths of an electron-volt, fully consistent with the $\sqrt{m_e/M}$ scaling estimate from the earlier order-of-magnitude derivation ($\kappa^2 E_{el}$ with $\kappa\approx 1/43$ for a proton and $E_{el}$ of order several eV gives corrections of order a few meV, matching this direct calculation's order of magnitude).