Potential Energy Surfaces from Electronic Energy

Every geometry optimization run inside a quantum chemistry package — whether computing the structure of a new drug candidate, a battery electrolyte molecule, or a catalytic transition-metal complex — is, at its core, a search for a stationary point on a potential energy surface built exactly as described here, just in far more than one dimension. Modern electronic structure methods (Hartree-Fock, density functional theory, coupled cluster) supply $E_{el}$ at a given geometry; software then walks along the multidimensional surface, computing gradients and curvatures (the multidimensional analogues of $U'(R)$ and $U''(R)$), until it finds a minimum corresponding to a stable molecular structure or a saddle point corresponding to a transition state. In drug design, comparing the PES of a small molecule bound inside a protein pocket versus free in solution is the computational essence of estimating binding affinity — the depth and shape of local minima on that combined surface determine which candidate molecules are predicted to bind tightly. In materials science, PES minima and their curvatures predict the elastic constants and phonon spectra of crystals, feeding directly into predictions of thermal conductivity and mechanical strength for battery materials and semiconductors. In combustion and atmospheric chemistry, locating saddle points (transition states) on multidimensional PESs, connecting reactant-side minima to product-side minima, is precisely how activation energies for reactions like ozone depletion or soot formation are computed before any experiment is run. In spectroscopy and astrochemistry, fitting observed rotational-vibrational transitions of interstellar molecules to Morse-like potentials remains the standard route to inferring bond lengths and force constants of species that can never be brought into a laboratory, from distant molecular clouds to planetary atmospheres. The entire enterprise of computational chemistry — force fields in molecular dynamics, ab initio geometry optimization, transition-state theory in reaction kinetics — is built on the premise introduced here: that a well-defined function $E_{el}(R)$ (or its multidimensional generalization) exists, that nuclei move on it, and that its minima and saddle points carry all the chemical information a structure or a reaction mechanism requires.

:::quiz[In the Morse potential $U(R) = D_e[1-e^{-a(R-R_e)}]^2$, what does the parameter $a$ control, given fixed $D_e$?]

• The dissociation energy directly
• The equilibrium bond length directly
• The curvature at the minimum (and hence the force constant), via $k = 2D_ea^2$ ✓
• The nuclear charge of the atoms

Ask the tutor: How does the variational principle guarantee that a computed electronic energy at a fixed geometry is never lower than the true ground-state energy, and why does that matter for finding the minimum of a PES?

Imagine holding two hydrogen nuclei apart with an imaginary clamp at a fixed distance $R$. The two electrons in H₂ do not know or care that the nuclei are clamped; they instantaneously settle into whatever electronic wavefunction minimizes the energy for that particular $R$. Solving the electronic Schrödinger equation

at this one clamped distance gives one number: $E_{el}(R)$. Now move the clamp — shorten $R$ to squeeze the nuclei together, or lengthen it to pull them apart — and solve the electronic problem again. Each value of $R$ produces its own electronic energy. Collecting all of these numbers as a function of $R$ produces a smooth curve. Add the classical nuclear-nuclear repulsion $Z_AZ_Be^2/(4\pi\epsilon_0 R)$ (which costs nothing extra to compute once you know $R$) and you obtain the total potential energy that governs nuclear motion,

This $U(R)$ is the potential energy surface for a diatomic molecule. It is not a surface in the everyday sense — with only one internal coordinate ($R$) a diatomic's "surface" degenerates to a curve — but for a polyatomic molecule with $N$ atoms there are $3N-6$ internal coordinates (bond lengths, angles, dihedrals), and $U$ becomes a genuine multidimensional hypersurface. Every foundational feature, however, is already visible in the one-dimensional diatomic case, which is why it is the natural starting point. The shape of $U(R)$ tells a complete story. At very large $R$, the two atoms feel almost no interaction, and $U(R)$ flattens out to the energy of two isolated, non-interacting atoms — call this the dissociation asymptote. As $R$ decreases from infinity, bonding orbitals begin to overlap constructively, electron density accumulates between the nuclei, and the energy drops: this is chemical bond formation. The energy keeps falling until, at some optimal separation $R_e$, the stabilization from shared electron density is exactly balanced against the cost of nuclear-nuclear and electron-electron repulsion pressing the atoms together. Push the nuclei closer than $R_e$ and repulsion dominates violently — the curve rises steeply, almost like a wall, because squeezing two positively charged nuclear cores together, and cramming electron clouds into overlapping regions in violation of the Pauli principle, both cost enormous energy. This asymmetry — a slow flattening on the long-range side and a steep wall on the short-range side — is the single most universal feature of every bonding curve in chemistry. The minimum of $U(R)$ therefore has direct physical meaning. Its position, $R_e$, is the equilibrium bond length — the number reported in every table of bond lengths. Its depth below the dissociation asymptote, called $D_e$, is the electronic dissociation energy: how much energy would be released bringing the atoms from infinite separation down to the bottom of the well (this differs slightly from the experimentally measured $D_0$, which includes the fact that even the lowest vibrational state has zero-point energy and sits above the bottom of the well, so $D_0 = D_e - \tfrac{1}{2}\hbar\omega$). The curvature of the well at its minimum — how sharply it turns upward — determines how stiff the bond is to stretching and compressing, and that stiffness is precisely what a vibrational spectrum measures.

Near the bottom of the well, any smooth curve looks approximately like a parabola — this is just a Taylor expansion argument, true for essentially any smooth function near its minimum. That is why the simplest possible model of a chemical bond, the harmonic oscillator, works as well as it does for small-amplitude vibrations: locally, $U(R)$ really is close to $\tfrac12 k(R-R_e)^2$. But a parabola never turns over — it rises forever — so it cannot describe dissociation. Real bonds break. A model that keeps the harmonic behavior near the minimum but flattens out correctly at large $R$, and rises steeply (though not infinitely) at short $R$, captures the full bonding curve; the Morse potential, developed in the derivations below, is the classic example.

The starting point is the clamped-nuclei electronic Schrödinger equation itself:

where:

• $\hat{H}_{el}$ is the electronic Hamiltonian — kinetic energy of electrons plus electron-nuclear attraction plus electron-electron repulsion — evaluated with the nuclei held fixed at separation $R$,
• $\psi_{el}(\mathbf{r};R)$ is the electronic wavefunction, a function of electron coordinates $\mathbf{r}$ that depends parametrically on $R$,
• $E_{el}(R)$ is the electronic energy eigenvalue, a single number for each chosen $R$.

Reading it: for every fixed value of the "dial" $R$, this equation is solved as an ordinary time-independent eigenvalue problem in the electron coordinates only; $R$ never appears as a dynamical variable in this step, only as a parameter. The total potential felt by the nuclei adds the classical repulsion between nuclear point charges:

where:

• $Z_A, Z_B$ are the nuclear charges of the two atoms,
• $e$ is the elementary charge, $\epsilon_0$ the vacuum permittivity,
• $U(R)$ is the potential energy surface itself — the function nuclei actually move on.

Reading it: $E_{el}(R)$ already contains all electronic effects (bonding, antibonding, exchange, correlation); the extra term is pure classical electrostatics between the nuclei, always positive and always repulsive, growing without bound as $R\to0$. At the equilibrium geometry, the surface is stationary:

where:

• $R_e$ is the equilibrium bond length,
• $k$ is the force constant, the curvature of the well at its bottom.

Reading it: this is exactly the calculus condition for a minimum — zero slope, positive curvature — applied to the electronic-energy-plus-repulsion curve; it is the rigorous definition of "bond length" and "bond stiffness" in terms of the PES. Near $R_e$, Taylor-expanding $U(R)$ and keeping only the quadratic term gives the harmonic approximation:

where:

• the linear term vanishes because $dU/dR=0$ at the minimum,
• $k$ is the same second derivative defined above.

Reading it: this is why every introductory treatment of molecular vibration starts with a spring; it is the mathematically inevitable local shape of any smooth potential minimum, not a special property of chemical bonds. Solving the nuclear Schrödinger equation with this harmonic potential (a textbook result, not re-derived here) gives vibrational energy levels

where:

• $v = 0,1,2,\dots$ is the vibrational quantum number,
• $\mu = m_Am_B/(m_A+m_B)$ is the reduced mass of the two nuclei,
• $\omega$ is the classical vibrational angular frequency of the harmonic oscillator equivalent to that bond.

Reading it: the force constant extracted from the curvature of the PES at $R_e$ directly predicts the observed infrared or Raman stretching frequency — this is the empirical bridge between the abstract PES and a number you can read off a spectrometer.

Example 1 — Foundational: extracting $R_e$, $D_e$, and $k$ conceptually from a tabulated (hypothetical) electronic energy scan. Suppose an electronic structure calculation on a diatomic molecule XY (not a real, symmetric textbook case — take masses $m_X = 19.0$ u, $m_Y = 35.0$ u, chosen asymmetrically) produces the following clamped-nuclei electronic energies plus nuclear repulsion, i.e., values of $U(R)$ in units of $10^{-19}$ J, at five bond lengths:

Step 1: locate the approximate minimum by inspection — the lowest value, $-6.62\times10^{-19}$ J, occurs near $R=1.40$ Å, so $R_e \approx 1.40$ Å (a proper treatment would fit a parabola through the three points bracketing the minimum, which we do next). Step 2: fit a parabola through the three central points $(1.30,-6.55)$, $(1.40,-6.62)$, $(1.50,-6.40)$, all in the stated units. Writing $U(R) = A + B(R-1.40) + C(R-1.40)^2$ with $R$ in Å: At $R=1.40$: $U = A = -6.62$. At $R=1.30$ ($\Delta R = -0.10$): $-6.55 = -6.62 -0.10B + 0.01C$. At $R=1.50$ ($\Delta R = +0.10$): $-6.40 = -6.62 +0.10B + 0.01C$. Subtracting the first from the second:

Adding the two equations:

The nonzero $B$ means the true minimum is not exactly at $R=1.40$ Å but shifted slightly; the vertex of the parabola is at

so the refined equilibrium bond length is

The corresponding minimum energy is

Taking the dissociation asymptote (energy of separated, non-interacting X and Y atoms) as, say, $U(\infty) = -5.20\times10^{-19}$ J (a value that would itself come from separate atomic calculations), the dissociation energy is

The force constant comes directly from the curvature coefficient $C$, remembering $U(R)=A+B\Delta R + C\Delta R^2$ means $d^2U/dR^2 = 2C$:

Converting to SI (1 Å $=10^{-10}$ m, so 1 Å$^{-2}=10^{20}$ m$^{-2}$):

This is a physically reasonable single-bond force constant (typical range roughly 200–1000 N/m). This example shows explicitly how three raw numbers on a computed energy scan — a set of $U(R)$ values, nothing more — become the three headline quantities every chemist quotes: bond length, bond energy, and bond stiffness.

Example 2 — Applied: predicting a vibrational frequency for HCl from a literature force constant, and checking the isotope effect. The force constant of HCl, obtained by exactly this parabola-fitting procedure applied to its measured (or computed) PES, is $k = 516.3$ N/m, and $R_e = 1.2746$ Å. Compute the fundamental vibrational wavenumber for $^1$H$^{35}$Cl, using atomic masses $m_H = 1.00783$ u and $m_{Cl} = 34.9689$ u ($1$ u $=1.66054\times10^{-27}$ kg). Step 1 — reduced mass:

Convert to kg:

Step 2 — angular frequency:

Step 3 — convert to wavenumber ($\tilde\nu = \omega/2\pi c$, $c=2.9979\times10^{10}$ cm/s):

This is very close to the experimental HCl fundamental (~2886 cm$^{-1}$ observed, with the small discrepancy arising because the true anharmonic $\omega_e$ is slightly larger than the observed $v=0\to1$ transition, which already includes anharmonic correction — exactly the effect built into the Morse treatment above). Step 4 — the isotope check (limiting-case verification). Repeat with deuterium, $m_D = 2.0141$ u, keeping the same force constant $k$ (a good approximation since $k$ depends on electronic structure, essentially mass-independent):

The ratio $\tilde\nu_{HCl}/\tilde\nu_{DCl} = 2990.6/2144.6 = 1.3945$, close to the expected $\sqrt{\mu_{DCl}/\mu_{HCl}} = \sqrt{1.9042/0.97962}=\sqrt{1.9439}=1.3942$ (the tiny residual difference is rounding). This is exactly the isotope-shift argument used routinely to confirm vibrational mode assignments, and it works only because $k$ — the curvature of the electronic PES — does not depend on nuclear mass at all, a direct consequence of the Born–Oppenheimer separation.

Example 3 — Exam/competition level: locating a PES stationary point from an anharmonic (cubic-corrected) energy expression, and testing whether it is a true minimum. A computational PES scan around the minimum of a triatomic bond stretch (treated as an effective one-dimensional problem in $R$, with all other coordinates relaxed) is fit to a quartic-truncated Taylor series about a trial point $R_0 = 1.55$ Å:

with, in units of $10^{-19}$ J and Å, $U_0 = -8.400$, $b_1 = -0.960$, $b_2 = 11.20$, $b_3 = -6.500$. Note $b_1 \neq 0$: $R_0$ is not the true minimum, only a nearby trial geometry, so the task is to find the actual stationary point and classify it. Step 1 — differentiate and set to zero. Let $x = R-R_0$:

Setting this to zero:

Apply the quadratic formula:

Two roots:

Step 2 — classify each stationary point using the second derivative,

At $x_2 = 0.04459$:

a genuine local minimum. At $x_1 = 1.1041$:

a local maximum (physically, this cubic-truncated fit is only valid locally; far from $R_0$ it is not a trustworthy extrapolation of the true PES, but within the fitted region it correctly signals a turning point). Step 3 — report the corrected equilibrium bond length and force constant:

This example is deliberately built with a nonzero linear term and an explicit cubic term so that neither the location of $R_e$ nor the classification of the stationary point can be guessed by symmetry — both require solving the actual stationary-point equations, exactly as a real anharmonic PES fit demands, and exactly the kind of calculation that motivates the systematic geometry-optimization algorithms built on the variational principle.

The mathematical shape used almost universally to model diatomic bonding curves is named for Philip M. Morse, an American physicist whose 1929 paper, "Diatomic Molecules According to the Wave Mechanics," proposed the exponential functional form still taught today. Morse was working at Princeton in the earliest years after Schrödinger's and Heisenberg's papers had established quantum mechanics, at a moment when physicists and chemists were racing to see whether the new theory could actually reproduce the spectroscopically observed vibrational structure of real molecules — band spacings, anharmonicity, dissociation limits — rather than merely offer qualitative bonding pictures. His key insight was mathematical economy: rather than solving the nuclear Schrödinger equation numerically for an arbitrary tabulated PES (a formidable numerical task in 1929, done by hand), he asked what functional form was both physically sensible (correct dissociation limit, a genuine minimum, finite curvature) and exactly solvable. The exponential-squared form he chose satisfies both demands simultaneously, an elegant piece of applied mathematics that remains in every physical chemistry course nearly a century later. Morse's later career moved well beyond diatomic spectroscopy — he became a central figure in the development of theoretical acoustics, underwater sound propagation during and after the Second World War, and operations research, eventually helping found the field that would become known as management science. But it is the single-parameter exponential curve from his 1929 paper, immortalizing the shape of the chemical bond in a closed analytical formula, that chemists invoke every time they fit a spectroscopic band system or sketch a bonding curve on a blackboard.