Generalized Solid-State NEB Calculations

• Generalized solid-state NEB calculations are methods that find minimum-energy paths and saddle points in crystalline systems by incorporating both atomic and cell variables.
• They utilize advanced initialization techniques such as adaptive semi-rigid-body approximations, hybrid lattice–SLERP, and symmetry-adaptive perturbations to ensure realistic and convergent paths.
• The approach accelerates computation using surrogates like Gaussian processes and machine-learned force fields, significantly reducing ab initio calculation costs while maintaining accuracy.

A generalized solid-state nudged elastic band (NEB) calculation refers to the family of algorithms, workflows, and supporting models developed to efficiently and robustly search for minimum-energy paths (MEPs) and saddle points associated with transition states in crystalline and extended solid-state systems. These frameworks address practical and theoretical challenges specific to solids, including periodicity, cell-shape degrees of freedom, collective mechanisms, and the high computational cost of ab initio potentials. Over the last decade, methodological advances have established a comprehensive toolkit: semi-rigid-body initializations for rapid convergence, variable-cell algorithms, symmetry-adapted perturbations for pathway discovery, surrogate and machine-learned force fields for cost reduction, and high-dimensional optimizers for scalable application.

1. Theoretical Foundations and Force Projection in Solid-State NEB

The NEB method formulates the MEP search as a force-projection problem in configuration space. For solid-state systems, the relevant coordinate manifold extends naturally from atomic positions $\{\mathbf{R}_i\}$ to include cell-shape/volume and potentially additional order parameters such as magnetization. The band of images $\{(\mathbf{R}_i, \mathbf{h}_i)\}_{i=0}^{N}$, with $\mathbf{h}_i$ the cell matrix, is evolved by separating the total force on each image into physical forces perpendicular to the path and artificial (spring) forces along the band tangent. The projected forces for image $i$ are

$$\mathbf{F}_i^{\mathrm{NEB}} = \mathbf{F}_i^\perp + \mathbf{F}_i^\mathrm{spring,\parallel}$$

where

$$\mathbf{F}_i^\perp = -\nabla_{\mathbf{R},\mathbf{h}} V(\mathbf{R}_i, \mathbf{h}_i) + [\nabla_{\mathbf{R},\mathbf{h}} V(\mathbf{R}_i, \mathbf{h}_i) \cdot \hat{\tau}_i] \hat{\tau}_i,$$

and the elastic force,

$$\mathbf{F}_i^\mathrm{spring,\parallel} = k\left( \|\mathbf{R}_{i+1} - \mathbf{R}_i\| - \|\mathbf{R}_i - \mathbf{R}_{i-1}\| \right) \hat{\tau}_i,$$

maintains image spacing. For methods permitting variable cell or non-conserved order parameters, appropriate generalizations of $\nabla$ and $\hat{\tau}_i$ are adopted, e.g., block vectors and cell-filtered coordinates (Goncharova et al., 2024, Jia et al., 19 Oct 2025, Zarkevich et al., 2015).

2. Initialization Schemes: Semi-Rigid Body, Hybrid, and Symmetry-Adaptive Strategies

Initial guess construction is pivotal for convergence efficiency. Standard linear interpolation is prone to generating unrealistic or high-energy intermediate images in the solid state, especially with large unit cells or rigid molecular motifs.

Adaptive Semi-Rigid Body Approximation (ASRBA)

ASRBA initializes images by imposing near-invariance of local bond lengths along the path. Each atom and its nearest neighbors are grouped into "semi-rigid bodies," and adaptive atomic radii $\{R_i\}$ are determined by solving

$$R_i + R_j = b_{ij} = \tfrac12( d_{ij}^A + d_{ij}^B ),$$

over all bonded pairs in the endpoint structures. Inter-image forces penalize deviation from these equilibrium distances, ensuring physically plausible intermediates and suppressing artificial deformations at the initialization stage (Cai et al., 2022).

Hybrid Lattice–SLERP (Spherical Linear Interpolation) Interpolation

For molecular crystals and polymorphic transitions, initial images are generated by linear interpolation of lattice matrices and SLERP of molecular orientations encoded as quaternions. Interfragment bond distances and centers of mass are interpolated linearly; intermolecular collisions are detected and resolved by iterative cell expansion and redistribution. This produces collision-free, continuous starting paths for NEB (Goncharova et al., 2024).

Distortion Symmetry and Symmetry-Adaptive Perturbations

Symmetry constraints of the initial path can artificially restrict the accessible set of MEPs. The distortion symmetry group $G$ of a path is computed, and symmetry-adapted perturbations constructed via projection operators onto specific irreducible representations systematically lower path symmetry and expose otherwise hidden mechanisms (Munro et al., 2018). This systematic perturbative approach is essential for thorough exploration of transition pathways in solids with multiple competing mechanisms.

3. Acceleration via Surrogates: Gaussian Processes and Machine-Learned Force Fields

Ab initio NEB calculations incur high computational cost per image. Surrogate models—constructed via Gaussian process regression (GPR) or neural-network force fields—provide efficient, uncertainty-aware interpolations of the potential energy surface (PES) and forces.

GPR-Accelerated NEB

Atomic coordinates (and, as needed, cell parameters and spins) are encoded in descriptor vectors $x$. A GPR surrogate $E^\mathrm{GP}(x)$ is trained on observed energies and forces. At each NEB step, images are updated on the GP-surface, and new ab initio evaluations are triggered if prediction uncertainty $\sigma_F(x_i)$ exceeds a user-set threshold. The typical workflow leads to 3–10$\times$ reduction in first-principles calculations while retaining high-fidelity barriers (errors $<$ 0.05 eV) (Onyango et al., 9 Apr 2025, Koistinen et al., 2017, Koistinen et al., 2017).

$$\mu(x_*) = k_*^T [K + \sigma_n^2 I ]^{-1} y, \quad \sigma^2(x_*) = k(x_*, x_*) - k_*^T [K + \sigma_n^2 I ]^{-1} k_*$$

Machine-Learned Force Fields (MLFFs)

Equivariant neural-network-based MLFFs, such as SO3krates, are trained on DFT energies and forces sampled along initial MEPs (e.g., from NPT molecular dynamics and key images). The MLFF predicts energies, forces, and stresses for new images, with training protocols targeting MAE $<$ 0.4 kJ/mol per molecular unit and force MAE $\sim$1–2 kJ/mol/Å. This enables efficient global geometry optimization of large images and rapid NEB path refinement (Goncharova et al., 2024).

4. Specialization and Extension: Variable-Cell, Symmetry, and Order Parameter Methods

Variable-Cell and Cell-Filtered NEB

The solid-state NEB must accommodate cell degrees of freedom, especially for phase transitions and reconstructive transformations. Advanced schemes such as CFNEB express the collective coordinates of each image as $(\tilde{R}_i, B)$, where $\tilde{R}_i=R_i D^{-1}$ eliminates rigid cell rotations, and $B=D D^T$ (left Cauchy–Green tensor) represents lattice deformations. NEB force projections are performed in the combined $(3N+6)$- or $(3N+9)$-dimensional metric, with conjugate forces evaluated via virtual work and virial stress. Adaptive image insertion/deletion dynamically maintains image density along high-curvature regions of the path, and the scheme achieves GPU-accelerated scaling to $10^{5}$ atoms per band (Jia et al., 19 Oct 2025).

NEB with Non-Conserved Order Parameters

For magneto-structural transformations (e.g., bcc $\rightarrow$ hcp Fe) involving changes in a non-conserved order parameter $m$ (magnetization), NEB optimizes an enthalpy function $H(\{R\}, V, m)$, constraining $m$ within each image and allowing for discontinuous transitions (cusps) at the saddle point. Two climbing-image variants are used to bracket the TS from both sides, avoiding artifacts from false intermediate states (Zarkevich et al., 2015).

Reduced-Dimension and Constrained Relaxation NEB

For transitions governed by a few collective variables (e.g., defect migration), the NEB force projection can be restricted to a reduced set of reaction coordinates $\mathbf{q}$, with all other degrees of freedom relaxed fully at each NEB step. This mitigates “image clustering” artifacts and ensures uniform progression along the physical reaction coordinate (Gröger et al., 2011).

5. Numerical Optimization, Scaling, and Automation

Dynamic and Preconditioned Optimization

Optimization efficiency in high-dimensional NEB landscapes is improved through algebraic preconditioning (local Hessian or bond-stiffness preconditioners), P-norm force projections, and adaptive timestep selection (ode12r). These enhancements yield acceleration factors of 2–10× in force evaluations for MEP convergence—particularly on stiff solid-state landscapes (Makri et al., 2018).

Dynamic NEB (dyNEB) assigns tighter convergence thresholds to saddle-adjacent images and relaxed criteria to tail images, adaptively focusing effort where needed, reducing wall time and force-call count across diverse solid-state and electrochemical systems (Lindgren et al., 2019).

Automated Band Construction and Parallelization

Frameworks such as PASTA and AutoNEB automate image initialization, dynamic insertion of new images in high-curvature or high-energy regions, and flexible parallelization across images. These features enable systematic, high-throughput exploration of reaction barriers and phase transformations in solid-state and interface systems (Kundu et al., 2018).

6. Performance Benchmarks, Limitations, and Best Practices

Algorithmic innovations deliver significant computational acceleration and improved robustness. Performance comparisons demonstrate that GPR-NEB, MLFF-accelerated NEB, and preconditioned/dyNEB methods achieve order-of-magnitude speedups compared to classical NEB, with mean errors in predicted barriers typically $<$0.05 eV (GPR) or $<$0.4 kJ/mol (MLFF). Adaptive and symmetry-adapted initialization methods guarantee rapid convergence and correct identification of relevant pathways, provided the underlying physical assumptions (bond invariance, symmetry, or reaction coordinate completeness) are satisfied (Cai et al., 2022, Goncharova et al., 2024, Onyango et al., 9 Apr 2025, Zarkevich et al., 2014, Munro et al., 2018).

Limitations include the potential combinatorial complexity in large or low-symmetry systems, dependency of surrogate accuracy on descriptor completeness and training data, and pitfalls if reaction coordinates fail to capture all relevant transitions. In high-dimensional settings, dimensionality reduction techniques (PCA, ARD), sparse GPR, and dynamical image management are essential.

Parameter recommendations, as synthesized from recent implementations, include:

• Spring constants $k \sim$ 2–5 eV/Å$^2$ (higher for cell or collective DOF)
• Force convergence thresholds: $<$0.05 eV/Å for images, $<$0.01 eV/Å for TS
• GPR force-uncertainty threshold: 0.05–0.1 eV/Å
• Image count: 8–32 after DOF reduction, 16–64 for full-dimensional cell-variable bands
• MLFF training: $\sim$5000–10,000 structures sampled along path and far-from equilibrium ensembles

These best practices together establish the modern, scalable, and systematically improvable paradigm for solid-state NEB calculations across bulk, surface, defect, and phase-transition phenomena.