Minimum Energy Paths (MEPs)

• Minimum Energy Paths (MEPs) are continuous curves on potential energy surfaces where forces act tangentially, linking metastable minima with index-1 saddle points.
• They provide a precise framework to identify activation barriers in transition-state theory by ensuring the gradient aligns with the path at all points except at critical points.
• MEP computations integrate both continuous formulations, like the string method, and discrete approaches, such as NEB, enhanced by surrogate, invariant-set, and manifold-based acceleration techniques.

Minimum Energy Paths (MEPs) are curves on a potential-energy surface that connect metastable minima while remaining stationary with respect to all directions orthogonal to the curve. In the standard formulation, if $\gamma:[0,1]\to\mathbb{R}^d$ is a path with unit tangent $t(s)=\partial_s\gamma(s)/\|\partial_s\gamma(s)\|$, then an MEP satisfies

$$P_{\perp}(\gamma(s))\,\nabla V(\gamma(s)) = \big(I-t(s)t(s)^\top\big)\nabla V(\gamma(s)) = 0.$$

Equivalently, the gradient is everywhere parallel to the path tangent, except at critical points where $\nabla V=0$. The highest-energy point along an MEP is typically an index-1 saddle, and the associated barrier $\Delta E=\max_{s\in[0,1]}V(\gamma(s))-V(\gamma(0))$ enters transition-state theory through $k(T)=A\exp(-\Delta E/(k_BT))$ (Koten et al., 2018, Liu et al., 2022).

1. Geometric definition and energetic role

An MEP is defined for a smooth potential $V:\mathbb{R}^d\to\mathbb{R}$, usually with $V\in C^3$ or $C^4$, as a continuous curve connecting two local minima. Along the path, the force is tangent to the curve except at critical points. In the atomistic setting, the path typically connects minima through one or more index-1 saddles; in the formulation analyzed for the string method, MEPs may alternate minima and index-1 saddles, and kinks at minima are allowed because the parametrization need not be differentiable there (Koten et al., 2018).

The energetic significance of the MEP is twofold. First, it encodes the dominant transition mechanism between metastable states. Second, its maximum energy determines the activation barrier used in transition-state-theory rate estimates. In the discrete-operator formulation used for NEB analysis, the continuous barrier is

$$\delta E(\phi)=\sup_{\alpha\in[0,1]}E(\phi(\alpha))-E(\phi(0)),$$

with the discrete analogue

$t(s)=\partial_s\gamma(s)/\|\partial_s\gamma(s)\|$0

and the latter converges quadratically in the number of images under the assumptions of the discrete MEP theory (Liu et al., 2022).

A recurrent misconception is that vanishing perpendicular force on a discrete chain is automatically equivalent to an exact MEP. For finite systems, the quaternion-based NEB analysis emphasizes that a “zero gradient path” is not necessarily an MEP, and that verification at the highest-energy image should include confirming a first-order saddle and vanishing atomic forces at the saddle (Melander et al., 2015).

2. Continuous formulations and invariant-set viewpoints

One continuous viewpoint treats the MEP as an invariant set of the gradient flow

$t(s)=\partial_s\gamma(s)/\|\partial_s\gamma(s)\|$1

If $t(s)=\partial_s\gamma(s)/\|\partial_s\gamma(s)\|$2 denotes the associated flow map, then the gradient descent dynamics on curves (GDDC) evolves a curve pointwise by

$t(s)=\partial_s\gamma(s)/\|\partial_s\gamma(s)\|$3

Under this dynamics, an MEP is a fixed set because it is composed of stationary points and heteroclinic trajectories invariant under the flow (Koten et al., 2018).

The original continuous string method instead evolves a parametrized curve by a projected gradient flow,

$t(s)=\partial_s\gamma(s)/\|\partial_s\gamma(s)\|$4

together with an arclength constraint or explicit reparametrization step enforcing approximately uniform arclength density. This removes tangential drift, which would otherwise pull points toward minima, while preserving the orthogonal dynamics responsible for contraction toward the MEP (Koten et al., 2018).

A complementary variational formulation recasts the path as a family of constrained minimizers on hyperplanes normal to the curve. In the optimization-based string formulation, if $t(s)=\partial_s\gamma(s)/\|\partial_s\gamma(s)\|$5 is the hyperplane normal to the path at $t(s)=\partial_s\gamma(s)/\|\partial_s\gamma(s)\|$6, then the MEP satisfies

$t(s)=\partial_s\gamma(s)/\|\partial_s\gamma(s)\|$7

This converts MEP computation into a sequence of constrained minimization problems and permits projected CG, BFGS, L-BFGS, or FIRE updates in the normal subspace rather than pure steepest descent (Samanta et al., 2010).

A more operator-theoretic formulation introduces an admissible class of regular curves together with an operator $t(s)=\partial_s\gamma(s)/\|\partial_s\gamma(s)\|$8 enforcing both perpendicular-force stationarity and normalized arclength parametrization. In this framework, the MEP is the solution of $t(s)=\partial_s\gamma(s)/\|\partial_s\gamma(s)\|$9, and the weighted functional setting is designed to control endpoint and saddle singularities (Liu et al., 2022).

In gradient systems with small noise, MEPs also coincide, up to reparametrization, with minimum-action paths. The PINN-MEP formulation makes this connection through the Onsager–Machlup action and uses it to motivate continuous optimization of an entire path as a neural function rather than a finite image chain (Petersen et al., 23 Apr 2025).

3. Discrete chain-of-states methods

Discrete MEP algorithms approximate the path by images and differ mainly in how they enforce orthogonal relaxation and image distribution. In the standard NEB formulation, the force on image $$P_{\perp}(\gamma(s))\,\nabla V(\gamma(s)) = \big(I-t(s)t(s)^\top\big)\nabla V(\gamma(s)) = 0.$$0 is

$$P_{\perp}(\gamma(s))\,\nabla V(\gamma(s)) = \big(I-t(s)t(s)^\top\big)\nabla V(\gamma(s)) = 0.$$1

with

$$P_{\perp}(\gamma(s))\,\nabla V(\gamma(s)) = \big(I-t(s)t(s)^\top\big)\nabla V(\gamma(s)) = 0.$$2

and

$$P_{\perp}(\gamma(s))\,\nabla V(\gamma(s)) = \big(I-t(s)t(s)^\top\big)\nabla V(\gamma(s)) = 0.$$3

At stationarity, the discrete MEP is defined by $$P_{\perp}(\gamma(s))\,\nabla V(\gamma(s)) = \big(I-t(s)t(s)^\top\big)\nabla V(\gamma(s)) = 0.$$4 (Liu et al., 2022).

The string method uses the same perpendicular relaxation idea but controls image distribution by reparametrization rather than springs. In the “simplified and improved” variant analyzed rigorously, each image is first advanced by a full gradient-flow integrator and reparametrization is invoked only when the maximum segment length exceeds a threshold. This formulation replaces explicit tangent-projected force updates by full-gradient evolution plus redistribution, with the stability theory showing that tangential drift is neutralized by reparametrization (Koten et al., 2018).

The modified string method formulates each image update as constrained minimization on the hyperplane normal to the current path. The acceleration-space chain-of-states method instead treats the second derivatives of the path coordinates as the variational variables and updates

$$P_{\perp}(\gamma(s))\,\nabla V(\gamma(s)) = \big(I-t(s)t(s)^\top\big)\nabla V(\gamma(s)) = 0.$$5

then reconstructs the path by double integration, obtaining a built-in $$P_{\perp}(\gamma(s))\,\nabla V(\gamma(s)) = \big(I-t(s)t(s)^\top\big)\nabla V(\gamma(s)) = 0.$$6 damping of high-frequency path oscillations (Hernandez et al., 2015).

Method	Orthogonal relaxation	Path distribution
NEB	Projected true force	Springs along tangent
String	Projected flow or full-gradient step	Explicit reparametrization
Modified string	Minimization on normal hyperplanes	Periodic reparametrization
Acceleration method	Update normal accelerations	No springs or explicit redistribution required

For first-order finite-difference tangents, the discrete NEB theory proves that, for sufficiently small mesh spacing $$P_{\perp}(\gamma(s))\,\nabla V(\gamma(s)) = \big(I-t(s)t(s)^\top\big)\nabla V(\gamma(s)) = 0.$$7, there exists a stationary discrete MEP $$P_{\perp}(\gamma(s))\,\nabla V(\gamma(s)) = \big(I-t(s)t(s)^\top\big)\nabla V(\gamma(s)) = 0.$$8 satisfying

$$P_{\perp}(\gamma(s))\,\nabla V(\gamma(s)) = \big(I-t(s)t(s)^\top\big)\nabla V(\gamma(s)) = 0.$$9

and

$\nabla V=0$0

The same paper shows that the stationary discrete path of a first-order string method with equal-spacing redistribution coincides with the NEB discrete MEP (Liu et al., 2022).

The discrete string convergence theory is complementary rather than identical. For the simplified and improved string method, the asymptotic distance to the exact MEP is bounded by a function $\nabla V=0$1 that tends to zero as both $\nabla V=0$2 and $\nabla V=0$3 vanish, with contributions of order $\nabla V=0$4 from reparametrization, $\nabla V=0$5 from integrator truncation, and $\nabla V=0$6 from commutation of evolution and interpolation (Koten et al., 2018).

4. Stability, convergence, and perturbation theory

A major theoretical development is the recognition of the MEP as a stable invariant object. Under assumptions excluding degeneracies, higher-index saddles, and direct saddle-to-saddle heteroclinics, the MEP $\nabla V=0$7 is uniformly stable in Hausdorff distance under GDDC, and asymptotically stable with uniform convergence on a compact forward-invariant set containing a neighborhood of $\nabla V=0$8 (Koten et al., 2018).

The same analysis constructs a Lyapunov function $\nabla V=0$9 on curves satisfying

$\Delta E=\max_{s\in[0,1]}V(\gamma(s))-V(\gamma(0))$0

where $\Delta E=\max_{s\in[0,1]}V(\gamma(s))-V(\gamma(0))$1 denotes Hausdorff distance between compact curves. This device makes the long-time error analysis of the discrete string method possible because it prevents numerical perturbations from accumulating indefinitely (Koten et al., 2018).

In the operator formulation of MEP stability, the key result is that the linearization $\Delta E=\max_{s\in[0,1]}V(\gamma(s))-V(\gamma(0))$2 is an isomorphism under two structural assumptions: first, the path contains only two strong minima and a unique index-1 saddle; second, the endpoint tangents align with simple isolated minimal Hessian eigenvalues. Under these conditions, the MEP depends Lipschitz-continuously on perturbations of the energy landscape, with a quantitative estimate of the form

$\Delta E=\max_{s\in[0,1]}V(\gamma(s))-V(\gamma(0))$3

This provides an error-to-residual principle underlying subsequent convergence analyses of discrete MEP approximations (Liu et al., 2022).

These stability results also delimit the domain of validity of standard algorithms. Degenerate endpoint eigenvalues, higher-index saddles, direct saddle-to-saddle connections, or near-degenerate pathways can destroy local stability, reduce convergence rates, or invalidate uniqueness. The discrete NEB theory shows that dropping the endpoint spectral condition may degrade the path convergence rate even when barrier estimates remain comparatively robust (Liu et al., 2022).

5. MEPs on constrained and non-Euclidean configuration spaces

The MEP concept extends naturally beyond Euclidean atomic coordinates. In magnetic systems, the configuration space is the product manifold $\Delta E=\max_{s\in[0,1]}V(\gamma(s))-V(\gamma(0))$4 of unit spins, and the geodesic nudged elastic band (GNEB) replaces Euclidean differences by per-spin geodesic angles and tangent-space projections. The resulting method uses geodesic distances, projected torques, and energy-weighted springs; when combined with orthogonal spin optimization, LBFGS-OSO outperforms dynamics-based optimizers by up to a factor of $\Delta E=\max_{s\in[0,1]}V(\gamma(s))-V(\gamma(0))$5 in the reported tests on magnetization reversal, skyrmion collapse, and chiral bobber annihilation (Ivanov et al., 2020).

In variable-cell solid-state transitions, the configuration space includes both atomic coordinates and lattice vectors. The solid-state NEB (SSNEB) therefore operates in a $\Delta E=\max_{s\in[0,1]}V(\gamma(s))-V(\gamma(0))$6-dimensional space and combines atomic forces with stress contributions. In the hybrid ML+DFT framework, the SSNEB force is written as

$\Delta E=\max_{s\in[0,1]}V(\gamma(s))-V(\gamma(0))$7

and pressure-dependent paths are evaluated using the enthalpy $\Delta E=\max_{s\in[0,1]}V(\gamma(s))-V(\gamma(0))$8. For CsPbI$\Delta E=\max_{s\in[0,1]}V(\gamma(s))-V(\gamma(0))$9, GaN, and TiO$k(T)=A\exp(-\Delta E/(k_BT))$0, the hybrid framework converged to the same pathways predicted by first-principles calculations while achieving up to a $k(T)=A\exp(-\Delta E/(k_BT))$1-fold speedup (Zhang et al., 11 Jun 2026).

Liquid-crystal applications supply an alternative continuum setting. In chiral nematic cells, MEPs connect metastable helical states on a free-energy landscape defined by Frank elasticity, dielectric coupling, and Rapini–Papoular anchoring. GNEB calculations reveal multiple MEPs, including director-slippage and anchoring-breaking scenarios, and identify first-order saddle points with planar or strongly tilted localized structures. In this setting, the electric-field dependence of the energy barrier can exhibit hysteresis, and a first-order saddle-mediated transition may occur below the Freedericksz thresholds of the endpoint minima (Tenishchev et al., 2019, Tenishchev et al., 2020).

Finite, non-periodic molecular systems introduce a different geometric complication: overall translation and rotation. Quaternion-based rigid-body alignment removes these external degrees of freedom before NEB or dimer updates. In the reported benchmarks, this led to shorter paths, fewer required images, and significantly reduced iteration counts; the method was implemented in the Atomic Simulation Environment (ASE) for both NEB and DIMER workflows (Melander et al., 2015).

6. Initialization, acceleration, and surrogate-based MEP discovery

A persistent practical issue is that chain-of-states methods are highly sensitive to the initial path. Several initialization strategies therefore operate before any expensive PES optimization. S-IDPP constructs an initial path gradually on the image-dependent pair-potential surface, introducing images sequentially from the endpoints and controlling spacing with nonuniform springs; in reactions involving large rotations, this avoids the unnecessary bond breaking seen in linear-interpolation IDPP seeds (Schmerwitz et al., 2023). The velocity-field method in redundant internal coordinates generates a continuous seed path by integrating

$k(T)=A\exp(-\Delta E/(k_BT))$2

then projecting this velocity into Cartesian space through the internal-coordinate Jacobian. In the reported test cases, the resulting zero-temperature string calculations converged in markedly fewer geometry and SCF cycles than seeds from Cartesian interpolation, IDPP, or geodesic interpolation in two of three benchmarks, and its midpoint provided the closest initial TS guess in the HONO elimination example (Palenik, 2020). ASRBA offers another physically informed initialization route: it exploits the empirical observation that bond lengths along several DFT MEPs typically remain within about $k(T)=A\exp(-\Delta E/(k_BT))$3–$k(T)=A\exp(-\Delta E/(k_BT))$4 Å of endpoint values and uses adaptive semi-rigid-body forces to generate initial paths that reduced cNEB CPU time by approximately $k(T)=A\exp(-\Delta E/(k_BT))$5–$k(T)=A\exp(-\Delta E/(k_BT))$6 across the reported examples (Cai et al., 2022).

Acceleration during the MEP solve itself has followed several directions. Metric-based preconditioning replaces the Euclidean projector by a $k(T)=A\exp(-\Delta E/(k_BT))$7-weighted geometry and combines it with the adaptive timestepper ode12r. Across vacancy and dislocation examples modeled with both interatomic potentials and DFT, the reported gains reach approximately $k(T)=A\exp(-\Delta E/(k_BT))$8 and also enable convergence at tolerances that unpreconditioned schemes fail to reach (Makri et al., 2018). Gaussian-process surrogate surfaces reduce the number of expensive energy and force evaluations by repeatedly relaxing the path on a GP posterior mean and adding new true evaluations only where needed; for the heptamer island benchmark, the count was reduced to less than one fifth of regular NEB (Koistinen et al., 2017).

Recent neural and active-learning approaches replace the discrete path by trainable surrogates of either the PES or the path itself. PINN-MEP represents the entire path as a continuous neural function $k(T)=A\exp(-\Delta E/(k_BT))$9 and optimizes it directly against differentiable molecular force fields. For an explicitly hydrated BPTI system with about $V:\mathbb{R}^d\to\mathbb{R}$0 atoms, the reported cost was about $V:\mathbb{R}^d\to\mathbb{R}$1 force-field evaluations and about $V:\mathbb{R}^d\to\mathbb{R}$2 minutes on a single NVIDIA A6000 GPU, compared with approximately $V:\mathbb{R}^d\to\mathbb{R}$3 evaluations in the referenced Anton trajectory (Petersen et al., 23 Apr 2025). NN-BAX instead runs NEB on an actively updated neural surrogate of the PES; in Lennard–Jones and EAM benchmarks it achieved one to two order-of-magnitude reductions in energy and force evaluations, with reported speedup factors of $V:\mathbb{R}^d\to\mathbb{R}$4, $V:\mathbb{R}^d\to\mathbb{R}$5, $V:\mathbb{R}^d\to\mathbb{R}$6, and $V:\mathbb{R}^d\to\mathbb{R}$7 on four LJ pathways and barrier errors within $V:\mathbb{R}^d\to\mathbb{R}$8 of classical NEB (Kakhandiki et al., 17 Dec 2025).

When reliable forces are unavailable at the target theory level, force-free MEP identification becomes possible through surrogate-Hessian line searches. In the QMC setting, MEP images are optimized in subspaces orthogonal to the current path tangent using only stochastic energies, while a cheaper deterministic method supplies the local Hessian directions. The same framework refines first-order saddles by minimizing along positive-curvature directions and maximizing along the unique negative-curvature direction, and it was validated on ammonia inversion and an $V:\mathbb{R}^d\to\mathbb{R}$9 reaction against DFT and coupled-cluster benchmarks (Iyer et al., 2024). A distinct stochastic route is the AtomREM walker method, which evolves an ensemble under a transformed Smoluchowski dynamics and can discover minimum-energy escape paths from a single known minimum without coarse graining or a prescribed endpoint (Nagornov et al., 2019).

Taken together, these developments show that the modern theory of MEPs comprises more than a geometric condition on a PES. It also includes invariant-set stability, discrete approximation theory, manifold generalizations, initialization strategies, and surrogate-assisted algorithms that reduce the cost of traversing high-dimensional energy landscapes while preserving the defining requirement that the force be tangent to the path.