Spectral Neighbor Analysis Potential (SNAP)

• SNAP is a machine-learning interatomic potential that decomposes total energy into local contributions using rotationally invariant bispectrum components derived from a 3-sphere mapping.
• It employs both linear and quadratic regression models to fit energies, forces, and stresses, enabling scalable molecular dynamics simulations with near quantum-mechanical accuracy.
• Extensions like EME-SNAP and eSNAP enhance its applicability by incorporating multi-element descriptors and electrostatic corrections to address complex materials and ionic systems.

Searching arXiv for SNAP foundational and extension papers to ground the article in current and original sources. arxiv_search(query="Spectral Neighbor Analysis Potential SNAP bispectrum 2014 Thompson quadratic SNAP Wood Thompson 2017 explicit multi-element SNAP 2020 electrostatic SNAP lithium nitride", max_results=10) Spectral Neighbor Analysis Potential (SNAP) is a machine-learning interatomic potential in which the total energy is decomposed into local atomic contributions, and each local environment is represented by rotationally invariant bispectrum components derived from a spectral expansion of an atom-centered neighbor density on the 3-sphere. In its original formulation, SNAP was introduced as a short-ranged, linear-scaling potential for solids and liquids that aims to approach quantum-mechanical accuracy while remaining suitable for large-scale molecular dynamics; subsequent work extended the same descriptor machinery to quadratic regressions, explicit multi-element descriptors, electrostatic corrections for ionic systems, and spectral-neighbour representations for vector fields including spin (Thompson et al., 2014).

1. Historical emergence and conceptual role

SNAP was introduced as a new interatomic potential for solids and liquids in order to interpolate between two established limits of atomistic simulation: conventional empirical potentials, which are computationally efficient but often too restrictive, and direct quantum molecular dynamics, which is accurate but too expensive for large systems and long trajectories. Its original motivating use case was tantalum, especially across crystals, defects, liquid structure, and screw-dislocation motion in bcc Ta (Thompson et al., 2014).

A defining conceptual feature of SNAP is that it separates descriptor construction from regression form. The descriptor is a rotationally invariant spectral fingerprint of the local atomic environment, while the original regression model is linear in those descriptors. This distinguishes SNAP from Gaussian Approximation Potential (GAP), which uses the same bispectrum representation but models local energy with Gaussian-process regression rather than an explicit linear expansion (Thompson et al., 2014).

Later work positioned SNAP more broadly as a DFT-trained interatomic potential capable of carrying ab initio fidelity into classical molecular dynamics at scales relevant to defects, phase boundaries, nucleation, and extreme-pressure thermodynamics. In the carbon literature, for example, SNAP is presented as a route from small-cell QMD/AIMD to simulations of millions of atoms, with applications to melting curves, shock Hugoniots, and phase boundaries under multi-TPa conditions (Willman et al., 2022).

2. Descriptor construction and invariance structure

For an atom $i$, SNAP begins from a local neighbor density defined within a cutoff sphere. In the original formulation this density is written as

$$\rho_i ({\bf r}) = \delta({\bf r}) + \sum_{r_{ii'} < R_{cut}} f_c(r_{ii'})\, w_{i'} \,\delta({\bf r}-{\bf r}_{ii'}),$$

where the central atom has unit weight, $w_{i'}$ distinguishes species, and $f_c(r)$ smoothly turns off neighbor contributions at the cutoff (Thompson et al., 2014).

The 3D neighborhood is then mapped onto a subset of the 3-sphere by promoting the radial coordinate to a fourth angular coordinate,

$\theta_0 = \theta_0^{max}\frac{r}{R_{cut}},$

and the density is expanded in 4D hyperspherical harmonics $U^j_{m,m'}(\theta_0,\theta,\phi)$: $$\rho({\bf r}) = \sum_{j=0,\frac{1}{2},\ldots}^{\infty}\sum_{m=-j}^{j}\sum_{m'=-j}^{j} u^j_{m,m'} U^j_{m,m'}(\theta_0,\theta,\phi).$$ The coefficients $u^j_{m,m'}$ are not themselves invariant under rotation, so SNAP forms rotationally invariant bispectrum components as scalar triple products of the hyperspherical-harmonic coefficients (Thompson et al., 2014).

In one explicit multicomponent statement of the formalism, the local density may be written as

$$\rho_i(\vec r)=\delta(\vec r-\vec R_i)+\sum_j w_{Z_j}\,\delta(\vec r-\vec R_j)\,f_c(R_{ij}),$$

followed by the hyperspherical-harmonic expansion

$$\rho_i(\vec r)=\sum_{J=0}^{\infty}\sum_{m,m'=-J}^{J} c_{m',m}^{J}\,U^{J}_{m',m}(\theta,\phi,\theta_0),$$

and the bispectrum invariant

$$\rho_i ({\bf r}) = \delta({\bf r}) + \sum_{r_{ii'} < R_{cut}} f_c(r_{ii'})\, w_{i'} \,\delta({\bf r}-{\bf r}_{ii'}),$$0

In practice the expansion is truncated at $$\rho_i ({\bf r}) = \delta({\bf r}) + \sum_{r_{ii'} < R_{cut}} f_c(r_{ii'})\, w_{i'} \,\delta({\bf r}-{\bf r}_{ii'}),$$1, and the retained bispectrum components form the descriptor vector $$\rho_i ({\bf r}) = \delta({\bf r}) + \sum_{r_{ii'} < R_{cut}} f_c(r_{ii'})\, w_{i'} \,\delta({\bf r}-{\bf r}_{ii'}),$$2 (Minotakis et al., 2023).

The resulting descriptor set is explicitly presented as invariant under translation, rotation, and permutation of identical atoms, while retaining angular and radial information in a compact local representation. The original SNAP paper also emphasizes a bispectrum symmetry,

$$\rho_i ({\bf r}) = \delta({\bf r}) + \sum_{r_{ii'} < R_{cut}} f_c(r_{ii'})\, w_{i'} \,\delta({\bf r}-{\bf r}_{ii'}),$$3

which reduces redundant components and lowers the asymptotic cost of bispectrum-derivative evaluation from $$\rho_i ({\bf r}) = \delta({\bf r}) + \sum_{r_{ii'} < R_{cut}} f_c(r_{ii'})\, w_{i'} \,\delta({\bf r}-{\bf r}_{ii'}),$$4 to $$\rho_i ({\bf r}) = \delta({\bf r}) + \sum_{r_{ii'} < R_{cut}} f_c(r_{ii'})\, w_{i'} \,\delta({\bf r}-{\bf r}_{ii'}),$$5 (Thompson et al., 2014).

3. Linear SNAP, quadratic SNAP, and the regression problem

The original SNAP energy model is linear in the bispectrum descriptors. In the elemental form,

$$\rho_i ({\bf r}) = \delta({\bf r}) + \sum_{r_{ii'} < R_{cut}} f_c(r_{ii'})\, w_{i'} \,\delta({\bf r}-{\bf r}_{ii'}),$$6

so that

$$\rho_i ({\bf r}) = \delta({\bf r}) + \sum_{r_{ii'} < R_{cut}} f_c(r_{ii'})\, w_{i'} \,\delta({\bf r}-{\bf r}_{ii'}),$$7

Because the energy is linear in $$\rho_i ({\bf r}) = \delta({\bf r}) + \sum_{r_{ii'} < R_{cut}} f_c(r_{ii'})\, w_{i'} \,\delta({\bf r}-{\bf r}_{ii'}),$$8, the forces and stress are also linear in the same coefficients: $$\rho_i ({\bf r}) = \delta({\bf r}) + \sum_{r_{ii'} < R_{cut}} f_c(r_{ii'})\, w_{i'} \,\delta({\bf r}-{\bf r}_{ii'}),$$9

$w_{i'}$0

Training is therefore posed as a weighted least-squares regression against QM energies, forces, and stresses, assembled as a linear system $w_{i'}$1 and solved by QR factorization in the foundational work (Thompson et al., 2014).

A later extension introduced quadratic SNAP, in which the atomic energy includes pairwise products of bispectrum components: $w_{i'}$2 This preserves the same descriptor space but augments the regression basis. The quadratic form is explicitly described as structurally similar to the embedded atom method (EAM), with bispectrum components acting as counterparts to the two-body density functions in EAM. In this interpretation, linear SNAP is a four-body potential, whereas quadratic SNAP can involve up to seven distinct atoms (Wood et al., 2017).

The increase in expressivity is accompanied by a large increase in parameter count. For example, the quadratic-SNAP paper reports that at $w_{i'}$3, $w_{i'}$4, the linear model has $w_{i'}$5 coefficients while the quadratic form has $w_{i'}$6 quadratic coefficients; at $w_{i'}$7, $w_{i'}$8, the quadratic coefficient count rises to $w_{i'}$9. The same work reports that the quadratic model produces about an order-of-magnitude lower training energy error than the linear form for tantalum at high $f_c(r)$0, but also requires substantially more training data to prevent overfitting (Wood et al., 2017).

Quadratic forms have subsequently been adopted in production potentials. For carbon under extreme conditions, the atomic energy was written as

$f_c(r)$1

and for cobalt the q-SNAP formulation was defined as “SNAP bispectrum descriptors + quadratic regression terms,” with the main distinction from standard SNAP stated explicitly as: SNAP is linear in descriptors, q-SNAP adds a quadratic feature expansion (Willman et al., 2022, Bideault et al., 2024).

4. Training databases, fitting workflows, and domain of applicability

SNAP fitting is inseparable from training-set design. The original tantalum work emphasized a training set spanning displaced crystals, strained cells, generalized stacking faults, surfaces, and liquid snapshots rather than only equilibrium structures (Thompson et al., 2014). Later studies made this dependence on coverage still more explicit.

For cobalt, the q-SNAP training set contains 1049 cobalt structures spanning bulk unit cells of hcp, fcc, bcc, and $f_c(r)$2 phases; strained unit cells; bulk MD snapshots; vacancies and interstitials; surfaces with adatoms, vacancies, and steps; nanoparticles; small nanoclusters; and liquid/amorphous structures. The same study states that inclusion of very small clusters (2–30 atoms) was essential because excluding them caused failure for highly undercoordinated surface and cluster configurations. Its final weighted least-squares fit used descriptor scans over $f_c(r)$3 to 6 and cutoff radii from 4 to 6 Å, with the best cobalt model obtained at band limit 8 ($f_c(r)$4) and radial cutoff 5.0 Å; the final fitting weights were energies $f_c(r)$5, forces $f_c(r)$6, and stresses $f_c(r)$7 (Bideault et al., 2024).

Carbon at extreme conditions was trained with a different but equally explicit workflow: a database of 636 structures, 124,907 unique atomic environments, and 388,077 regression equations, followed by an iterative cycle of database construction, parameter training, validation against QMD and experiment, deficiency identification, and retraining. Hyperparameters $f_c(r)$8 Å and $f_c(r)$9 were selected to balance accuracy and computational cost (Willman et al., 2022).

Other studies emphasized data-reduction and diversity selection. The molybdenum work used PCA of bispectrum features and their derivatives to inspect descriptor-space coverage and excluded a vacancy dataset judged redundant in that space (Chen et al., 2017). The niobium study used Vendi Score and D-optimal selection, trained on 550 total configurations, and then validated on elastic, thermal, and defect properties including collision cascades (Bhardwaj et al., 5 Feb 2025).

A recurrent limitation is that SNAP is reliable primarily when prediction remains interpolative in descriptor space. In binary-trained ternary screening, an ensemble of SNAPs trained only on binary phases predicted ternary energies near equilibrium reasonably well but failed for unrelaxed prototype geometries and for structural relaxation itself. For Cu–Ag–Au, the reported RMSE was 17.06 meV/atom on the ternary test set, remained 36.98 meV/atom for relaxed ternary prototypes, but rose to about 1.2 eV/atom for the same prototypes in their initial unrelaxed geometries. Feature-space analysis by PCA showed that the relaxation path traversed regions poorly represented in the training set (Minotakis et al., 2023).

The same ternary study compared three extrapolation metrics and found that the ensemble standard deviation $\theta_0 = \theta_0^{max}\frac{r}{R_{cut}},$0 correlated best with actual error: in the Cu–Ag–Au test set the reported correlation coefficients were $\theta_0 = \theta_0^{max}\frac{r}{R_{cut}},$1 for $\theta_0 = \theta_0^{max}\frac{r}{R_{cut}},$2, $\theta_0 = \theta_0^{max}\frac{r}{R_{cut}},$3 for the extrapolation grade $\theta_0 = \theta_0^{max}\frac{r}{R_{cut}},$4, and much worse for a globalized minimum Euclidean distance $\theta_0 = \theta_0^{max}\frac{r}{R_{cut}},$5. This suggests that domain-of-applicability control is an operational part of SNAP deployment rather than a secondary diagnostic (Minotakis et al., 2023).

5. Formal extensions beyond the original single-channel local model

Several extensions retain the spectral-neighbour logic of SNAP while modifying the local density, the descriptor channels, or the total-energy decomposition.

Extension	Mechanism	Representative result
Quadratic SNAP	Adds $\theta_0 = \theta_0^{max}\frac{r}{R_{cut}},$6 terms to the linear bispectrum model	Order-of-magnitude training-energy improvement for Ta; deployed for C and Co (Wood et al., 2017)
EME-SNAP	Partitions neighbor density into element-resolved partial densities	Descriptor count scales as $\theta_0 = \theta_0^{max}\frac{r}{R_{cut}},$7, force cost as $\theta_0 = \theta_0^{max}\frac{r}{R_{cut}},$8 (Cusentino et al., 2020)
eSNAP	Adds screened Ewald electrostatics to local SNAP energy	For $\theta_0 = \theta_0^{max}\frac{r}{R_{cut}},$9-$U^j_{m,m'}(\theta_0,\theta,\phi)$0, $U^j_{m,m'}(\theta_0,\theta,\phi)$1 with $U^j_{m,m'}(\theta_0,\theta,\phi)$2 (Deng et al., 2019)

The explicit multi-element extension (EME-SNAP) replaces the weighted-density treatment of chemistry by partitioning the neighbor density into partial densities for each chemical element: $U^j_{m,m'}(\theta_0,\theta,\phi)$3 with each $U^j_{m,m'}(\theta_0,\theta,\phi)$4 expanded separately in hyperspherical harmonics. The resulting explicit multi-element bispectrum is indexed by ordered element triplets,

$U^j_{m,m'}(\theta_0,\theta,\phi)$5

For $U^j_{m,m'}(\theta_0,\theta,\phi)$6 elements, the number of descriptors grows as $U^j_{m,m'}(\theta_0,\theta,\phi)$7, while the force calculation in LAMMPS scales as $U^j_{m,m'}(\theta_0,\theta,\phi)$8. In InP, this chemically resolved formulation improved defect energetics relative to weighted-density SNAP while retaining good bulk-property accuracy (Cusentino et al., 2020).

The electrostatic SNAP formulation addresses the inability of a finite-cutoff local descriptor to represent long-range Coulomb interactions in ionic systems. For $U^j_{m,m'}(\theta_0,\theta,\phi)$9-$$\rho({\bf r}) = \sum_{j=0,\frac{1}{2},\ldots}^{\infty}\sum_{m=-j}^{j}\sum_{m'=-j}^{j} u^j_{m,m'} U^j_{m,m'}(\theta_0,\theta,\phi).$$0, the total potential energy was written as

$$\rho({\bf r}) = \sum_{j=0,\frac{1}{2},\ldots}^{\infty}\sum_{m=-j}^{j}\sum_{m'=-j}^{j} u^j_{m,m'} U^j_{m,m'}(\theta_0,\theta,\phi).$$1

with forces

$$\rho({\bf r}) = \sum_{j=0,\frac{1}{2},\ldots}^{\infty}\sum_{m=-j}^{j}\sum_{m'=-j}^{j} u^j_{m,m'} U^j_{m,m'}(\theta_0,\theta,\phi).$$2

Using formal charges $$\rho({\bf r}) = \sum_{j=0,\frac{1}{2},\ldots}^{\infty}\sum_{m=-j}^{j}\sum_{m'=-j}^{j} u^j_{m,m'} U^j_{m,m'}(\theta_0,\theta,\phi).$$3 for Li and $$\rho({\bf r}) = \sum_{j=0,\frac{1}{2},\ldots}^{\infty}\sum_{m=-j}^{j}\sum_{m'=-j}^{j} u^j_{m,m'} U^j_{m,m'}(\theta_0,\theta,\phi).$$4 for N, the optimized screening factor was $$\rho({\bf r}) = \sum_{j=0,\frac{1}{2},\ldots}^{\infty}\sum_{m=-j}^{j}\sum_{m'=-j}^{j} u^j_{m,m'} U^j_{m,m'}(\theta_0,\theta,\phi).$$5, interpreted as an effective dielectric screening (Deng et al., 2019).

A separate line of development generalized the local scalar density to a local vector density for systems with vector fields such as spin: $$\rho({\bf r}) = \sum_{j=0,\frac{1}{2},\ldots}^{\infty}\sum_{m=-j}^{j}\sum_{m'=-j}^{j} u^j_{m,m'} U^j_{m,m'}(\theta_0,\theta,\phi).$$6 Spatial and vector angular momenta are coupled with Clebsch–Gordan coefficients, and the invariant descriptor is a power spectrum

$$\rho({\bf r}) = \sum_{j=0,\frac{1}{2},\ldots}^{\infty}\sum_{m=-j}^{j}\sum_{m'=-j}^{j} u^j_{m,m'} U^j_{m,m'}(\theta_0,\theta,\phi).$$7

The total energy is again decomposed into local contributions and may be fit by linear Ridge regression “as in conventional spectral neighbour analysis potentials,” or by a Gaussian approximation for more nonlinear magnetic energy surfaces (Domina et al., 2022).

6. Applications, computational realization, and limits

SNAP has been applied to a broad range of materials problems that are explicitly represented in the literature. In tantalum, the foundational model reproduced solid-phase energetics, liquid structure, melting behavior, surface and defect energies, and the screw-dislocation Peierls barrier with the correct single-hump shape seen in DFT (Thompson et al., 2014). In tungsten, a SNAP model optimized with a parallel genetic algorithm reproduced elastic constants and defect formation energies closely and could be combined with a helium pair potential to model W–He defect energetics relevant to plasma-facing materials (Wood et al., 2017).

For carbon, a quadratic SNAP trained on QMD/DFT data up to 5 TPa and 20,000 K reproduced phase-diagram features, melting curves, and the shock Hugoniot within about 3% of QMD, and was used for 1 million atom polycrystalline diamond melting simulations (Willman et al., 2022). A later review/commentary presented SNAP as an enabling technology for revisiting the carbon phase diagram, including the diamond-to-BC8 transition, metastability, and nucleation windows under dynamic compression (Correa et al., 2 Jun 2025).

For cobalt, q-SNAP closely reproduced phonon dispersions of hcp and fcc Co, surface energies, and nanoparticle stability ordering. The work emphasized its numerical stability in long MD: a 100 ns simulation of a 2406-atom truncated octahedron at 300 K remained stable throughout, the nanoparticle heat capacity converged to less than $$\rho({\bf r}) = \sum_{j=0,\frac{1}{2},\ldots}^{\infty}\sum_{m=-j}^{j}\sum_{m'=-j}^{j} u^j_{m,m'} U^j_{m,m'}(\theta_0,\theta,\phi).$$8 after 100 ns, and extrapolation of size-dependent melting temperatures yielded a bulk melting point of 1785 K, only 15 K above experiment (Bideault et al., 2024).

For niobium, a linear SNAP with ZBL short-range repulsion was trained to resolve a long-standing defect-physics failure of classical Nb potentials. It reproduced the DFT ordering of self-interstitial dumbbells, predicting the $$\rho({\bf r}) = \sum_{j=0,\frac{1}{2},\ldots}^{\infty}\sum_{m=-j}^{j}\sum_{m'=-j}^{j} u^j_{m,m'} U^j_{m,m'}(\theta_0,\theta,\phi).$$9 dumbbell as the lowest-energy SIA configuration, and in 5 keV collision cascade simulations it produced 98% $u^j_{m,m'}$0 single dumbbells in the primary-damage state, whereas FS and EAM models overwhelmingly yielded $u^j_{m,m'}$1 dumbbells (Bhardwaj et al., 5 Feb 2025).

SNAP has also been integrated into surrogate workflows rather than used only as a final force field. In ternary convex-hull construction, ensembles of binary-trained SNAPs were used to pre-relax and screen very large pools of ternary derivative structures before DFT refinement. In the Mo–Ta–W case, the workflow evaluated 331,734 ternary structures using SNAP energies and forwarded only the 15 best candidates per stoichiometry to DFT (Rossignol et al., 2023).

The formalism carries a clear computational cost. The original SNAP paper states that it is more expensive than simpler empirical potentials by roughly one to two orders of magnitude in that implementation because bispectrum and bispectrum-derivative evaluation are costly (Thompson et al., 2014). Performance engineering has therefore been a major parallel development. In the Kokkos CUDA backend of LAMMPS, algorithmic refactorization, adjoint reformulation, hierarchical parallelism, and memory-layout redesign produced a roughly 22× time-to-solution improvement on an NVIDIA Tesla V100-16GB relative to an existing implementation (Gayatri et al., 2020).

Several limitations are recurring rather than incidental. SNAP is fundamentally local and finite-cutoff, so short-range collision physics is often delegated to a ZBL reference term and long-range electrostatics require explicit augmentation in ionic systems (Wood et al., 2017, Deng et al., 2019). The reliability of a given fit depends on whether relevant local environments are represented in the training set; the ternary-screening studies show that poor feature-space coverage can make energies and forces unreliable along relaxation paths outside the sampled domain (Minotakis et al., 2023). Quadratic SNAP can improve accuracy substantially, but it requires substantially more data and exhibits a stronger overfitting risk (Wood et al., 2017). In extreme-condition carbon simulations, additional limitations were identified from the classical-potential side: no explicit electronic entropy and a hydrodynamic approximation in Hugoniot calculations (Willman et al., 2022).

Taken together, these results define SNAP not as a single immutable model but as a family of bispectrum-based local interatomic potentials whose central architectural choice is stable across applications: encode the neighborhood spectrally, enforce geometric invariances at the descriptor level, and fit energies, forces, and stresses through a regression model that remains explicit enough to support large-scale molecular dynamics.