Docs are under construction. Content may be incomplete or change.

Nudged Elastic Band Method

Page type: Concept
On this page

The nudged elastic band (NEB) method approximates a minimum-energy path between two known states on a potential-energy surface. It is useful when the endpoint structures are known but the intervening mechanism and barrier are not.

NEB does not discover arbitrary products, prove that a mechanism is kinetically dominant, or by itself establish a transition state. It optimizes a discretized path under a chosen energy-and-force model. The resulting path is only as meaningful as the endpoints, atom correspondence, model chemistry, and numerical convergence.

The path on a potential-energy surface

For nuclear coordinates R\mathbf R, a calculator defines a potential energy V(R)V(\mathbf R) and force

F(R)=V(R).\mathbf F(\mathbf R)=-\nabla V(\mathbf R).

A path R(s)\mathbf R(s) connects the initial state RA\mathbf R_A and final state RB\mathbf R_B, with s[0,1]s\in[0,1]. Along a minimum-energy path, the physical force perpendicular to the path vanishes:

V ⁣(R(s))=0.\nabla V\!\left(\mathbf R(s)\right)_{\perp}=\mathbf 0.

The energy generally varies along the path. A local maximum on a converged minimum-energy path is a candidate first-order saddle, but that statement is local to the path and calculator. A transition-state assignment additionally requires refinement and examination of the Hessian or an equivalent vibrational analysis.

Discretizing the path

NEB represents the continuous path by NN images,

R0,R1,,RN1,\mathbf R_0,\mathbf R_1,\ldots,\mathbf R_{N-1},

where R0=RA\mathbf R_0=\mathbf R_A and RN1=RB\mathbf R_{N-1}=\mathbf R_B. Tako holds the endpoint images fixed and optimizes only the N2N-2 interior images.

Each image is a full copy of the atomic system. Consequently, atom ii in every image must refer to the same physical atom. Element-by-index validation catches mismatched species, but it cannot detect the exchange of two chemically equivalent atoms.

The images are joined by artificial springs. Those springs distribute images along the band; they are not physical bonds and do not contribute to the reported electronic energy profile.

@startuml skinparam backgroundColor white skinparam defaultFontName Helvetica skinparam titleFontColor #24292f skinparam titleBackgroundColor White skinparam titleBorderColor White title Discrete NEB image and nudged-force decomposition rectangle "Image i-1" as A rectangle "Image i" as B rectangle "Image i+1" as C rectangle "Physical force\nperpendicular to tangent" as P A -right-> B : spring force parallel to tangent B -right-> C : spring force parallel to tangent B -up-> P note bottom of B NEB removes the physical force parallel to the band and the spring force perpendicular to the band. end note @enduml

Nudging the forces

Let τ^i\hat{\boldsymbol\tau}_i be the unit tangent to the band at interior image ii. The physical force is decomposed into components parallel and perpendicular to the band:

Fi=(Fiτ^i)τ^i,Fi=FiFi.\mathbf F_i^{\parallel} =\left(\mathbf F_i\cdot\hat{\boldsymbol\tau}_i\right)\hat{\boldsymbol\tau}_i, \qquad \mathbf F_i^{\perp} =\mathbf F_i-\mathbf F_i^{\parallel}.

The NEB force retains the perpendicular physical force and the parallel spring force:

FiNEB=Fi+Fi,spring.\mathbf F_i^{\mathrm{NEB}} =\mathbf F_i^{\perp}+\mathbf F_{i,\mathrm{spring}}^{\parallel}.

For the common equal-spring form,

Fi,spring=k(Ri+1RiRiRi1)τ^i.\mathbf F_{i,\mathrm{spring}}^{\parallel} =k\left( \lVert\mathbf R_{i+1}-\mathbf R_i\rVert -\lVert\mathbf R_i-\mathbf R_{i-1}\rVert \right)\hat{\boldsymbol\tau}_i.

This projection is the defining “nudge,” and it exists to remove two characteristic failures of an unprojected elastic band. The physical force parallel to the path drives images downhill toward the endpoint minima — sliding. The perpendicular component of the spring force pulls a curved chain toward the chord between endpoints — corner cutting. NEB removes exactly those two offending components: the physical force determines the band in transverse directions, while the spring force only distributes images along it.

The spring constant has units of energy per squared length; Tako accepts it in eV Å2^{-2} and converts it for the atomic-unit backend. A very weak spring allows poor spacing and bunching in low-gradient regions. A very strong spring stiffens the artificial longitudinal modes, restricts useful redistribution around a climbing image, and can make the combined optimization ill-conditioned. Because the spring is not part of the physical potential, a converged, adequately resolved calculation should not depend on it: strong sensitivity of the apparent mechanism or barrier to the spring constant is a diagnostic of poor resolution, initialization, or convergence.

Tangent construction

The tangent determines which part of the physical force is removed. A simple central difference,

τ^iRi+1Ri1,\hat{\boldsymbol\tau}_i\propto\mathbf R_{i+1}-\mathbf R_{i-1},

can become unstable near extrema and corners, letting the tangent switch arbitrarily around the barrier and kink the band.

The improved-tangent construction uses the neighboring energies instead. Define the forward and backward displacement vectors

τi+=Ri+1Ri,τi=RiRi1.\boldsymbol\tau_i^+=\mathbf R_{i+1}-\mathbf R_i, \qquad \boldsymbol\tau_i^-=\mathbf R_i-\mathbf R_{i-1}.

When the energies increase monotonically through image ii,

Ei+1>Ei>Ei1,τi=τi+,E_{i+1}>E_i>E_{i-1}, \qquad \boldsymbol\tau_i=\boldsymbol\tau_i^+,

and when they decrease monotonically,

Ei+1<Ei<Ei1,τi=τi.E_{i+1}<E_i<E_{i-1}, \qquad \boldsymbol\tau_i=\boldsymbol\tau_i^-.

At a local extremum, the segments are combined with energy-difference weights. Let

ΔEi+=Ei+1Ei,ΔEi=Ei1Ei,\Delta E_i^+=|E_{i+1}-E_i|, \qquad \Delta E_i^-=|E_{i-1}-E_i|,

and let ΔEimax\Delta E_i^{\max} and ΔEimin\Delta E_i^{\min} be the larger and smaller magnitudes. Then

τi={τi+ΔEimax+τiΔEimin,Ei+1>Ei1,τi+ΔEimin+τiΔEimax,Ei+1<Ei1.\boldsymbol\tau_i= \begin{cases} \boldsymbol\tau_i^+\Delta E_i^{\max} +\boldsymbol\tau_i^-\Delta E_i^{\min}, & E_{i+1}>E_{i-1},\\[4pt] \boldsymbol\tau_i^+\Delta E_i^{\min} +\boldsymbol\tau_i^-\Delta E_i^{\max}, & E_{i+1}<E_{i-1}. \end{cases}

The selected vector is normalized before projection. The weighting biases the tangent toward the higher-energy neighbor, which matters most near a maximum: the tangent should represent the ascent toward the saddle rather than cut diagonally between two descending branches. Nearly equal neighbor energies remain numerically delicate because small energy noise can change the case selection; dense images and stable calculator forces help more there than the formal case distinction.

Tangent quality affects both components of the NEB force. An erroneous tangent removes the wrong component of the physical force and applies the spacing force in the wrong direction. The resulting band can kink, cut a corner, or drift into an adjacent valley while still reducing its reported residual.

Band formulations

Tako exposes several formulations because no tangent and spacing scheme is universally optimal:

MethodCharacterAppropriate use
Improved tangentEnergy-weighted local tangent with conventional springsDefault for most reaction and diffusion paths
ASENEBASE-style tangent and spring projectionCompatibility and comparison with established ASE workflows
Full springRetains a fuller spring-force treatmentDiagnostic alternative when spacing is problematic
SplineRepresents the band through a smooth interpolantSmooth paths where spline geometry is useful
StringReparametrizes images along a path coordinateAlternative spacing dynamics without conventional NEB springs

The full-spring formulation retains a more complete elastic-band spring treatment than the conventional projected form; it can help when band spacing behaves poorly, but it also changes artificial transverse influences and is not a more physical force. ASENEB’s value is compatibility when reproducing an established ASE workflow. Spline and string variants address spacing through a continuous representation: the spline method constructs tangents from a smooth band representation, and the string variant reparametrizes images along the fitted path rather than relying on a longitudinal spring. For periodic systems, the backend internally uses consecutive minimum-image unwrapping for spline coordinates and string reparametrization so a boundary crossing can remain a continuous path — but that internal treatment does not repair incorrect endpoint atom mapping or make independently wrapped endpoints safe.

No spacing method prevents convergence to a wrong mechanism. A smooth spline through chemically meaningless images remains chemically meaningless; a uniformly parametrized string can follow the wrong side of a bifurcation. Two numerically converged formulations can follow different mechanisms when the initialization lies near competing pathways, so changing the formulation is never a substitute for examining the band.

Initializing the band

Linear interpolation

Cartesian interpolation is componentwise and deterministic. For atom aa with endpoint positions raA\mathbf r_a^A and raB\mathbf r_a^B,

ra,i=(1λi)raA+λiraB,λi=iN1.\mathbf r_{a,i} =(1-\lambda_i)\mathbf r_a^A +\lambda_i\mathbf r_a^B, \qquad \lambda_i=\frac{i}{N-1}.

The construction assumes that the endpoint coordinate difference represents the intended displacement. It does not align independently oriented molecules, permute identical atoms, infer bond changes, or select minimum-image displacements. Rigid rotation represented as independent straight atomic chords can shrink or collide a molecule; fragment exchange can place atoms on top of one another. It is transparent and appropriate for small displacements, and unreliable during bond exchange, molecular rotation, diffusion around an obstacle, or concerted rearrangement.

In a periodic cell, equivalent coordinates differing by a lattice vector can generate an artificial traversal across the entire cell. If equivalent endpoint atoms are stored in different periodic images, wrap and remap the endpoints consistently before constructing the band — for a diffusion hop this may require placing the final atom in a translated periodic image rather than applying a generic wrap independently to both endpoints. These steps define the path problem; they are not optional preprocessing.

IDPP interpolation

The image-dependent pair potential (IDPP) method improves the seed band by interpolating pair distances rather than Cartesian components. If dab(Ri)d_{ab}(\mathbf R_i) is the distance between atoms aa and bb in image ii, and

dabtarget(λi)=(1λi)dab(RA)+λidab(RB),d_{ab}^{\mathrm{target}}(\lambda_i) =(1-\lambda_i)d_{ab}(\mathbf R_A) +\lambda_i d_{ab}(\mathbf R_B),

then a representative IDPP objective is

SIDPP(Ri)=a<bwab,i[dab(Ri)dabtarget(λi)]2,S_{\mathrm{IDPP}}(\mathbf R_i) =\sum_{a<b} w_{ab,i} \left[d_{ab}(\mathbf R_i) -d_{ab}^{\mathrm{target}}(\lambda_i)\right]^2,

with weights that emphasize short distances, commonly proportional to an inverse fourth power of distance. The interior images are relaxed against this fictitious objective to the requested IDPP threshold or step limit; the physical calculator is not evaluated.

Pair-distance interpolation often avoids direct overlaps because close contacts receive large weight, and it represents rotations and collective rearrangements more naturally than independent Cartesian chords. Nevertheless, the objective knows neither electronic structure nor reaction chemistry: it can preserve an unwanted contact, break a covalent framework, or choose a smooth path through an implausible topology. Atom correspondence remains fixed, and chemically equivalent permutations remain invisible to validation. IDPP is an initialization procedure, not a mechanism.

Interpolation-only execution

Running the interpolation branch without optimization is valuable precisely because it exposes these assumptions without invoking the calculator. In Tako that branch returns the constructed frames, reports neb.converged: true and zero steps to indicate successful seed generation, and omits the barrier, saddle, residual, and vibration result. The zero energy and force-like fields in its shared frame schema are placeholders, not calculated values. Scientific interpretation begins only after a physical band optimization.

Climbing-image NEB

Ordinary NEB distributes images along the minimum-energy path but does not require an image to sit at the path maximum. Climbing-image NEB selects the highest-energy interior image, removes its spring force, and reverses the physical force parallel to the path:

FiCI=Fi2(Fiτ^i)τ^i.\mathbf F_i^{\mathrm{CI}} =\mathbf F_i-2\left(\mathbf F_i\cdot\hat{\boldsymbol\tau}_i\right) \hat{\boldsymbol\tau}_i.

The perpendicular component continues to pull the image toward the path, while the inverted parallel component drives it uphill toward the local maximum.

Tako stages climbing by default, because the highest image of an unrelaxed seed is often an interpolation artifact. With a positive reactionPreclimbSteps budget the runtime first runs ordinary NEB with climbing disabled, then starts a fresh climbing optimization with the full reactionNebMaxSteps budget. This two-phase behavior applies to every exposed optimizer. The reported neb.steps and optimizer trace describe the final run, not the sum of the phases; with the registered defaults, up to 30 ordinary steps can precede up to 150 climbing steps. Ordinary NEB uses the final budget directly, and a zero pre-climb setting makes CI-NEB climb in its only optimizer run.

The identity of the climbing image may change while the band relaxes, and the fresh second phase permits it. Switching is not automatically an error — a neighboring image can become the better representation of the maximum. Persistent alternation, a flat barrier region, or a climbing image at the edge of the band indicates that the saddle region is not yet robustly resolved: possible flat ridges, inadequate resolution, calculator noise, or two nearby mechanisms. A stable image identity is supporting evidence, not an independent convergence criterion.

The calculation in Tako

The following diagram describes the implementation boundary. It shows which validation and numerical stages belong to the operation; it is not a replacement for the equations above.

@startuml skinparam backgroundColor white skinparam defaultFontName Helvetica skinparam titleFontColor #24292f skinparam titleBackgroundColor White skinparam titleBorderColor White title Tako reaction-path execution branches start :Read initial and final structures; :Validate atom count and element identity by index; if (Endpoints compatible?) then (yes) :Construct linear or IDPP seed band; if (Interpolation only?) then (yes) :Return seed images with unevaluated zero-valued energies; else (no) :Relax band with projected NEB forces; if (Climbing image enabled?) then (yes) :Invert the parallel physical force\non the highest-energy image; endif :Evaluate residual and step budget; :Return band, barrier, saddle candidate,\noptimizer traces, and optional vibration result; endif :Publish path frames and structured artifacts; else (no) :Reject before numerical execution; endif stop @enduml

Tako performs the following sequence:

  1. Validate endpoint atom count, element-by-index identity, image count, and normalized settings. Endpoint cell equality is a user invariant and is not checked by the engine (Known Limitations).
  2. Copy endpoint structures and propagate constraints from the initial structure.
  3. Construct the interior images using linear or IDPP initialization.
  4. Create an independent calculator instance for each active image.
  5. Evaluate energies and forces.
  6. Construct tangents and projected NEB forces.
  7. Update the interior images with the selected optimizer.
  8. For staged CI-NEB with any exposed optimizer, run the ordinary band for up to reactionPreclimbSteps, then start a fresh climbing optimizer run with its own reactionNebMaxSteps budget.
  9. Stop when the residual satisfies the threshold or the step budget is exhausted.
  10. Store the full band, energy profile, normalized settings, residual, convergence state, and saddle candidate after successful return. Cancelled runs do not produce a restartable partial band.

All images use the same calculator definition. Separate calculator instances prevent mutable electronic state from leaking between images. The cost of a band step is dominated by the force evaluations on its interior images; increasing image count raises cost approximately in proportion to the number of active images.

For periodic structures, Tako currently treats the cell as fixed across the band. Variable-cell NEB is not implemented. Constraints propagate from the initial endpoint, so a fully fixed atom remains fixed at its initial coordinates throughout the band. Fixing an atom that occupies a different position in the final endpoint silently conflicts with the intended path.

Optimization and convergence

Tako reports convergence against the maximum per-image residual from the backend. For improved-tangent, ASENEB, and full-spring bands, each image residual is the largest atomic projected-force vector norm. Spline and string variants instead use the preconditioner’s largest absolute transformed-force component. Only the final scalar, in eV Å1^{-1}, is exposed as neb.residual_ev_per_angstrom. Numerical thresholds with the same displayed unit therefore refer to closely related but not identical aggregations across formulations; comparing residual numbers between methods without acknowledging the metric can be misleading.

The projected NEB force is generally nonconservative — it is not the gradient of a single scalar potential for the band — so optimizer behavior must be judged from the residual and the path, not from an assumption that a conventional energy decreases monotonically. BFGS is the default because accumulated curvature information converges efficiently on smooth bands; poor initial images or abrupt tangent changes can make its curvature model temporarily unreliable, so the maximum image displacement remains an important stabilizer. FIRE and MDMin use damped dynamical ideas and can dissipate large initial forces robustly. The ODE optimizer adapts its integration using a relative tolerance and the c1c_1, c2c_2, and minimum-step controls; the static optimizer applies a direct scaled update. These are numerical tools for a diagnosed convergence problem: changing optimizer does not change the calculator, and materially different converged paths from different optimizers should be investigated rather than averaged.

Tako also exposes dynamic relaxation. When enabled, the backend can suppress updates for nonclimbing images whose residuals already satisfy a distance-scaled threshold relative to the highest image. This can reduce unnecessary movement of locally settled images, but it changes which parts of the band are active during a step; inspect a dynamic-relaxation result for spacing and continuity, particularly when the highest-image identity changes.

Convergence means that the discretized band is stationary under the selected NEB formulation and tolerance. It does not establish that:

  • the endpoints are the relevant chemical states;
  • the band follows the globally lowest path;
  • the calculator describes the changing bonding correctly;
  • the band contains enough images to resolve the saddle region;
  • the climbing image is a first-order saddle in the full coordinate space.

Tightening the force threshold reduces numerical residual but cannot correct model error, endpoint error, or a wrong mechanism. Likewise, a longer step budget only helps when the optimizer is progressing toward an appropriate path.

Selecting the number of images

Image count controls spatial resolution, not accuracy in isolation.

  • A simple, smooth, single-event path often begins with seven total images.
  • A path with a shoulder, intermediate, concerted rearrangement, or long diffusion distance may require 11–15 images.
  • Adding images is useful when adjacent structures near the maximum differ too much to identify the saddle region.
  • Adding images is not useful when the endpoint correspondence is wrong or the seed follows the wrong topology.

Convergence with respect to image count should compare the mechanism, saddle geometry, and barrier — not the barrier alone. If increasing the count reveals a new intermediate, the lower-resolution calculation did not represent the same path.

Constraints

Constraints can remove irrelevant motion, such as relaxation of deep substrate layers during surface diffusion. They can also create false convergence.

The NEB condition concerns force perpendicular to the path in the available coordinate space. If constraints remove those perpendicular degrees of freedom, an interpolated band may have no direction in which to relax. It can then satisfy a numerical criterion without approaching the unconstrained minimum-energy path.

Use constraints only when they represent a declared physical approximation. Apply the same constraint meaning to both endpoints, and confirm that no constrained atom is expected to move during the reaction.

Diagnosing a band

ObservationNumerical interpretationScientific question
Images bunch in one regionSpacing control, tangent, or resolution is inadequateIs there a shoulder or intermediate requiring more images?
Energies oscillate between neighboring imagesBand is kinked, optimizer is unstable, or step size is too largeDoes the seed attempt incompatible structural changes?
Interior images lie below both endpointsBand found a lower-energy basinIs the supposed transition actually mediated by an intermediate?
Highest image changes repeatedlyBarrier region is flat or poorly resolvedAre multiple mechanisms or saddle regions competing?
Residual stalls above thresholdOptimizer, tangent, or calculator is not producing progressIs the model smooth and valid for every image geometry?
Converged band is chemically implausibleNumerical stationarity occurred on an unintended pathAre endpoints, mapping, model, and spin state correct?
Climbing image has multiple imaginary modesCandidate is not a first-order saddleDoes it require further refinement or a different path?

Always inspect the image trajectory. An energy plot alone hides atomic crossing, fragment rotation, periodic wrapping errors, and unintended changes in coordination.

Energies and barriers

For image energies E0,,EN1E_0,\ldots,E_{N-1}, Tako defines the forward electronic barrier as

ΔEf=maxi(Ei)E0\Delta E_{\mathrm f}^{\ddagger}=\max_i(E_i)-E_0

and the reaction energy as

ΔErxn=EN1E0.\Delta E_{\mathrm{rxn}}=E_{N-1}-E_0.

When the same path maximum applies in reverse,

ΔEr=ΔEfΔErxn.\Delta E_{\mathrm r}^{\ddagger} =\Delta E_{\mathrm f}^{\ddagger}-\Delta E_{\mathrm{rxn}}.

These quantities are model-dependent electronic energies. Activation enthalpies and free energies require additional vibrational, thermal, standard-state, and sometimes tunneling treatment. The same level of theory, charge state, spin treatment, cell, and numerical settings must be maintained across every energy entering a comparison.

A zero barrier can be physically meaningful, or it can indicate that the band descended into an intermediate, the endpoints are not minima, the model stabilizes unintended bonding, or the maximum is one of the fixed endpoints. Inspect the structures and energy ordering before assigning an interpretation.

The stored reaction_barrier.coordinate is the normalized image index i/(N1)i/(N-1). It provides a convenient plotting axis but is not cumulative Cartesian distance, mass-weighted distance, or an intrinsic reaction coordinate.

NEB and intrinsic reaction coordinates

NEB and an intrinsic reaction coordinate (IRC) answer related but different questions. NEB begins from two supplied endpoints and relaxes a discrete band between them. An IRC begins from a verified transition state, displaces along its imaginary Hessian mode, and follows the mass-weighted steepest-descent path in both directions. The IRC is therefore a post-transition-state connectivity test rather than an endpoint interpolation method. Tako does not currently expose an IRC operation; use saddle refinement and mode inspection without describing the resulting NEB band as an IRC.

Minimum-energy paths as variational objects

The condition V=0\nabla V_\perp=0 defines a minimum-energy path geometrically, but it is useful to state what is and is not minimized. A path is not obtained by minimizing the arithmetic sum of its image energies without qualification. Such an objective would drive every unconstrained image into a nearby minimum and destroy the representation of the connecting route. Instead, an MEP is a curve whose local tangent is parallel to the potential-energy gradient wherever that gradient is nonzero. Equivalently,

V ⁣(R(s))=λ(s)τ^(s),\nabla V\!\left(\mathbf R(s)\right) =\lambda(s)\,\hat{\boldsymbol\tau}(s),

where λ(s)\lambda(s) is a scalar and τ^(s)\hat{\boldsymbol\tau}(s) is the unit path tangent. The equation leaves motion along the curve unconstrained while requiring stationarity in every direction normal to it. This local condition explains both the force projection in NEB and the role of the springs: the physical force determines the curve in transverse directions, while an artificial parametrization mechanism distributes samples along it.

The word “minimum” refers to energy with respect to transverse variations of the path, not necessarily to a global comparison among all topologically distinct connections. Given two endpoints, several locally stable paths can satisfy the perpendicular-force condition. They may pass on different sides of a steric obstacle, break different bonds, visit different intermediates, or differ in collective surface motion. An NEB calculation normally converges to the basin of attraction of the initialized band in path space. It does not enumerate those basins or certify which converged path has the smallest maximum energy.

The continuous path also has a parametrization freedom. If ss is replaced by any smooth monotonic function q(s)q(s), the geometric curve is unchanged even though the speed dR/dq\lVert d\mathbf R/dq\rVert changes. A discrete band does not enjoy this freedom automatically because image positions determine both resolution and tangent estimates. Conventional springs, spline parametrization, and string reparametrization are different numerical responses to the same issue: they attempt to keep the discrete samples useful without changing the intended geometric curve.

This distinction matters when comparing outputs. A normalized image index is a plotting coordinate, not a physical reaction coordinate. Cartesian arc length is more geometric but depends on alignment, periodic unwrapping, and the weighting assigned to different atoms. Mass-weighted distance is appropriate for dynamical constructions such as an IRC but is not the coordinate stored by Tako’s barrier payload. Bond distances, coordination numbers, or collective variables may be chemically illuminating, yet each is a projection of the band rather than the definition of the MEP.

Path topology and stationary points

A first-order saddle is a stationary point with exactly one negative Hessian eigenvalue — the classification, its zero-mode caveats, and vibrational verification are developed in Transition States and the Dimer Method. For the band, the important statement is weaker: an image that is the maximum sampled energy along a band is not thereby a stationary point of the full surface. On an exact MEP through an isolated first-order saddle, the saddle is a maximum along the path and a minimum perpendicular to it; a finite band only approximates that construction. Ordinary NEB may place the true maximum between images. CI-NEB improves the highest image, but its direction is still defined by neighboring images. A small projected NEB residual therefore does not prove that the full Cartesian gradient is small enough for a transition-state claim, nor does it determine the Hessian index.

Flat ridges and bifurcations complicate the simple picture. If the Hessian has a very small eigenvalue transverse to the path, modest changes in the model or initialization can redirect the band into another valley. Near a valley-ridge inflection, a single saddle-centered path can split toward more than one product basin. A band may also contain an intermediate minimum and two maxima rather than one elementary barrier. In these situations a single scalar “barrier” compresses important topology. The image structures and the complete energy profile must be retained, and separate path segments may be needed to characterize distinct elementary steps.

Endpoint stationarity is similarly important. NEB fixes the endpoint coordinates. If an endpoint has a significant force under the selected calculator, its energy includes strain that the band cannot remove there. The largest image energy relative to that strained reference may not represent a barrier from a metastable state. Endpoint relaxation is therefore normally part of a defensible protocol, subject to the caveat that weakly bound or symmetry-related endpoints can collapse into the same basin. In that case the scientific definition of the states, not merely the optimizer settings, must be reconsidered.

Discretization, parametrization, and error sources

A band with NN images replaces an infinite-dimensional curve by 3M(N2)3M(N-2) adjustable Cartesian coordinates for a system of MM atoms, before constraints are applied. This discretization introduces several coupled errors. The tangent is estimated from finite differences between neighboring images. The maximum energy is sampled rather than known continuously. Curvature between images is unresolved. Image spacing is controlled by an artificial numerical mechanism. Optimizer termination adds a residual error on top of these representation errors.

Increasing NN reduces the distance between neighboring images when the path length is fixed, but it does not guarantee monotonic improvement. More degrees of freedom can slow optimization, expose additional roughness in the calculator, or allow different parts of the band to enter competing valleys. A meaningful image-count study compares paths after convergence and asks whether the same atomic rearrangement, intermediate structure, and saddle region are represented. If a calculation with 15 images finds a new intermediate absent from a seven-image band, the two results are not simply different numerical resolutions of one scalar barrier.

The highest sampled energy in ordinary NEB generally underestimates the maximum of a smooth curve when the peak lies between images. CI-NEB reduces that sampling error only if the climbing image is in the correct saddle basin and the neighboring images define a reliable tangent. Conversely, a poorly initialized climbing image can be pushed to an irrelevant high-energy region. Convergence of barrier height with image count should therefore be assessed together with the climbing-image geometry, adjacent-image separation, and full-force behavior after independent saddle refinement.

Discretization error and optimizer error must be separated in reporting. A tight force threshold on seven badly spaced images does not compensate for missing spatial resolution. Likewise, a dense band stopped at a large residual is not a converged high-resolution path. A practical hierarchy is to establish a chemically plausible seed, converge a moderate band, refine the image count around regions of high curvature or multiple events, and then test tighter residual thresholds. Each stage answers a different question.

Endpoint preparation and state definition

The endpoints determine both the energy reference and the topological question posed to NEB. Ideally each is a local minimum under the same calculator, charge, spin treatment, cell, constraints, and numerical settings used for the band. Consistent relaxation removes arbitrary strain and makes the reported forward barrier relative to a defined metastable state. The final endpoint supplies the reaction energy and, together with the first endpoint, determines the reverse barrier when one common path maximum is appropriate.

This rule has exceptions that must be treated scientifically rather than procedurally. Independently relaxing weakly associated fragments can lead both endpoints into the same aggregate. Symmetry-related exchange endpoints can collapse to one representation. A surface adsorbate may diffuse away from a constructed site if the chosen model does not support that site as a minimum. In such cases blindly enforcing endpoint optimization erases the process under study. Possible responses include defining more separated fragments, using physically justified constraints, choosing a model that supports the intended states, or accepting that the proposed states are not metastable on that potential-energy surface.

Constraints must describe the same approximation across the path. Tako copies constraints from the initial structure into every image. It does not merge a separate constraint definition from the final endpoint. A fixed atom is therefore governed by its initial constraint and reference position throughout the band. If that atom occupies a different intended position in the supplied final structure, the band construction and constraint conflict. Partial Cartesian constraints can be useful for substrate layers or symmetry planes, but they also change the coordinate space in which the MEP condition is satisfied.

Atom mapping is part of endpoint preparation. Equal formulas are insufficient, and element-by-index validation catches only the simplest mismatch. Exchanging two hydrogens, permuting equivalent surface atoms, or importing independently sorted structures can define a long or self-crossing path without triggering an error. Useful mapping checks compare local environments, fragment identities, and displacements, not only element symbols. For periodic structures, this includes choosing consistent lattice translations for each mapped atom.

Electronic state consistency is equally important. A bond-breaking path can traverse geometries with different multireference character, charge localization, or spin ordering even when the endpoint labels are simple. Tako uses one calculator definition for all images; this consistency is necessary for a continuous profile but does not guarantee that the model follows the physically relevant electronic state. Discontinuous SCF solutions, unsupported elements, or model extrapolation can appear as noisy forces, abrupt energy steps, or convergence failure. A reaction path is an especially demanding validation set because it intentionally samples configurations unlike equilibrium training or parametrization data.

Competing mechanisms and calculator validation

A single converged band establishes only a local numerical solution. Mechanistic comparison requires deliberately different seeds. For a substitution, one might distinguish front-side and back-side approach. For surface diffusion, paths can pass through bridge, hollow, or exchange configurations. A conformational change can rotate different torsions in different orders. Each hypothesis should be encoded in its own atom-consistent seed, converged independently, and compared using the same calculator and endpoint convention.

Barrier height alone is not enough to decide whether two bands represent the same mechanism. Compare bond changes, coordination, fragment orientation, intermediates, path length, and saddle geometry. Two profiles with nearly identical maxima may correspond to different chemistry, while two image counts can yield different maxima because one resolves an intermediate and the other does not. If paths share an intermediate, splitting the problem into elementary segments often produces clearer convergence and transition-state verification.

Calculator validation should test quantities relevant to the path. Endpoint energies and geometries are necessary but weak evidence for a barrier. Where feasible, compare selected path images or the refined saddle with a more defensible electronic-structure method, check reaction energies against known thermochemistry, and examine whether bond-breaking forces are smooth. For an MLIP, assess whether the path leaves its training domain; model disagreement or uncertainty estimates are valuable when available. For semiempirical models, known failures in charge transfer, radicals, unusual coordination, or spin energetics should be treated as model uncertainty rather than NEB error.

Energy conservation is not a relevant check because NEB is an optimization, not dynamics. Smoothness is useful but not decisive: an inaccurate model can produce a beautifully smooth wrong surface. Conversely, an electronically difficult but physically meaningful region can generate noisy convergence. Independent calculations at fixed images help distinguish optimizer behavior from calculator instability. Changing the band optimizer cannot repair inconsistent energies and forces returned by the model.

Barrier uncertainty and interpretation

The reported forward barrier is a deterministic summary of the sampled electronic energies, but its scientific uncertainty has several components. Model error includes the electronic-structure approximation, parametrization, basis or numerical settings, dispersion and solvation treatment, and electronic-state selection. Path error includes endpoint definition, atom mapping, mechanism choice, discretization, tangent formulation, and incomplete convergence. Thermodynamic corrections add vibrational, temperature, entropy, standard-state, and possibly tunneling approximations.

These components should not be collapsed into the final NEB residual. A residual of 0.010.01 eV Å1^{-1} does not imply a barrier uncertainty of 0.010.01 eV. Nor does agreement to two decimal places between image counts establish model accuracy. A defensible uncertainty discussion states which dimensions were varied: for example, image count and residual threshold for numerical sensitivity, alternative seeds for mechanistic sensitivity, and a second calculator for model sensitivity.

Electronic barriers are also direction-dependent. If the initial and final energies differ, the forward and reverse barriers differ even though they share a maximum. When the path contains multiple maxima or intermediates, using the global maximum to infer every elementary rate is inappropriate. Each metastable basin and adjacent saddle define its own elementary barrier. Free-energy kinetics further require a consistent treatment of partition functions and standard states; the NEB electronic profile alone does not supply those quantities.

The saddle-image index is a discrete label, not an uncertainty measure. A change from image 3 to image 4 after increasing NN may simply reflect denser sampling of the same geometry. Compare aligned structures and local coordinates rather than indices. Likewise, the normalized coordinate has no kinetic interpretation. It is suitable for plotting image order but not for calculating a transmission coefficient or assigning a physical time.

From band maximum to transition state

CI-NEB improves the spatial location of the path maximum, but the climbing image remains constrained by the band tangent and the finite image representation. Treat it as a high-quality starting geometry for a dedicated saddle refinement, not as a verified transition state.

The refinement itself — dimer translation, the exactly-one-imaginary-mode requirement, mode inspection, and forward/backward connectivity evidence — is developed in Transition States and the Dimer Method. Within the band’s own evidence, the NEB neighbors on either side of the candidate provide contextual support for the tangent direction, but they do not replace the local Hessian and downhill connectivity checks. If the refined structure has no imaginary mode it is probably a minimum; multiple substantial imaginary modes indicate a higher-order saddle or incomplete refinement; a sole mode pointing to unintended products means the band located a different mechanism. These outcomes are scientific results about the chosen surface and initialization, not reasons to relabel the climbing image as verified.

Computational cost and scaling

Each NEB force evaluation requires calculator work on the active images, with the two endpoints sometimes benefiting from cached energies while the interior images dominate repeated cost. For NN total images, the leading work per band step is approximately proportional to N2N-2 calculator force evaluations. Wall time also depends on whether those evaluations are scheduled serially or concurrently, calculator initialization cost, memory pressure, and electronic convergence at difficult geometries.

Tako creates an independent calculator instance for each image. This prevents mutable state from leaking between images but increases memory use, particularly for MLIP runtime sessions or electronic-structure state. Increasing image count therefore affects both computation and memory. IDPP adds a comparatively inexpensive fictitious optimization, while an optional vibration calculation can be much more expensive because finite differences require many force evaluations for the selected structure.

The two CI-NEB phases affect cost accounting. The configured final step limit is not the total possible number of optimizer steps when pre-climbing is enabled. Failed calculator evaluations can terminate before either budget is exhausted. Dynamic relaxation may reduce coordinate updates but does not make image-count scaling disappear.

Performance choices must preserve scientific resolution. Reducing image count is reasonable after demonstrating that the same mechanism and saddle are retained. Loosening the residual is reasonable for scouting alternative seeds, but final comparisons should use consistent criteria. A cheap calculator can generate a seed for a higher-level calculation only if the handoff preserves atom mapping and the higher-level model does not redirect the mechanism unexpectedly.

Reproducibility and reporting

The shared record skeleton and reporting standard are defined in Reproducibility and Reporting; the NEB-specific requirements follow. A reproducible NEB record contains both endpoint structures, their atom order and cell, constraints, calculator identity and parameters, charge and spin treatment, image count, interpolation method and controls, NEB formulation, optimizer, spring constant, residual threshold, maximum displacement, pre-climb and final budgets, dynamic-relaxation settings when used, software version, convergence state, final residual, and the full ordered band. Reporting only the scalar barrier discards the evidence needed to assess the mechanism.

The energy profile should state whether energies are absolute or relative and identify the reference endpoint. Include reaction energy alongside the forward barrier so the reverse relationship is visible. State explicitly that values are electronic energies unless vibrational or thermodynamic corrections were performed. If a band is unconverged, label every derived barrier as provisional and preserve the residual and step history.

Sensitivity calculations should be organized as comparable cases rather than overwritten reruns. Record which quantity changed and keep all other settings fixed. Compare structures as well as numbers. When two seeds converge to distinct mechanisms, report them as alternatives; do not average their barriers. When a model comparison changes endpoint ordering or creates an intermediate, explain that the potential-energy surfaces are qualitatively different.

Figures should show enough images or structural descriptors to establish the atomic event, not only a smooth line through energies. A profile plotted against normalized image index should be labeled accordingly. Interpolation-only frames should never be mixed into an electronic-energy plot as though their zero placeholders were calculated energies. For a claimed transition state, report the refined geometry, gradient criterion, imaginary frequency or curvature evidence, mode interpretation, and connectivity test separately from the NEB convergence statement.

Worked reasoning patterns

Consider first a compact intramolecular rearrangement with relaxed endpoints and one expected bond exchange. A seven-image IDPP seed is inspected and contains no overlaps. A converged CI-NEB band has a stable highest image, smooth neighboring structures, and a maximum that changes little at 11 images. The evidence supports a resolved local path, but not yet a transition-state claim. Dimer refinement and a single chemically aligned imaginary mode complete the local saddle argument; comparison with another calculator addresses model sensitivity.

Now consider surface diffusion across a periodic boundary. Independent wrapping sends the adsorbate across the full cell, producing a long seed with distorted substrate contacts. IDPP cannot repair the incorrect lattice branch because the endpoint mapping defines that displacement. Remapping the final adsorbate by a lattice translation creates the intended short hop. Fixed bottom-layer atoms are defensible if declared, but fixing the moving surface environment could create an artificially narrow barrier. Here endpoint mapping and constraint sensitivity dominate before force convergence matters.

Finally, consider a profile with an interior minimum below both endpoints and two shoulders. Treating the global maximum as one elementary barrier hides the likely intermediate. Additional images may confirm two separate saddle regions. The scientifically useful response is to characterize the intermediate, split the path, and refine each candidate saddle. If another calculator removes the intermediate entirely, the difference is model dependence, not an image-count correction. These patterns illustrate why path structures, topology, and model validation carry more information than a single converged flag.

User interfaces and Script API

This chapter owns the scientific and numerical description. The operational instructions are deliberately separate:

References

  1. G. Mills, H. Jónsson, and G. K. Schenter, “Reversible work transition state theory: application to dissociative adsorption of hydrogen,” Surface Science 324, 305–337 (1995), doi:10.1016/0039-6028(94)00731-4.
  2. H. Jónsson, G. Mills, and K. W. Jacobsen, “Nudged elastic band method for finding minimum energy paths of transitions,” in Classical and Quantum Dynamics in Condensed Phase Simulations, 385–404 (1998), doi:10.1142/9789812839664_0016.
  3. G. Henkelman and H. Jónsson, “Improved tangent estimate in the nudged elastic band method for finding minimum energy paths and saddle points,” Journal of Chemical Physics 113, 9978–9985 (2000), doi:10.1063/1.1323224.
  4. G. Henkelman, B. P. Uberuaga, and H. Jónsson, “A climbing image nudged elastic band method for finding saddle points and minimum energy paths,” Journal of Chemical Physics 113, 9901–9904 (2000), doi:10.1063/1.1329672.
  5. S. Smidstrup, A. Pedersen, K. Stokbro, and H. Jónsson, “Improved initial guess for minimum energy path calculations,” Journal of Chemical Physics 140, 214106 (2014), doi:10.1063/1.4878664.
  6. D. Sheppard, R. Terrell, and G. Henkelman, “Optimization methods for finding minimum energy paths,” Journal of Chemical Physics 128, 134106 (2008), doi:10.1063/1.2841941.
  7. W. E, W. Ren, and E. Vanden-Eijnden, “String method for the study of rare events,” Physical Review B 66, 052301 (2002), doi:10.1103/PhysRevB.66.052301.
  8. A. Heyden, A. T. Bell, and F. J. Keil, “Efficient methods for finding transition states in chemical reactions: comparison of improved dimer method and partitioned rational function optimization method,” Journal of Chemical Physics 123, 224101 (2005), doi:10.1063/1.2104507.
  9. K. Fukui, “The path of chemical reactions—the IRC approach,” Accounts of Chemical Research 14, 363–368 (1981), doi:10.1021/ar00072a001.