Potential Energy Surfaces - Physical Chemistry 2

SECTION 1

Historical Context & Motivation

The concept of a potential energy surface (PES) arose from a fundamental question in early quantum chemistry: given that molecules are quantum mechanical systems, how can we predict the energetics and dynamics of chemical reactions in a systematic way? Classical kinetic theory, embodied by Arrhenius's empirical rate law, provided parameters like the activation energy but offered no molecular-level picture of why a barrier exists or how nuclear motion drives bond-breaking and bond-forming events. The PES concept bridged quantum mechanics and chemical kinetics by providing a multidimensional energy landscape on which nuclear dynamics unfold.

1927 • Born–Oppenheimer Approximation

Max Born and J. Robert Oppenheimer publish their landmark paper separating electronic and nuclear motion in molecules. By treating the nuclei as slow, classical-like particles moving on an electronic energy landscape, they lay the theoretical foundation for the potential energy surface.

1929 • London–Eyring–Polanyi Surface

Fritz London extends Heitler–London valence bond theory to the H + H₂ exchange reaction. Henry Eyring and Michael Polanyi then construct the first semi-empirical PES for a chemical reaction, revealing a saddle-point transition state connecting reactant and product valleys.

1935 • Transition State Theory

Eyring, and independently Evans and Polanyi, formulate transition state theory (TST), which directly uses the PES concept to derive rate constants by assuming quasi-equilibrium between reactants and the activated complex at the saddle point.

1968 • First Ab Initio PES for H₃

Bowen Liu publishes the first fully ab initio PES for the H + H₂ system, demonstrating that accurate quantum chemical calculations can reproduce the essential features of reactive surfaces without empirical parameters.

1980s–Present • Computational Revolution

Advances in density functional theory, coupled-cluster methods, and machine-learned potentials enable the construction of global PESs for polyatomic systems, powering accurate molecular dynamics simulations and the prediction of reaction mechanisms.

The central question the PES framework addresses is deceptively simple: given a collection of atoms, how does the total electronic energy depend on the positions of all nuclei, and how does this energy landscape dictate the pathways, rates, and products of chemical transformations? The remainder of this lesson develops the theoretical machinery to answer that question rigorously.

SECTION 2

Core Principles & Definitions

A potential energy surface is a hypersurface that expresses the electronic energy of a molecular system as a continuous function of its nuclear coordinates. For a system of N atoms, the PES is a function of 3N − 6 internal degrees of freedom (3N − 5 for linear molecules), after removing translations and rotations. Understanding the topology of this surface—its minima, saddle points, and gradient paths—is essential for predicting molecular structure, spectroscopy, and reaction dynamics.

Born–Oppenheimer Approximation

Electrons adjust instantaneously to nuclear positions because m_e ≪ m_N. The electronic Schrödinger equation is solved at fixed nuclear geometries, yielding an energy eigenvalue E(R) that defines the PES.

Stationary Points

Points where the gradient ∇E = 0. Minima (all positive curvatures) correspond to stable species; first-order saddle points (one negative curvature) correspond to transition states connecting minima.

Minimum Energy Path (MEP)

The steepest-descent path in mass-weighted coordinates from the transition state down to reactant and product minima. The MEP defines the intrinsic reaction coordinate (IRC) and serves as the backbone of reaction dynamics.

Dimensionality

For a nonlinear N-atom system, the PES lives in 3N − 6 dimensions. A triatomic (e.g., H₃) has 3 internal coordinates; a pentatomic has 9. Full PES construction becomes exponentially expensive with system size.

Adiabatic vs. Diabatic Surfaces

Adiabatic surfaces are eigenvalues of the electronic Hamiltonian and may exhibit avoided crossings. Diabatic surfaces are smooth, state-specific representations connected by coupling terms, useful near conical intersections.

✦

KEY TAKEAWAY

Think of a potential energy surface as a topographic map for molecules. Just as a hiker traversing mountain terrain seeks the lowest pass between two valleys, reacting molecules follow the minimum energy path across the PES from the reactant valley, over the saddle-point pass (the transition state), and down into the product valley. The height of the pass determines the activation energy, and the steepness of the terrain controls how fast the system moves.

SECTION 3

Visualizing the Potential Energy Surface

For the simplest reactive system—three hydrogen atoms undergoing the collinear H + H₂ → H₂ + H exchange—the PES can be represented as a contour map in two dimensions: the distance r_AB between atoms A and B, and the distance r_BC between atoms B and C. The resulting diagram reveals the fundamental topological features that govern all reactive PESs: reactant and product valleys separated by a saddle point.

Contour map of the collinear A + BC potential energy surface. The reactant valley (lower right) and product valley (upper left) are separated by a saddle point (transition state). The dashed minimum energy path (MEP) traces the lowest-energy route from reactants to products. Concentric contours represent equipotential lines, with the corners indicating high repulsive energy (both bonds compressed simultaneously).

Several features of this diagram deserve careful attention. In the reactant valley, atom A is far from the BC molecule, so r_AB is large while r_BC is near its equilibrium value. In the product valley, the situation is reversed: the AB bond has formed and C is departing. The saddle point represents the transition state geometry where both bonds are partially formed. The corners of the contour map, where both internuclear distances are small, correspond to strongly repulsive configurations—attempting to compress all three atoms together simultaneously requires enormous energy. The minimum energy path smoothly connects the two valleys through the saddle point and defines the intrinsic reaction coordinate used in transition state theory.

SECTION 4

Mathematical Framework

The mathematical backbone of the PES is the electronic Schrödinger equation solved within the Born–Oppenheimer approximation. At each fixed set of nuclear coordinates R, the electronic Hamiltonian Ĥ_el is diagonalized to obtain adiabatic electronic energies. These energies, parametrically dependent on R, define the potential energy surface on which nuclei move.

ELECTRONIC SCHRÖDINGER EQUATION

Ĥ_el(r; R) ψ_n(r; R) = E_n(R) ψ_n(r; R)

Here r denotes all electronic coordinates, R denotes all nuclear coordinates (treated as parameters), ψ_n is the nth adiabatic electronic wavefunction, and E_n(R) is the corresponding electronic energy eigenvalue—the PES for electronic state n.

Near a stationary point, the PES can be expanded as a Taylor series in the nuclear displacement coordinates. This expansion connects the PES to observable molecular properties: the gradient gives forces on the nuclei, and the second-derivative matrix (the Hessian) determines harmonic vibrational frequencies.

TAYLOR EXPANSION OF THE PES

E(R) ≈ E(R₀) + gᵀΔR + ½ ΔRᵀ H ΔR + ⋯

Here R₀ is the reference geometry, g = ∂E/∂R is the gradient vector (force), H = ∂²E/∂R_i∂R_j is the Hessian matrix, and ΔR = R − R₀. At a stationary point, g = 0, so the Hessian governs the local curvature.

HESSIAN EIGENVALUE CRITERION

Minimum: all λᵢ > 0 | Transition State (1st-order saddle): exactly one λᵢ < 0

The eigenvalues λ_i of the mass-weighted Hessian matrix determine the character of a stationary point. Positive eigenvalues correspond to restoring forces (bound vibrations); a single negative eigenvalue indicates the direction of imaginary frequency along the reaction coordinate.

HARMONIC FREQUENCY FROM HESSIAN

ωᵢ = √(λᵢ / μᵢ)

Each normal mode frequency ω_i is obtained from the corresponding Hessian eigenvalue λ_i divided by the reduced mass μ_i. At a transition state, the single negative eigenvalue yields an imaginary frequency, confirming the saddle-point nature.

SECTION 5

Key Features & Classification of PES Topology

The topology of a potential energy surface determines the chemical behavior of the system. Different types of stationary points and connecting pathways encode information about equilibrium structures, isomerization, and reaction mechanisms. Beyond stationary points, features like conical intersections (points where two adiabatic PESs become degenerate) play critical roles in photochemistry by providing ultrafast non-radiative decay channels between electronic states.

One-dimensional energy profile along the intrinsic reaction coordinate (IRC). The reactant minimum and product minimum are separated by the transition state. The forward and reverse activation energies (E_a,fwd and E_a,rev) are marked, along with the reaction energy ΔE_rxn.

Feature	Gradient	Hessian Eigenvalues	Physical Meaning
Global Minimum	∇E = 0	All λᵢ > 0	Most stable molecular geometry
Local Minimum	∇E = 0	All λᵢ > 0	Metastable conformer or isomer
1st-Order Saddle Point	∇E = 0	One λᵢ < 0	Transition state for a single reaction step
2nd-Order Saddle Point	∇E = 0	Two λᵢ < 0	Hilltop connecting multiple saddle points; not a true TS
Conical Intersection	Undefined (degenerate)	N/A (two PESs touch)	Non-adiabatic decay funnel between electronic states

Classification of critical points on a potential energy surface

ℹ

Dimensionality Note

For a six-atom molecule like ethane (C₂H₆), the PES has 3(8) − 6 = 18 internal coordinates. Visualizing or computing such a high-dimensional surface directly is intractable; in practice, one explores local regions around stationary points or projects the surface onto one or two key coordinates (e.g., a dihedral angle for conformational analysis).

SECTION 6

Worked Example: Characterizing a Stationary Point

Consider a triatomic molecule with three internal coordinates (two bond lengths and one bond angle). An ab initio geometry optimization locates a stationary point, and the mass-weighted Hessian matrix is computed. The task is to determine whether the stationary point is a minimum or a transition state, and, if it is a transition state, to identify the reaction coordinate.

Classifying a Stationary Point from Hessian Eigenvalues

1
Step 1 — State the Problem

A geometry optimization on a bent triatomic ABC at the CCSD(T)/cc-pVTZ level converges to a stationary point with ∇E = 0. The mass-weighted Hessian is diagonalized, yielding eigenvalues λ₁ = 5.82 mdyn/Å, λ₂ = 1.04 mdyn/Å, and λ₃ = −0.37 mdyn/Å. Classify this stationary point.
2
Step 2 — Count Negative Eigenvalues

The Hessian has three eigenvalues. Two are positive (λ₁ = 5.82, λ₂ = 1.04) and one is negative (λ₃ = −0.37). The number of negative eigenvalues (the index of the stationary point) is 1.

Index = 1 → This is a first-order saddle point (transition state).
3
Step 3 — Identify the Reaction Coordinate

The eigenvector associated with the single negative eigenvalue λ₃ defines the direction of the imaginary vibrational mode. This eigenvector points along the reaction coordinate—the direction in which the system descends toward reactants on one side and products on the other. In practice, one would visualize this eigenvector to determine which bonds are stretching or compressing at the transition state.

The reaction coordinate is defined by the eigenvector of λ₃.
4
Step 4 — Compute the Imaginary Frequency

Since λ₃ < 0, the harmonic frequency for this mode is imaginary: ω₃ = √(|λ₃|/μ₃) × i, where μ₃ is the reduced mass for that normal mode. Force constants from a Hessian in mdyn/Å convert to SI units using 1 mdyn/Å = 100 N/m, so |λ₃| = 0.37 mdyn/Å = 37 N/m. Suppose μ₃ = 1.67 × 10⁻²⁷ kg (approximately one atomic mass unit). Then ω₃ = √(37 / 1.67 × 10⁻²⁷) = √(2.22 × 10²⁸) ≈ 1.49 × 10¹⁴ s⁻¹. Converting to wavenumbers: ν̃ = ω₃/(2πc) ≈ 1.49 × 10¹⁴/(2π × 3.0 × 10¹⁰) ≈ 790 cm⁻¹ × i.

The transition state exhibits one imaginary frequency (≈ 790i cm⁻¹), a magnitude typical of real transition states and confirming its saddle-point character.
5
Step 5 — Physical Interpretation

The stationary point is a transition state connecting two minima on the PES. Following the intrinsic reaction coordinate (IRC) in both directions from this saddle point—by integrating the steepest-descent path in mass-weighted coordinates—leads to the reactant and product geometries. The two positive eigenvalues correspond to real vibrational modes orthogonal to the reaction path, representing oscillations that the molecule undergoes as it traverses the transition state.

An IRC calculation from this TS will reveal the full reaction path and identify the reactant and product minima.

SECTION 7

Methods for Constructing PESs: Strengths & Limitations

Constructing an accurate PES requires solving the electronic Schrödinger equation at many nuclear geometries, then fitting or interpolating between these points. The choice of electronic structure method and fitting strategy involves fundamental trade-offs between accuracy, computational cost, and the dimensionality of the system.

Method	Strengths	Limitations
Ab initio (CCSD(T), MRCI)	High accuracy (sub-kcal/mol for small systems); systematically improvable; no empirical parameters	Steep computational scaling (N⁷ for CCSD(T)); limited to small molecules for global PES construction
DFT (B3LYP, ωB97X-D)	Favorable N³ scaling; applicable to large molecules; reasonable barriers with hybrid functionals	Functional-dependent errors; poor for dispersion, charge transfer, and multireference systems; no systematic improvement path
Semi-empirical (PM7, GFN2-xTB)	Very fast; suitable for conformational searches and pre-screening; thousands of geometries per minute	Parameterized for specific element sets; quantitatively unreliable for barriers and reaction energies
Analytic Fitting (Polynomial, Spline)	Continuous, smooth surfaces; fast evaluation for dynamics; can enforce correct asymptotic behavior	Requires careful choice of functional form; poor extrapolation; labor-intensive for high dimensionality
Machine-Learned Potentials (NNP, GAP)	Near ab initio accuracy at force-field cost; flexible; excellent for high-dimensional systems; active learning reduces training data needs	Black-box nature; may fail outside training distribution; requires careful validation; energy conservation issues in MD

Comparison of methods for PES construction

⚗️

PRACTICAL GUIDANCE

In modern computational chemistry, a common workflow involves using DFT to explore the PES topology (locating minima and transition states), then refining energetics at the CCSD(T) or multireference level at key stationary points. For dynamics simulations requiring millions of energy evaluations, machine-learned potentials trained on high-level ab initio data are increasingly the method of choice, offering the accuracy of quantum chemistry at near-classical force-field speeds.

SECTION 8

Connection to Reaction Dynamics & Advanced Theory

The PES serves as the foundation for virtually all theories of reaction dynamics, from classical trajectory methods to full quantum scattering calculations. The level of dynamical theory applied on top of the PES determines the sophistication and accuracy of the predicted observables—rate constants, product-state distributions, angular distributions, and branching ratios.

Dynamical Method	PES Requirement	Key Capability
Transition State Theory (TST)	Only stationary-point energies and frequencies (local PES information)	Thermal rate constants via partition functions; no dynamical detail
Classical Trajectory (CT)	Global, smooth, analytic PES with continuous gradients	State-resolved cross sections, angular distributions; misses quantum effects (tunneling, ZPE)
Quantum Wave Packet	Global PES on a numerical grid	Full quantum dynamics including tunneling, resonances, interference; computationally demanding
Ab Initio Molecular Dynamics (AIMD)	On-the-fly electronic structure (no pre-fitted PES)	Explore unknown PES regions; expensive but avoids fitting errors; classical nuclei
Non-Adiabatic Dynamics (Surface Hopping)	Multiple coupled PESs with non-adiabatic coupling vectors	Photochemical reactions, conical intersections, excited-state dynamics

Dynamical methods and their PES requirements

Looking forward, the interplay between machine learning and quantum chemistry is transforming PES construction. Neural network potentials and Gaussian approximation potentials, trained on datasets of thousands of ab initio energies and forces, can now represent global PESs for systems with dozens of atoms at near-CCSD(T) accuracy. Coupled with advanced sampling methods (replica exchange, metadynamics), these representations are enabling the study of complex reaction networks, enzyme catalysis, and materials chemistry at unprecedented levels of detail. The PES concept, introduced nearly a century ago, remains the conceptual cornerstone upon which these modern developments rest.

SECTION 9

Practice Problems

Question 1: CONCEPTUAL

Explain the physical basis for the Born–Oppenheimer approximation. Under what circumstances does it break down, and what features of the PES signal this breakdown?

Show answer

The Born–Oppenheimer approximation rests on the large mass disparity between electrons and nuclei (m_p/m_e ≈ 1836), which means electrons respond nearly instantaneously to nuclear displacements. This allows the electronic problem to be solved at fixed nuclear geometries, generating a PES. The approximation breaks down when two adiabatic PESs approach each other closely (avoided crossing) or become degenerate (conical intersection). At such points, the non-adiabatic coupling terms—which couple nuclear motion to changes in the electronic wavefunction—become large and cannot be neglected. The PES features signaling breakdown include conical intersections and regions of small energy gaps between electronic states.

Question 2: BASIC CALCULATION

A nonlinear molecule contains 5 atoms. How many internal degrees of freedom define its PES? If you wished to compute the PES on a grid with 10 points along each coordinate, how many single-point energy calculations would be required?

Show answer

For a nonlinear N-atom molecule, the PES dimensionality is 3N − 6. With N = 5, we have 3(5) − 6 = 9 internal degrees of freedom. A grid with 10 points in each of 9 dimensions requires 10⁹ = 1,000,000,000 (one billion) single-point calculations. This exponential scaling with dimensionality—the so-called 'curse of dimensionality'—is precisely why full grid-based PES construction is impractical for systems beyond triatomics or tetraatomics, motivating the development of interpolation schemes and machine-learned potentials.

Question 3: INTERMEDIATE

You perform a frequency calculation on a stationary point of a tetraatomic molecule (4 atoms, nonlinear) and obtain the following mass-weighted Hessian eigenvalues in mdyn/Å: λ₁ = 8.21, λ₂ = 3.45, λ₃ = 1.12, λ₄ = 0.89, λ₅ = −0.56, λ₆ = 0.74. (a) How many internal degrees of freedom should there be? (b) Classify this stationary point. (c) What does the eigenvector of λ₅ represent physically?

Show answer

(a) For a nonlinear tetraatomic, 3(4) − 6 = 6 internal degrees of freedom, consistent with the six eigenvalues provided. (b) Among the six eigenvalues, exactly one (λ₅ = −0.56 mdyn/Å) is negative, so the stationary point is a first-order saddle point, i.e., a transition state. (c) The eigenvector corresponding to λ₅ defines the direction of the reaction coordinate in the 6-dimensional internal coordinate space. This is the single mode with imaginary frequency; following it in both directions via an IRC calculation traces the minimum energy path from the transition state to the connected reactant and product minima. The five positive eigenvalues correspond to real vibrational frequencies orthogonal to the reaction path.

Question 4: APPLIED

In a computational study of the SN2 reaction Cl⁻ + CH₃Br → ClCH₃ + Br⁻, you calculate the energies (relative to separated reactants) of the following stationary points on the PES: reactant ion-dipole complex = −42 kJ/mol, transition state = +12 kJ/mol, product ion-dipole complex = −50 kJ/mol, separated products = −8 kJ/mol. (a) Sketch a qualitative 1D energy profile along the reaction coordinate. (b) What is the activation energy for the forward reaction, measured from the reactant complex? (c) Is the overall reaction exothermic or endothermic?

Show answer

(a) The energy profile is a double-well potential: the energy drops from separated reactants (0 kJ/mol) to the reactant ion-dipole complex (−42 kJ/mol), rises to the central barrier at the SN2 transition state (+12 kJ/mol), descends to the product ion-dipole complex (−50 kJ/mol), then rises slightly to separated products (−8 kJ/mol). (b) The forward activation energy from the reactant complex to the transition state is E_a = (+12) − (−42) = 54 kJ/mol. Note that the activation energy relative to separated reactants is only +12 kJ/mol—the ion-dipole complexation provides significant stabilization that must be overcome. (c) The overall reaction is exothermic by 8 kJ/mol (products at −8 kJ/mol relative to reactants at 0 kJ/mol), consistent with the greater leaving-group ability of Br⁻ compared to Cl⁻.

Question 5: CRITICAL THINKING

Consider two electronic states of a diatomic molecule whose adiabatic PES curves undergo an avoided crossing at internuclear distance R = 2.5 Å. (a) Explain why the adiabatic curves do not actually cross for a homonuclear diatomic where both states have the same symmetry. (b) Under what kinematic conditions would you expect the Born–Oppenheimer approximation to fail most severely at this avoided crossing? (c) How would you set up a computational study to properly treat the dynamics in this region? Discuss the choice between adiabatic and diabatic representations.

Show answer

(a) The non-crossing rule (von Neumann–Wigner theorem) states that two adiabatic states of the same spatial and spin symmetry cannot cross in a one-dimensional parameter space (here, R). Instead, the electronic Hamiltonian matrix elements mix the two states, producing an avoided crossing where the energy gap is minimized but non-zero. This gap is proportional to the off-diagonal coupling matrix element. (b) The BO approximation fails most severely when the nuclear velocity is high relative to the energy gap at the avoided crossing—specifically, when the Massey parameter ξ = ΔE × Δx/(ℏv) is of order unity, where ΔE is the energy gap, Δx is the spatial width of the coupling region, and v is the nuclear velocity. Fast-moving, light nuclei (e.g., hydrogen) traversing a narrow, small-gap avoided crossing will exhibit the greatest non-adiabatic transition probability. (c) For dynamics near the avoided crossing, the diabatic representation is preferred because the diabatic surfaces and coupling terms vary smoothly, avoiding the singular non-adiabatic coupling vectors that arise in the adiabatic representation. One would compute diabatic states and their coupling using methods like block diagonalization or the fourfold-way, then propagate nuclear dynamics using trajectory surface hopping (Tully's fewest-switches algorithm) or full quantum scattering on the coupled diabatic surfaces.

SUMMARY

Summary

A potential energy surface is the multidimensional function E(R) that maps nuclear geometry to electronic energy, made possible by the Born–Oppenheimer approximation, which separates electronic and nuclear motion. For an N-atom nonlinear molecule, the PES has 3N − 6 internal degrees of freedom. The topology of the PES—its minima (stable species), saddle points (transition states), and connecting minimum energy paths—encodes the complete energetic and mechanistic information about a chemical reaction.

Stationary points are classified by the Hessian eigenvalue signature: all positive for minima, exactly one negative for transition states. PESs are constructed using electronic structure methods ranging from ab initio quantum chemistry (highest accuracy) to machine-learned potentials (greatest efficiency), and serve as the foundation for dynamical theories including transition state theory, classical trajectory simulations, and quantum wave packet dynamics. When the Born–Oppenheimer approximation fails—near conical intersections or avoided crossings—non-adiabatic dynamics methods on coupled PESs become essential.