# Compressing and forecasting atomic material simulations with descriptors

Thomas. D. Swinburne<EMAIL_ADDRESS>Aix-Marseille Université, CNRS,
CINaM UMR 7325, Campus de Luminy, 13288 Marseille, France

###### Abstract

Atomic simulations of material microstructure require significant resources to
generate, store and analyze. Here, atomic descriptor functions are proposed as
a general latent space to compress atomic microstructure, ideal for use in
large-scale simulations. Descriptors can regress a broad range of properties,
including character-dependent dislocation densities, stress states or radial
distribution functions. A vector autoregressive model can generate
trajectories over yield points, resample from new initial conditions and
forecast trajectory futures. A forecast confidence, essential for practical
application, is derived by propagating forecasts through the Mahalanobis
outlier distance, providing a powerful tool to assess coarse-grained models.
Application to nanoparticles and yielding of dislocation networks confirms low
uncertainty forecasts are accurate and resampling allows for the propagation
of smooth microstructure distributions. Yielding is associated with a collapse
in the intrinsic dimension of the descriptor manifold, which is discussed in
relation to the yield surface.

Materials evolve via complex, non-intuitive atomic mechanisms spanning a wide
range of time and length scalesWales (2003); Perez  _et al._ (2009). Atomic
simulations (MD) with empirical force fields offer exceptional insight, but
although spatial decomposition schemes give excellent (weak) parallel scaling
with system sizePlimpton (1995), serial time integration limits trajectory
duration, irrespective of available processorsPerez  _et al._ (2015). The
ubiquity yet high cost of MD means development of predictive techniques to
coarse-grain (CG) in space or time is an active research area Laio and
Parrinello (2002); Perez  _et al._ (2009); Mardt  _et al._ (2018); Van Der
Giessen  _et al._ (2020); Bonati  _et al._ (2021); Swinburne and Perez (2022).
Material microstructure requires large system sizes, necessitating efficient
and scalable CG techniques. Whilst many structural analysis tools
existHoneycutt and Andersen (1987); Kelchner  _et al._ (1998); Ackland and
Jones (2006); Lazar  _et al._ (2012, 2015); Stukowski  _et al._ (2012), none
provide generic compression of atomic data with a clear metric for similarity
or diversity, nor is it clear a priori how to select CG properties, leading to
massive storage requirements at scaleWu  _et al._ (2018); Zepeda-Ruiz  _et
al._ (2017). A further challenge is that simulations of materials are
typically non-equilibrium and exhibit, in part due to timescale limitations,
partially disordered structures with a dense kinetic spectrum and an unknown
steady state, often with external drivingSetyawan  _et al._ (2015); Priezjev
(2018); Mason  _et al._ (2020). To harness modern parallel computers there is
thus a recognized need to resample sparse simulation data and to forecast
simulation futures, both for physical insight and to maximize the information
yield of additional computational effortPerez  _et al._ (2015); Swinburne and
Perez (2020); Garmon  _et al._ (2022); Schaarschmidt  _et al._ (2022); Zhu
_et al._ (2021); Andrews  _et al._ (2022).
However, the complexities of material deformation limit the applicability of
current CG and acceleration schemes, which require identification of a clear
timescale separationLe Bris  _et al._ (2012) to allow parallel time
accumulationVoter (1998); Henkelman and Jónsson (1999); So and Voter (2000);
Perez  _et al._ (2009); Le Bris  _et al._ (2012); Perez  _et al._ (2015);
Swinburne and Perez (2020, 2022) or the design of low rank (typically 1-4)
collective variables (CV) which can be used to bias dynamics Voter (1997);
Henkelman  _et al._ (2000); Laio and Parrinello (2002); Darve  _et al._
(2008); Lelièvre  _et al._ (2010); Bonati  _et al._ (2021); Swinburne and
Marinica (2018). Despite many recent advancesRogal (2021); Bonati  _et al._
(2021); Baima  _et al._ (2022) general CVs for microstructure remain
elusiveSwinburne and Marinica (2018); Baima  _et al._ (2022), instead
requiring specialized simulation setups with only a few active mechanisms such
as nucleationBonati  _et al._ (2021) or the migration of isolated defectsRogal
(2021); Swinburne and Marinica (2018). Exploring unseen regions of
configuration space is known to be uncontrolled as low rank CVs may not remain
descriptiveBussi and Laio (2020). These issues extend to the powerful post-
mortem analysis toolsNoé and Nuske (2013); Mardt  _et al._ (2018); Huang  _et
al._ (2017, 2018); Xie  _et al._ (2019); Soltani  _et al._ (2022); Mardt  _et
al._ (2018); Xie  _et al._ (2019), which learn collective variables that obey
a discrete state Markov model in order to identify kinetically important
configurations with implied transition timescales. Whilst all-atomWang and
Gómez-Bombarelli (2019); Klein  _et al._ (2023); Han  _et al._ (2021);
Leimkuhler  _et al._ (2013) or coarse-beadedFu  _et al._ (2022) generative
models may provide a route for accelerated time-stepping, they are currently
only competitive to direct time integration for fairly small equilibrium
systems with a static or slowly varying bonding topology and so cannot be
applied to large-scale simulations of material deformation where a highly
transient, heterogeneous atomic connectivity is fundamental.
In this contribution, atomic descriptor functions Thompson  _et al._ (2015);
Shapeev (2016); Goryaeva  _et al._ (2021); Allen  _et al._ (2021); Lysogorskiy
_et al._ (2021) are proposed as an efficient, general and uncertainty-aware
coarse-graining approach, mapping atomic positions ${\bf
X}\in\mathbb{R}^{N\times 3}$ to a global vector $\bar{\bf
D}\in\mathbb{R}^{\sim 100}$, Eq. 1. The main results are that 1) Descriptors
can classify and regress a remarkable range of microstructural properties (see
figures) and permit a data-driven model extrapolation measureMahalanobis
(1936), transferring advances in active learning Podryabinkin and Shapeev
(2017); Bernstein  _et al._ (2019); Goryaeva  _et al._ (2020) to atomic CG.
This generality means CG targets need not be specified a priori, giving huge
compression in storage and efficiencies in analysis at scale. 2) Descriptor
trajectories can be efficiently resampled and forecasted via a vector
autoregressive (VAR) model Lütkepohl (2005), with, crucially, a robust
forecast uncertainty derived from the descriptor outlier measure (5). This
allows rapid assessment of when forecasts can be trusted or when additional
training is needed, essential for practical usage but typically missing in
existing schemes. The approach is applied to analyze and forecast systems
essentially untreatable with existing methods, the annealing of large
nanoparticles and dislocation yielding under cyclic shear and uniaxial
tensionZepeda-Ruiz  _et al._ (2017). Yielding is identified with a collapse in
the intrinsic dimensionGlielmo  _et al._ (2022) of the descriptor manifold.

Figure 1: Coarse graining of dislocation networks in Al under cyclic shear,
detailed in the SMSupplementary Material . Left: Global descriptor vectors
$\bar{\bf D}$ (1) are stored every 1-10ps and positions ${\bf X}$ every
100-500ps. $\\{{\bf X},\bar{\bf D}\\}$ data is used to train estimators
$\hat{\mathcal{O}}(\bar{\bf D})$ of observables $\mathcal{O}({\bf X})$ and a
VAR forecaster (8). a)-f) : $\mathcal{O}$ vs $\hat{\mathcal{O}}$ from over 20
targets, including d) dislocation junctions and total length of e) screw or f)
$\langle 112\rangle/6$ dislocations. Non-scalar a) $\sigma_{\rm xy}$,
estimated with $\bar{\bf D}\oplus\bar{\bf V}$Supplementary Material . See also
$g(r)$ in Fig. 2.

Descriptor coarse-graining of microstructure DescriptorsDeringer  _et al._
(2019); Goryaeva  _et al._ (2020); Bartók  _et al._ (2017); Batatia  _et al._
(2022) map atomic coordinates ${\bf X}\in\mathbb{R}^{N\times 3}$ to ${\bf
D}({\bf X})\in\mathbb{R}^{N\times D}$, where each element $[{\bf D}({\bf
X})]_{ij}$ takes the local atomic environment of an atom $i$ as input and
returns a permutation-invariant scalar, vector or tensor depending on the
regression target (e.g. energies, forces)Chmiela  _et al._ (2017); Batzner
_et al._ (2022)). Descriptors which approximate a many-body atomic
basisThompson  _et al._ (2015); Shapeev (2016); Goryaeva  _et al._ (2021);
Allen  _et al._ (2021); Lysogorskiy  _et al._ (2021) have found use in linear
estimators $\hat{\mathcal{O}}\simeq{\bf\Theta}_{\mathcal{O}}\cdot{\bf
D}+\Theta^{0}_{\mathcal{O}}$ of some target observable $\mathcal{O}({\bf X})$,
where ${\bf\Theta}_{\mathcal{O}}\in\mathbb{R}^{D}$ and
$\Theta^{0}_{\mathcal{O}}\in\mathbb{R}$ are parameters. For $\mathcal{O}=E$,
the atomic potential energy, these can reach state of the art accuracyBatzner
_et al._ (2022); Bartók  _et al._ (2010, 2017), often with lower computational
cost and simpler fittingGoryaeva  _et al._ (2021); Thompson  _et al._ (2015);
Allen  _et al._ (2021); Lysogorskiy  _et al._ (2021); Shapeev (2016). The
first main result of this contribution is that linear estimators can capture
essentially any microstructural property which could be of relevance to a
coarse grained model. The widely usedPlimpton (1995) bSO(4) descriptorsBartók
_et al._ (2010); Thompson  _et al._ (2015); Goryaeva  _et al._ (2021);
Plimpton (1995) are used, detailed in the supplementary material
(SM)Supplementary Material , summing over all atoms to give the global
descriptors

$\bar{\mathbf{D}}=\sum_{i}{\bf D}_{i}\in\mathbb{R}^{D},\bar{\bf
V}=\sum_{i}{\bf X}_{i}\otimes{\bm{\nabla}}_{\bf X}{\bf
D}_{i}\in\mathbb{R}^{D\times 3\times 3},$ (1)

where $\otimes$ is the outer (dyadic) product111 In practice, the tensor of
symmetric matricies $\bar{\bf V}_{s}\equiv\bar{\bf V}+\bar{\bf V}^{\top}$ is
calculated in current implementationsPlimpton (1995); Goryaeva  _et al._
(2021). Figure (1) shows linear estimators

$\hat{\mathcal{O}}(\bar{\mathbf{D}})={\bf\Theta}_{\mathcal{O}}\cdot\bar{\mathbf{D}}+\Theta^{0}_{\mathcal{O}},$
(2)

applied to dislocation networks in aluminumSupplementary Material , accurately
capturing a broad range of properties including dislocation junction
densities, character-dependent line densities and crystal structure content.
Similar results were found for the nanoparticle ensemble and a range of
dislocated solids in fcc and bcc materials. Dislocation properties were
extracted with OVITO-DXAStukowski  _et al._ (2012) which has some intrinsic
noise due to the discretization parameters. It is also possible to capture the
radial distribution function (RDF) $g(r)$ by estimating coefficients
$\hat{a}_{l}({\bf D})$ of a basis expansion $g(r)\equiv\sum_{l}a_{l}u_{l}(r)$,
as shown in Fig. 2. As found in previous work targeting vibrational
entropiesLapointe  _et al._ (2020, 2022), all predictions were stable under
widely varying test/train ratios and truncation of training data range.
Matrix-valued observables such as the stress $\mathcal{O}({\bf
X})={\bm{\sigma}}\in\mathbb{R}^{3\times 3}$ can be estimated by building
equivariant estimators with $\bar{\bf V}$; the simplest ($l=0$Batzner  _et
al._ (2022)) example is simply
$\hat{\mathcal{O}}(\bar{\mathbf{D}})={\bf\Theta}_{\mathcal{O}}\cdot\bar{\bf
V}\in\mathbb{R}^{3\times 3}$. Examples for the non-scalar shear stress
$\sigma_{xy}$ are shown in figure (1) and the SMSupplementary Material .
However, in the following only $\bar{\bf D}$ is used for forecasting,
targeting the scalar pressure ${\rm Tr}({\bm{\sigma}})$, as model parameters
are scalars and $\bar{\bf D}$ has a metric distanceThompson  _et al._ (2015).
Whilst (2) is trained on the global descriptor signal (1), the same procedure
can be applied to spatially-dependent signal from atoms in some voxel
discretization. Further investigation of this spatially dependent signal and
constraints required for any forecast will be the subject of a future
contribution.
The accuracy and scope of (2) has particular relevance for massively parallel
workflows, as only $\bar{\mathbf{D}},\bar{\bf V}$ need to be stored to later
extract almost any global observable of interest a posteriori after training
on a small database of stored positions, offering massive data compression.

Unimodality and generation of descriptor data As the descriptors have a metric
distance, similar microstructures will be close in descriptor space. In
addition, their distribution in sufficiently high dimension can be expected to
be unimodal, routinely invoked in active learning schemes Podryabinkin and
Shapeev (2017); Bernstein  _et al._ (2019) and more recently in the analysis
of defect structuresGoryaeva  _et al._ (2020). Evidence for nanoparticle and
dislocation ensembles is provided in the SMSupplementary Material . It is then
simple to generate plausible descriptor vectors by fitting and sampling a
multivariate normal distribution $\mathcal{N}({\bm{\mu}},{\bm{\Sigma}})$ to
the descriptor dataset. An example of this is shown below in figure 3, where
the observed descriptor initial conditions are densely interpolated, allowing
the evolution of observable distributions to be monitored.

Resampling and forecasting of descriptor trajectories At regular intervals
$t_{n}=n\delta\tau$, $\delta\tau\simeq 10$ps, a ‘snapshot’ is taken by time
averaging $\bar{\bf X}_{n}=\tau^{-1}_{D}\int_{0}^{\bar{\tau}}{\bf
X}(t_{n}+t){\rm d}t$ over a period $\bar{\tau}\simeq 20-50$fs to reduce noise
from thermal fluctuationsSwinburne  _et al._ (2014), then calculating
descriptor vectors $\bar{\mathbf{D}}_{n}=\bar{\mathbf{D}}(\bar{\bf X}_{n})$. A
small database of positions $\bar{\bf X}_{n}$ is built by recording $1-5\%$ of
snapshots, though positions could be selected adaptively to maximise training
diversity. An ensemble of $M$ simulations thus produces $M$ discrete time
trajectories $\\{\bar{\mathbf{D}}_{n}\\}$, which are used to train a $P$-state
vector autoregressive VAR(P) modelLütkepohl (2005); Karlsson (2013)

$\bar{\mathbf{D}}_{n+1}=\sum_{p=0}^{p=P-1}{\bf
T}_{p}\bar{\mathbf{D}}_{n-p}+{\bf c}+{\bf w}_{n},\,\langle{\bf
w}_{n}^{\top}{\bf w}_{m}\rangle={\bf S}\delta_{nm}.$ (3)

For $P>1$ a Wold transformationWold (1948) ${\bf
Z}_{n}=1\oplus\bar{\mathbf{D}}_{n}\dots\oplus\bar{\mathbf{D}}_{n-p}\in\mathbb{R}^{1+P\tilde{D}}$
casts (8) as a Markovian Ornstien-Uhlenbeck equationCoffey and Kalmykov (2012)
${\bf Z}_{n+1}={\bf T}{\bf Z}_{n}+{\bf W}_{n}$. The maximum likelihood
estimator of ${\bf T}$ is simply the least squares solution, with ${\bf S}$
determined from the residual covarianceKarlsson (2013). To minimize
generalization error a baggingBreiman (1996); Lakshminarayanan  _et al._
(2017); Petropoulos  _et al._ (2018) approach was developed, applying Bayesian
ridge regressionPedregosa  _et al._ (2011) to random overlapping subsets.
Results were stable under 10-40 subsets each with 10-40% coverage, giving
epistemic uncertainties $\delta{\bf T},{\delta{\bf S}}$ from the covariance
across subsets. Training is robust and requires only a few CPU minutes, a key
advantage over (RNN/LSTM) neural networksYu  _et al._ (2019); Vlachas  _et
al._ (2021) or neural differential equationsKidger  _et al._ (2020) which
require significant resources, regularisation/correction schemesFu  _et al._
(2022), and limited in practice to data dimension $\tilde{D}<10$Liu  _et al._
(2019); Kidger  _et al._ (2020). A Chapman-Komologorov testMardt  _et al._
(2018) for the transfer matrix ${\bf T}$ is provided in the SMSupplementary
Material , but in practice the light computational demand also permits a
convergence test of model architecture by increasing $P$Supplementary Material
.
Deriving a forecast uncertainty Practical application of (8) requires a robust
measure of forecast uncertaintyPerez  _et al._ (2015); Swinburne and Perez
(2020); Garmon  _et al._ (2022); Schaarschmidt  _et al._ (2022); Zhu  _et al._
(2021); Andrews  _et al._ (2022), which should be larger for configurations
further from the training data independent of epistemic errors. This is
particularly relevant to the non-stationary dynamics of material deformation.
As uncertainty to previously unseen macroscopic changes is clearly not
quantifiableSwinburne and Perez (2020),the following bound is conditional on
the simulation ensemble remaining unimodal and not undergoing macroscopic
changes. Many extrapolation grade estimators have been developed for active
learning of energy modelsBartók  _et al._ (2010); Podryabinkin and Shapeev
(2017); Goryaeva  _et al._ (2020); Glielmo  _et al._ (2022); Zeni  _et al._
(2022); here, the Mahalanobis outlier distanceMahalanobis (1936) is used for
the unimodal descriptor distributionSupplementary Material ; Goryaeva  _et
al._ (2020). With training data mean ${\bm{\mu}}_{\rm tr}$ and covariance
${\bm{\Sigma}}_{\rm tr}$ estimated via a shinkage estimatorChen  _et al._
(2010), the squared Mahalanobis distance reads

$\mathcal{M}(\bar{\mathbf{D}})=\left[\bar{\mathbf{D}}-{\bm{\mu}}_{\rm
tr}\right]{\bm{\Sigma}}_{\rm tr}^{-1}\left[{\bf D}-{\bm{\mu}}_{\rm
tr}\right]/{\tilde{D}}.$ (4)

Importantly, (12) is independent of the VAR(P) forecast model (8); points
drawn from a low density region of $\rho_{\rm tr}$ will have a large
Mahalanobis distance, even if epistemic uncertainties $\delta{\bf T}$ are
small. At long forecasting times, (8) will reach its high dimensional steady
stateSupplementary Material , with $\langle\mathcal{M}\rangle$ constant.
However, model parameters cannot be assumed static, with a time dependence
bounded from below by $1/\tau_{M}=1/(M\tau_{\rm tr})$Swinburne and Perez
(2018), where $M$ is the ensemble size and $\tau_{\rm tr}$ training duration.
This drift can be estimated by propagating epistemic uncertainty in the steady
state to an uncertainty $\sigma^{2}_{\mathcal{M}}$ in
$\mathcal{M}(\bar{\mathbf{D}})$, which should be accumulatedSupplementary
Material , leading to an additional linear growth in (12) of

$\mathcal{M}(t_{n})=\langle\mathcal{M}(\bar{\mathbf{D}}_{n})\rangle+\mathcal{M}_{\sigma}(t_{n}),\quad\mathcal{M}_{\sigma}(t_{n})\geq\mathcal{M}_{0}(t_{n}).$
(5)

where $\mathcal{M}_{\sigma}(t)=\sigma^{2}_{\mathcal{M}}t/\delta\tau$ and
$\mathcal{M}_{0}(t)=t/\tau_{M}$. Equation (5) is the main theoretical result
of this contribution, an uncertainty metric for forecasting via (8). An
approximate parallel efficiency is implied by $\eta=\tau_{\rm
pred}/(M\tau_{\rm tr})$, where $\mathcal{M}(\tau_{\rm
pred})\equiv\mathcal{M}_{0}(M\tau_{\rm tr})=2$, giving $\eta=1$ when
$\mathcal{M}(t_{n})=\mathcal{M}_{0}(t)$.

Figure 2:  Annealing of Pt nanoparticles. a) Representative structure at
0,11ns. b) The average RDF $g(r)$ and the corresponding descriptor
estimationSupplementary Material $g(r;\bar{\bf D})$. c)-h) Ensemble data with
$M=60$, $\tau_{\rm tr}=0.5,1.0,1.5$ns, training starting at $t=$1ns (left-
right, red shade). Mean is solid line, with standard deviation as bands.
Black: MD data. Orange: VAR forecasts from 1ns. Blue: VAR forecast with
epistemic errors. c)-e) Potential energy change from 1ns mean. f)-h)
Mahalanobis distances, MD: $\mathcal{M}({\bf D}_{t})$, eq. (12), forecasts:
$\mathcal{M}(t)$, eq. (5). The theoretical lower bounds $\mathcal{M}_{0}$ (red
dash) and $\mathcal{M}_{\sigma}$ (purple dash) are also shown.

Annealing of Pt Nanoparticles Metallic nanoparticles are important functional
materials for catalysis; 50-150 atom clusters have been extensively studied in
simulationsWales (2003); Baletto  _et al._ (2002); Huang  _et al._ (2017,
2018), but for large sizes and high temperatures the landscape of energy
minimia is vast and insufficiently metastable for current acceleration
methodsHuang  _et al._ (2018). The current application to $M=$60 4000-atom
EAM-PtLiu  _et al._ (1991) nanoparticles at 900K is thus untreatable with
existing methods.

The initial structure was formed by quenching from the liquid state and
annealing for 100ps to give a highly disordered but predominantly fcc
structure ($c_{\rm FCC}\simeq 0.5$). Descriptor trajectories were extracted
every 1.5 ps, with a full structural analysis undertaken every 100ps, though
the dataset was sparsified by taking $\delta{\tau}=15$ps and removing
intermediate snapshots. Autoregressive models (8) were constructed with
$P=1-3$ and $\tau_{\rm tr}$=0.5, 1.0 or 1.5ns, with $P=1$ shown. Generated
trajectories were launched from the start of the training stage, meaning the
observed trajectories were resampled and then forecasted. Figure 2 displays
the ensemble simulation data, model predictions and epistemic errors for the
formation energy, the RDF $g(r)$ and the Mahalanobis uncertainty (5). The RDF
reflects the significant growth in FCC crystal structure, as can also be
directly extracted through estimation of $c_{\rm FCC}$, as shown in the
SMSupplementary Material . MD data used $\mathcal{M}({\bf D})$, eq. (12),
which closely follows the theoretical lower bound $\mathcal{M}_{\sigma}(t)$.
Whilst the prediction error systematically improves as training data increases
in duration and diversity, of central importance is that this is reflected in
the magnitude of $\mathcal{M}(\bar{\mathbf{D}},t)$, confirming that
predictions are reliable.

Figure 3:  Yielding of dislocation networks under uniaxial tension in Al.
Color scheme follows figure 2. Forecasts are from 0.2ns, with $\tau_{\rm
tr}=1.2$ns and $P=5$. Observable plots a) Pressure, b) $\langle 112\rangle/6$
dislocation density and c) amorphous content. d) Mahalanobis distance. e)-h)
Trajectory resampling from 0.3ns (pre-yield), with $100\times$ larger
ensemble. i) PCA analysis of the ensemble mean $\langle{\bar{\bf D}}\rangle$,
clearly showing a localization on yield. Individual trajectories shown as
histogram in grayscale. j) ID of the descriptor manifold, estimated via the
TwoNNFacco  _et al._ (2017) and DANCoCeruti  _et al._ (2014) methods. Both
show a collapse on yielding, but remain above the Von Mises lower bound.

Yielding of Al under uniaxial tension Dislocations the agents of plastic
deformation and form dense networks under irradiationMason  _et al._ (2020) or
extended loadingLavenstein  _et al._ (2020), with atomic simulations still
producing unexpected atomic mechanisms even in extensively studied pure
materialsBertin  _et al._ (2022). Understanding yield is a central goal of
metal physicsBaggio  _et al._ (2019), but atomic analysis is challenging as
collective variables are elusiveBaima  _et al._ (2022) and the dynamics are
insufficiently metastable or stationary for the methods discussed above. An
ensemble of $M=60$ dense dislocation networks were formed in an EAM model of
AlOnat and Durukanoğlu (2013) by creating simulation boxes of around
$1.5\times 10^{5}$ atoms, orientated to $[10\bar{1}],[1{1}1],[1\bar{2}1]$,
with populations of interstitial loops. The initial dislocation densities
spanned the typical MD rangeZepeda-Ruiz  _et al._ (2017); Bertin  _et al._
(2022) of $\rho_{\rm dis}\in[10^{11},10^{13}]{\rm cm}^{-2}$. Uniaxial tension
was applied at a rate $\dot{\epsilon}_{xx}=1\times 10^{8}{\rm s}^{-1}$ along
$[10\bar{1}]$, allowing other supercell dimensions to relaxZepeda-Ruiz  _et
al._ (2017). The SMSupplementary Material  shows application to cyclic shear
loading. The results are summarized in figure 3a)-d), using the linear
estimators (2). Increasing $\tau_{\rm tr}$ decreased error and uncertainty-
convergence tests showed optimal results with $P=5$Supplementary Material .
Training only on pre-yield structures led to unstable forecasts as yield is
characterized by a qualitative change in the descriptor manifold as detailed
below. However, $\mathcal{M}(\bar{\bf D},t)$ also diverged at the yield point,
clearly indicating that more training data is required. This again
demonstrates the utility and critical importance of a forecast uncertainty to
asses data-driven predictions. Resampling allows for ensembles to be increased
by orders of magnitude for negligible CPU effort, giving smooth microstructure
distributions as shown in figure 3e)-h). Initial descriptor states were
generated as described above from $\rho_{\rm
tr}\simeq\mathcal{N}({\bm{\mu}}_{0},{\bm{\Sigma}}_{0})$, fit from the
descriptor ensemble at times 0.3-0.31ns. The forecasted ensemble captures
multiple important microstructural evolutions that, whilst known for this
well-studied systemZepeda-Ruiz  _et al._ (2017), confirm the accuracy of the
VAR approach. Forecasts correctly predict the growth of amorphous atomic
environments due to defect production under continued loadingMarian  _et al._
(2004), the expected sharp peak in HCP content at yield, accompanied by a
growth, peak then steady state in the number of dislocation junctions (see
SMSupplementary Material ). Distributions can tighten or widen, here
indicating the evolution in dislocation character- initial populations of
$\langle 100\rangle/3$ Hirth dislocation loops decay to a tight distribution
close to zero upon loading, accompanied by an emergence of a broad, stable
distribution of $\langle 112\rangle/6$ dislocation lines that carry the
plastic flowZepeda-Ruiz  _et al._ (2017). The joint stability of junctions,
dislocation density and stress is consistent with a Kocks-Mecking steady
stateSills  _et al._ (2018). The descriptor data is highly sensitive to
yielding; figure 3i) shows the first two PCA components of the ensemble mean
$\langle D\rangle$, revealing a clear localization after yielding which can be
easily classified (see SMSupplementary Material ). Theoretical models for
yielding invoke the concept of a yield surface in 5D stress spaceHirth and
Lothe (1991), which for metallic systems is typically the Von Mises yield
surface, isosurfaces of the $J_{2}$ invariant with intrinsic dimension
(ID)Ceruti  _et al._ (2014); Glielmo  _et al._ (2022) of 2. Yielding is thus
expected to be accompanied by an abrupt drop in the intrinsic dimensionality
of the stress trajectory. As stresses are essentially deterministic from
descriptors, the descriptor ID is an upper bound to the yield surface ID. Two
empirical ID estimatorsCeruti  _et al._ (2014); Facco  _et al._ (2017); Bac
_et al._ (2021) were applied to the full descriptor data $\bar{\bf
D}\oplus\bar{\bf V}$. Figure 3j) shows the estimated ID collapses from around
5-7 to around 2-3 on yield. Although these typically underestimateCeruti  _et
al._ (2014), this is consistent with a Von Mises ID lower bound of 2.
Furthermore, this indicates the existence of a yield manifold in descriptor
space, generalizing the yield surface concept to a much richer description of
microstructure than stress alone. It is speculated that the yield manifold
provides a route for data-driven construction of advanced structure-property
relationships.
Conclusions This contribution has promoted scalar and matrix-valued
descriptors as a compressed representation of atomic microstructure ideal for
analysis, resampling and forecasting of simulations, with a robust forecast
uncertainty derived using outlier distances. Analysis of the descriptor
manifold indicates the existence of a generalized yield surface, a promising
direction for future research, alongside the use of forecasting in autonomous
resource allocationPerez  _et al._ (2015); Swinburne and Perez (2020); Garmon
_et al._ (2022); Schaarschmidt  _et al._ (2022); Zhu  _et al._ (2021); Andrews
_et al._ (2022) and extension to a spatially dependent, fully equivariant
descriptor framework.

##  I Acknowledgments

I thank M-C Marinica and L Truskinovsky for stimulating discussions, an
anonymous referee for careful reading of the manuscript, the ANR grant
ANR-19-CE46-0006-1, IDRIS allocations A0090910965, A0120913455, and Euratom
grant No 633053.

##  II Supplementary material for ‘Coarse-graining and forecasting atomic
simulations of solids with descriptor functions‘

###  II.1 Choice of descriptor function basis

Descriptors which approximate a many-body atomic basisThompson  _et al._
(2015); Shapeev (2016); Goryaeva  _et al._ (2021); Allen  _et al._ (2021);
Lysogorskiy  _et al._ (2021) have found use in linear (LML) models
$\mathcal{O}={\bf\Theta}_{\mathcal{O}}\cdot\sum_{i}{\bf D}_{i}$ of some target
$\mathcal{O}$, where ${\bf\Theta}_{\mathcal{O}}\in\mathbb{R}^{D}$ is a vector
of parameters 222offsets are absorbed into the first entry of the per-atom
descriptor ${\bf D}_{i}$ (i.e. ${\rm D}_{i1}=1$). A number of possible basis
functions have been proposed; an approximate classification can be made into
‘compact‘ spectral expansions of 50-200 termsGoryaeva  _et al._ (2021);
Thompson  _et al._ (2015) and ‘non-compact‘ polynomial expansions of
1000-10000 terms, which offer greater accuracy but must be carefully
regularised to avoid overfittingAllen  _et al._ (2021); Lysogorskiy  _et al._
(2021); Shapeev (2016). In this work, the bSO(4) bispectral descriptor
functions are usedBartók  _et al._ (2010); Thompson  _et al._ (2015); Goryaeva
_et al._ (2021), as implemented in the MILADY potential package, first
introduced as part of the SNAP family of LML potentials Thompson  _et al._
(2015). Briefly, let ${\rm B}_{ji}({\bf X}),j\in[0,{\rm N_{B}}]$ be the ${\rm
N_{B}}$ bispectral components for an atom $i$, along with a constant component
${\rm B}_{0i}\equiv 1$. Only neighboring atoms within the cutoff distance
(here 4.7 Å) are included in the descriptor function calculation, which for
the atom $i$ reads

${\bf D}_{i}({\bf X})=\sum_{i}{\bigoplus}_{j}{\rm B}_{ji}({\bf
X})\in\mathbb{R}^{{\rm N_{B}}},$ (6)

where $\oplus$ indicates concatenation, giving $({\rm N_{B}}+1)$ components.
The number of bispectrum components ${\rm N_{B}}$ is determined by an angular
moment parameter $j_{max}$, here set to $4$, giving ${\rm N_{B}}=55$.

###  II.2 Descriptor estimation of correlation functions

As discussed in the main text, it is also possible to predict the value of
correlation functions such as the radial distribution function $g(r)$, which
could be written as $g(r;{\bf X})$, as it is deterministic from a given set of
positions ${\bf X}$. Our goal is to make a basis expansion

$g(r;{\bf X})=\sum_{l}a_{l}({\bf X})u_{l}(r),$ (7)

where $u_{l}(r)$ are linearly independent. As the coefficients $a_{l}({\bf
X})$ are deterministic scalar functions of $\bf X$ they are valid targets for
estimation by a descriptor estimator $\hat{a}_{l}({\bf D})$. In practice,
$g(r;{\bf X})$ is evaluated up to a multiplicative constant as a histogram
with a vector of $H$ counts ${\bf g}({\bf X})\in\mathbb{Z}^{H}_{+}$ in $H$
bins with centers ${\bf r}$. Fixing ${\bf r}$, an appropriate discrete basis
${\bf u}_{l}\in\mathbb{R}^{H}$ can be found through singular value
decomposition of the rectangular matrix of $T$ training points ${\bf G}=[{\bf
g}({\bf X}_{1}),{\bf g}({\bf X}_{2}),\dots]\in\mathbb{R}^{T\times H}$. It is
then simple to learn the coefficient estimators $\hat{a}_{l}({\bf D})$ for as
many singular vectors as required, which in practice was found to be 4 or 5.
It is noted that $g(r)$ was only considered as an observable after simulations
has been run; this is a key advantage of the current approach, as adding
$g(r)$ did not require any new data generation, only estimator training on the
1% of retained positions then prediction on the full descriptor dataset. This
ability to select coarse-graining targets a posteriori is a central advantage
of the current method.

###  II.3 Regression of matrix-valued observables

As discussed in the main text, the descriptor vectors are then summed to give
the extensive vectors $\bar{\bf D}=\sum_{i}{\bf D}_{i}$, $\bar{\bf
V}=\left[\sum_{i}{\bf X}_{i}\otimes{\bm{\nabla}}_{\bf X}{\bf D}_{i}\right]$
where $\hat{\bf e}_{\alpha},\hat{\bf e}_{\beta}\in\mathbb{R}^{3}$ are
Cartesian axes and $\otimes$ is the outer (dyadic) product. In principle, the
scalar $\bar{\bf D}$ are sufficient to regress scalar projections of matrix
quantities such as shear stresses $\sigma_{xy}={\bf
e}_{x}\cdot{\bm{\sigma}}\cdot{\bf e}_{y}$. However, whilst not used in the
main text, it is possible to form scalar projections $\bar{\bf
V}_{\alpha\beta}=\bar{\bf V}:{\bf e}_{\alpha}\otimes{\bf e}_{\beta}$. It is
then possible to regress shear stress components $\sigma_{xy}={\bf
e}_{x}\cdot{\bm{\sigma}}\cdot{\bf e}_{y}$. The effect of this can be seen in
figure 4.

Figure 4: Linear models with 50 components from PCA of $\bar{D}$ (left) and
150 components from PCA on $\bar{\bf D}\oplus\bar{\bf V}$. Data is for the Al
under uniaxial loading presented in the main text, with the same training
regime.

It is possible to form matrix-equivariant linear models using only scalar
coefficients $\hat{\mathcal{O}}={\bf\Theta}^{\top}\bar{\bf V}$, shown in
figure 5. Use of $\bar{\bf D}\oplus\bar{\bf V}$ as in figure 4 gave slightly
superior results; future work will investigate equivariant use of $\bar{\bf
D}$ e.g. $\hat{\mathcal{O}}={\bf\Theta}^{\top}\bar{\bf V}+\bar{\bf
D}^{\top}{\bm{\Phi}}\bar{\bf V}$.

Figure 5: Matrix-equivariant ($l=0$ coefficients) model for the stress
$\hat{\bf\sigma}={\bf\Theta}_{\sigma}^{\top}\bar{\bf V}$. Data is for the Al
under uniaxial loading presented in the main text.

###  II.4 $P$-convergence of VAR models

In addition to Chapman-Komologorov tests (see below) it is also simple to
study convergence in model predictions with the memory $P$. As can be seen in
figure 6 , the predictions are stable and mildly improve with increasing $P$.

Figure 6: Convergence test of VAR approach for Dislocation networks in Al
under uniaxial tension, as presented in the main text. Left-right: VAR models
with $P=1,3,5$. $P=5$ was used in the main text.

###  II.5 Chapman-Komologorov tests

VAR models with $P=1$ are equivalent to extended dynamic mode
decompositionWilliams  _et al._ (2015), for which a Chapman-Komologorov test
is possible. In general, in implied timescales
$\omega=\ln|\lambda_{\tau}|/\tau$ should be a constant. As can be seen, this
is satisfied for $\tau>10$ps.

Figure 7: Chapman-Komologorov test, $\omega=\ln|\lambda_{\tau}|/\tau$, for the
Pt nanoparticle and Al dislocation ensembles. The implied timescales
$1/\omega$ are constant for $\tau\geq$10ps, which is the employed interval.

###  II.6 Steady state of autoregressive models

The VAR(P) vector autoregressive model readsLütkepohl (2005)

$\tilde{\bf D}_{n}=\sum_{p=1}^{P}{\bf T}_{p}\tilde{\bf D}_{n-p}+{\bf
c}+\delta{\bf W}_{n},\quad\langle\delta{\bf W}_{n}^{\top}\delta{\bf
W}_{m}\rangle={\bf S}\delta_{nm}.$ (8)

Consider the vector ${\bf Z}_{n-1}=[\tilde{\bf D}_{n-1},\dots,\tilde{\bf
D}_{n-p}]$. It satisfies the extended VAR(1) model via a Wold
transformationWold (1948)

${\bf Z}_{n}=\left[\begin{matrix}{\bf T}_{1}&\dots&{\bf T}_{p}\\\
\mathbb{I}&{\bf 0}&\dots\\\ {\bf 0}&\mathbb{I}&\dots\\\ {\bf
0}&\vdots&\dots\\\ \dots&\mathbb{I}&{\bf 0}\end{matrix}\right]{\bf
Z}_{n-1}+\left[\begin{matrix}{\bf c}+\delta{\bf W}_{n}\\\ {\bf 0}\\\ {\bf
0}\\\ \vdots\\\ {\bf 0}\\\ \end{matrix}\right]={\bf T_{Z}}{\bf
Z}_{n}+\delta{\bf Y}_{n}.$ (9)

We can now solve for the steady state as for any other VAR(1), or Ornstien-
Uhlenbeck equationCoffey and Kalmykov (2012), with a Gaussian steady state

$\lim_{t\to\infty}\rho({\bf Z},t)=\mathcal{N}({\bf Z}_{\infty},{\bf
Q}_{\infty})$ (10)

where

$\left[\mathbb{I}-{\bf T_{Z}}\right]{\bf Z}_{\infty}\equiv{\bf c},\quad{\bf
S}={\bf L}^{\top}{\bf L},\,\left[\mathbb{I}-{\bf T_{Z}}\right]{\bf V}={\bf
L},\quad{\bf Q}_{\infty}={\bf V}^{\top}{\bf V},$ (11)

where we express the mean and variance in a form which allows for a null space
of $\left[\mathbb{I}-{\bf T_{Z}}\right]$ in the case of inadequate training
data or underregularisationRasmussen (2004).

###  II.7 Unimodality of descriptor distributions

As discussed in the text, the error metrics we employ rely on approximately
unimodal descriptor distrinbutions. Figure (8) shows the descriptor datasets
used for the forecasting model, the raw descriptor data reduced to 50
dimensions and whitened (decorrelated) through application of the PCA
techniqueRasmussen (2004). As can be seen, the histograms are unimodal to a
high degree of approximation. The only slight deviation came from one
component for the Pt nanoparticles, but we attribute this primarily to the
relatively small sample size of 60 simulation workers in a small atomic
system. For the Al dislocation system, the unimodality is clear.

Figure 8: Descriptor histograms for the Pt nanoparticles (top) and Al
histograms (bottom).

###  II.8 Uncertainty propagation to the Mahalanobis distance

The bagging regression proceedure detailed in the main text returns
expectation values and covariances for model parameters ${\bf T_{Z}},{\bf
c},{\bf S}$. In practice, we decompose ${\bf S}={\bf X}^{\top}{\rm
Diag})({\bm{\lambda}}){\bf X}$ and only consider covariance of the positive
eigenvalues ${\bm{\lambda}}$. Epistemic uncertainties in ${\bf T_{Z}}$ where
found to have negligable influence and thus are not considered here. These
parameter uncertainties, expressed by the covariances ${\bm{\Sigma}}_{\bf c}$
and ${\bm{\Sigma}}_{\bm{\lambda}}$ of ${\bf c}$ and ${\bm{\lambda}}$, will
induce a random perturbation on the square Mahalanobis distance

$\mathcal{M}(\bar{\mathbf{D}})=\left[\bar{\mathbf{D}}-{\bm{\mu}}_{\rm
tr}\right]{\bm{\Sigma}}_{\rm tr}^{-1}\left[{\bf D}-{\bm{\mu}}_{\rm
tr}\right]/{\tilde{D}}.$ (12)

As discussed in the main text, as this uncertainty always increases
$\mathcal{M}$ it cannot be assumed to be a mean zero fluctuation and thus
should be accumulated, leading to a linear growth in time. It is simple to
propagate this uncertainty to $\mathcal{M}(\bar{\mathbf{D}})$:

$\sigma^{2}_{\mathcal{M}}={\rm Tr}\left(\left[\mathbb{I}-{\bf
T_{Z}}\right]^{-1}{\bm{\Sigma}}_{\bf c}{\bm{\Sigma}}_{\rm
tr}^{-1}\right)/{\tilde{D}}+{\rm Tr}\left(\left[\mathbb{I}-{\bf
T_{Z}}\right]^{-1}{\bf X}^{\top}{\bm{\Sigma}}_{\bm{\lambda}}{\bf
X}{\bm{\Sigma}}_{\rm tr}^{-1}\right)/{\tilde{D}}$ (13)

In practice, we ensure that $\sigma^{2}_{\mathcal{M}}>\delta\tau/\tau_{M}$
though the mapping $\sigma^{2}_{\mathcal{M}}\to
1/(1/\sigma^{2}_{\mathcal{M}}+\tau_{M}/\delta\tau)$, though as can be seen in
the figures the addition of this has only a very minor influence on the
expected uncertainty.

###  II.9 Yielding of dislocation networks under cyclic shear

An ensemble of $M=60$ dense dislocation networks were formed in an EAM model
of AlOnat and Durukanoğlu (2013) by creating simulation boxes of around
$1.5\times 10^{5}$ atoms, orientated to $[10\bar{1}],[1{1}1],[1\bar{2}1]$,
with populations of interstitial loops. The initial dislocation densities
spanned the range typically found in MDZepeda-Ruiz  _et al._ (2017); Bertin
_et al._ (2022), $\rho_{\rm dis}\in[10^{11},10^{13}]{\rm cm}^{-2}$. Ensembles
were subjected to cyclic shearing at a rate $|\dot{\epsilon}_{xy}|=1\times
10^{8}{\rm s}^{-1}$ up to $|\epsilon_{xy}|=0.1$ with a sawtooth profile for a
total of two cycles. As can be seen, the initially disperse dislocation
densities across the ensemble tightens after a few cycles, as the driven
steady state emerges.

Figure 9: Cyclic shearing of dislocation networks in Al. Left: Representative
simulation supercell. a)-c): Forecasts using 2.5ns (top) and 3ns (bottom) for
training and resampling, showing a) shear stress, b) total dislocation density
and c) outlier distance. For lower training data the outlier distance is
clearly seen to spike at yeild points, these being the least sampled
configurations. d)-e): resampling, increasing the effective ensemble size one
hundredfold. The total dislocation density distribution d) noticably tightens.
demonstrating a convergence to a driven steady state, whilst the vacancy
distribution e) slightly widens, most likely due to defect generation at high
strain rates.

###  II.10 Classification of yield

A simple linear support vector machine (SVM)Rasmussen (2004) classifier was
trained on 25$\%$ of trajectories around the yield point. As can be seen in
figure 10, the transferability to new trajectories shows excellent accuracy,
both in pointwise comparison (left) and distributionwise (right).

Figure 10: SVM classification of yielding. Left: scatter plot of time at yield
vs predicted time at 50$\%$ probability. Right: observed and predicted
proportion of yielded trajectories.

## References

  * Wales (2003) D. J. Wales,  _Energy Landscapes_ , edited by C. U. Press (Cambridge, 2003).
  * Perez  _et al._ (2009) D. Perez, B. P. Uberuaga, Y. Shim, J. G. Amar, and A. F. Voter, Annual Reports in computational chemistry 5, 79 (2009).
  * Plimpton (1995) S. Plimpton, Journal Computational Physics 117, 1 (1995).
  * Perez  _et al._ (2015) D. Perez, E. D. Cubuk, A. Waterland, E. Kaxiras, and A. F. Voter, Journal of chemical theory and computation 12, 18 (2015).
  * Laio and Parrinello (2002) A. Laio and M. Parrinello, Proceedings of the National Academy of Sciences 99, 12562 (2002).
  * Mardt  _et al._ (2018) A. Mardt, L. Pasquali, H. Wu, and F. Noé, Nature communications 9, 1 (2018).
  * Van Der Giessen  _et al._ (2020) E. Van Der Giessen, P. A. Schultz, N. Bertin, V. V. Bulatov, W. Cai, G. Csányi, S. M. Foiles, M. G. Geers, C. González, M. Hütter, _et al._ , Modelling and Simulation in Materials Science and Engineering 28, 043001 (2020).
  * Bonati  _et al._ (2021) L. Bonati, G. Piccini, and M. Parrinello, Proceedings of the National Academy of Sciences 118, e2113533118 (2021).
  * Swinburne and Perez (2022) T. D. Swinburne and D. Perez, Modelling and Simulation in Materials Science and Engineering 30, 034004 (2022).
  * Honeycutt and Andersen (1987) J. D. Honeycutt and H. C. Andersen, Journal of Physical Chemistry 91, 4950 (1987).
  * Kelchner  _et al._ (1998) C. L. Kelchner, S. J. Plimpton, and J. C. Hamilton, Phys. Rev. B 58, 11085 (1998).
  * Ackland and Jones (2006) G. J. Ackland and A. P. Jones, Phys. Rev. B 73, 054104 (2006).
  * Lazar  _et al._ (2012) E. A. Lazar, J. K. Mason, R. D. MacPherson, and D. J. Srolovitz, Phys. Rev. Lett. 109, 095505 (2012).
  * Lazar  _et al._ (2015) E. A. Lazar, J. Han, and D. J. Srolovitz, Proceedings of the National Academy of Sciences 112, E5769 (2015).
  * Stukowski  _et al._ (2012) A. Stukowski, V. V. Bulatov, and A. Arsenlis, Modelling and Simulation in Materials Science and Engineering 20, 085007 (2012).
  * Wu  _et al._ (2018) G. Wu, H. Song, and D. Lin, Computational Materials Science 144, 322 (2018).
  * Zepeda-Ruiz  _et al._ (2017) L. A. Zepeda-Ruiz, A. Stukowski, T. Oppelstrup, and V. V. Bulatov, Nature 550, 492 (2017).
  * Setyawan  _et al._ (2015) W. Setyawan, G. Nandipati, K. J. Roche, H. L. Heinisch, B. D. Wirth, and R. J. Kurtz, Journal of Nuclear Materials 462, 329 (2015).
  * Priezjev (2018) N. V. Priezjev, Journal of Non-Crystalline Solids 479, 42 (2018).
  * Mason  _et al._ (2020) D. R. Mason, S. Das, P. M. Derlet, S. L. Dudarev, A. J. London, H. Yu, N. W. Phillips, D. Yang, K. Mizohata, R. Xu, _et al._ , Physical Review Letters 125, 225503 (2020).
  * Swinburne and Perez (2020) T. Swinburne and D. Perez, npj Computational Materials 6, 190 (2020).
  * Garmon  _et al._ (2022) A. Garmon, V. Ramakrishnaiah, and D. Perez, Parallel Computing , 102936 (2022).
  * Schaarschmidt  _et al._ (2022) J. Schaarschmidt, J. Yuan, T. Strunk, I. Kondov, S. P. Huber, G. Pizzi, L. Kahle, F. T. Bölle, I. E. Castelli, T. Vegge, _et al._ , Advanced Energy Materials 12, 2102638 (2022).
  * Zhu  _et al._ (2021) L.-F. Zhu, J. Janssen, S. Ishibashi, F. Körmann, B. Grabowski, and J. Neugebauer, Computational Materials Science 187, 110065 (2021).
  * Andrews  _et al._ (2022) J. Andrews, O. Gkountouna, and E. Blaisten-Barojas, Chemical Science 13, 7021 (2022).
  * Le Bris  _et al._ (2012) C. Le Bris, T. Lelievre, M. Luskin, and D. Perez, Monte Carlo Methods and Applications 18, 119 (2012).
  * Voter (1998) A. F. Voter, Physical Review B 57, R13985 (1998).
  * Henkelman and Jónsson (1999) G. Henkelman and H. Jónsson, The Journal of chemical physics 111, 7010 (1999).
  * So and Voter (2000) M. So and A. Voter, The Journal of Chemical Physics 112, 9599 (2000).
  * Voter (1997) A. F. Voter, Physical Review Letters 78, 3908 (1997).
  * Henkelman  _et al._ (2000) G. Henkelman, B. P. Uberuaga, and H. Jonsson, The Journal of Chemical Physics 113, 9901 (2000).
  * Darve  _et al._ (2008) E. Darve, D. Rodríguez-Gómez, and A. Pohorille, The Journal of chemical physics 128, 144120 (2008).
  * Lelièvre  _et al._ (2010) T. Lelièvre, G. Stoltz, and M. Rousset,  _Free energy computations: A mathematical perspective_ (World Scientific, 2010).
  * Swinburne and Marinica (2018) T. D. Swinburne and M.-C. Marinica, Phys. Rev. Lett. 120, 135503 (2018).
  * Rogal (2021) J. Rogal, The European Physical Journal B 94, 1 (2021).
  * Baima  _et al._ (2022) J. Baima, A. S. Goryaeva, T. Swinburne, J.-B. Maillet, M. Nastar, and M. C. Marinica, Phys. Chem. Chem. Phys. , (2022).
  * Bussi and Laio (2020) G. Bussi and A. Laio, Nature Reviews Physics 2, 200 (2020).
  * Noé and Nuske (2013) F. Noé and F. Nuske, Multiscale Modeling & Simulation 11, 635 (2013).
  * Huang  _et al._ (2017) R. Huang, L.-T. Lo, Y. Wen, A. Voter, and D. Perez, The Journal of chemical physics 147, 152717 (2017).
  * Huang  _et al._ (2018) R. Huang, Y. Wen, A. Voter, and D. Perez, Physical Review Materials 2, 126002 (2018).
  * Xie  _et al._ (2019) T. Xie, A. France-Lanord, Y. Wang, Y. Shao-Horn, and J. C. Grossman, Nature communications 10, 1 (2019).
  * Soltani  _et al._ (2022) S. Soltani, C. W. Sinclair, and J. Rottler, Phys. Rev. E 106, 025308 (2022).
  * Wang and Gómez-Bombarelli (2019) W. Wang and R. Gómez-Bombarelli, npj Computational Materials 5, 125 (2019).
  * Klein  _et al._ (2023) L. Klein, A. Y. Foong, T. E. Fjelde, B. Mlodozeniec, M. Brockschmidt, S. Nowozin, F. Noé, and R. Tomioka, arXiv preprint arXiv:2302.01170 (2023).
  * Han  _et al._ (2021) C. Han, P. Zhang, D. Bluestein, G. Cong, and Y. Deng, Journal of Computational Physics 427, 110053 (2021).
  * Leimkuhler  _et al._ (2013) B. Leimkuhler, D. T. Margul, and M. E. Tuckerman, Molecular Physics 111, 3579 (2013).
  * Fu  _et al._ (2022) X. Fu, T. Xie, N. J. Rebello, B. D. Olsen, and T. Jaakkola, arXiv preprint arXiv:2204.10348 (2022).
  * Thompson  _et al._ (2015) A. P. Thompson, L. P. Swiler, C. R. Trott, S. M. Foiles, and G. J. Tucker, J. Comp. Phys. 285, 316 (2015).
  * Shapeev (2016) A. Shapeev, Multiscale Model. Sim. 14, 1153 (2016).
  * Goryaeva  _et al._ (2021) A. M. Goryaeva, J. Dérès, C. Lapointe, P. Grigorev, T. D. Swinburne, J. R. Kermode, L. Ventelon, J. Baima, and M.-C. Marinica, Phys. Rev. Materials 5, 103803 (2021).
  * Allen  _et al._ (2021) A. E. Allen, G. Dusson, C. Ortner, and G. Csányi, Machine Learning: Science and Technology 2, 025017 (2021).
  * Lysogorskiy  _et al._ (2021) Y. Lysogorskiy, C. van der Oord, A. Bochkarev, S. Menon, M. Rinaldi, T. Hammerschmidt, M. Mrovec, A. Thompson, G. Csányi, C. Ortner, _et al._ , npj Computational Materials 7, 1 (2021).
  * Mahalanobis (1936) P. C. Mahalanobis (National Institute of Science of India, 1936).
  * Podryabinkin and Shapeev (2017) E. V. Podryabinkin and A. V. Shapeev, Computational Materials Science 140, 171 (2017).
  * Bernstein  _et al._ (2019) N. Bernstein, G. Csányi, and V. L. Deringer, npj Computational Materials 5, 1 (2019).
  * Goryaeva  _et al._ (2020) A. M. Goryaeva, C. Lapointe, C. Dai, J. Dérès, J.-B. Maillet, and M.-C. Marinica, Nat. Commun. 11, 4691 (2020).
  * Lütkepohl (2005) H. Lütkepohl,  _New introduction to multiple time series analysis_ (Springer Science & Business Media, 2005).
  * Glielmo  _et al._ (2022) A. Glielmo, I. Macocco, D. Doimo, M. Carli, C. Zeni, R. Wild, M. d’Errico, A. Rodriguez, and A. Laio, Patterns 3, 100589 (2022).
  * (59) Supplementary Material, .
  * Deringer  _et al._ (2019) V. L. Deringer, M. A. Caro, and G. Csányi, Advanced Materials 31, 1902765 (2019).
  * Bartók  _et al._ (2017) A. P. Bartók, S. De, C. Poelking, N. Bernstein, J. R. Kermode, G. Csányi, and M. Ceriotti, Sci. Adv. 3, e1701816 (2017).
  * Batatia  _et al._ (2022) I. Batatia, D. P. Kovács, G. N. Simm, C. Ortner, and G. Csányi, arXiv preprint arXiv:2206.07697 (2022).
  * Chmiela  _et al._ (2017) S. Chmiela, A. Tkatchenko, H. E. Sauceda, I. Poltavsky, K. T. Schütt, and K.-R. Müller, Science advances 3, e1603015 (2017).
  * Batzner  _et al._ (2022) S. Batzner, A. Musaelian, L. Sun, M. Geiger, J. P. Mailoa, M. Kornbluth, N. Molinari, T. E. Smidt, and B. Kozinsky, Nature communications 13, 1 (2022).
  * Bartók  _et al._ (2010) A. P. Bartók, M. C. Payne, R. Kondor, and G. Csányi, Phys. Rev. Lett. 104, 136403 (2010).
  * Note (1) In practice, the tensor of symmetric matricies $\bar{\bf V}_{s}\equiv\bar{\bf V}+\bar{\bf V}^{\top}$ is calculated in current implementationsPlimpton (1995); Goryaeva  _et al._ (2021).
  * Lapointe  _et al._ (2020) C. Lapointe, T. D. Swinburne, L. Thiry, S. Mallat, L. Proville, C. S. Becquart, and M.-C. Marinica, Physical Review Materials 4, 063802 (2020).
  * Lapointe  _et al._ (2022) C. Lapointe, T. D. Swinburne, L. Proville, C. S. Becquart, N. Mousseau, and M.-C. Marinica, Under Review (2022).
  * Swinburne  _et al._ (2014) T. D. Swinburne, S. L. Dudarev, and A. P. Sutton, Physical Review Letters 113, 215501 (2014).
  * Karlsson (2013) S. Karlsson, Handbook of economic forecasting 2, 791 (2013).
  * Wold (1948) H. O. Wold, The Annals of Mathematical Statistics 19, 558 (1948).
  * Coffey and Kalmykov (2012) W. Coffey and Y. P. Kalmykov,  _The Langevin equation: with applications to stochastic problems in physics, chemistry and electrical engineering_ , Vol. 27 (World Scientific, 2012).
  * Breiman (1996) L. Breiman, Machine learning 24, 123 (1996).
  * Lakshminarayanan  _et al._ (2017) B. Lakshminarayanan, A. Pritzel, and C. Blundell, Advances in neural information processing systems 30 (2017).
  * Petropoulos  _et al._ (2018) F. Petropoulos, R. J. Hyndman, and C. Bergmeir, European Journal of Operational Research 268, 545 (2018).
  * Pedregosa  _et al._ (2011) F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay, Journal of Machine Learning Research 12, 2825 (2011).
  * Yu  _et al._ (2019) Y. Yu, X. Si, C. Hu, and J. Zhang, Neural computation 31, 1235 (2019).
  * Vlachas  _et al._ (2021) P. R. Vlachas, J. Zavadlav, M. Praprotnik, and P. Koumoutsakos, Journal of Chemical Theory and Computation 18, 538 (2021).
  * Kidger  _et al._ (2020) P. Kidger, J. Morrill, J. Foster, and T. Lyons, Advances in Neural Information Processing Systems 33, 6696 (2020).
  * Liu  _et al._ (2019) F. Liu, M. Cai, L. Wang, and Y. Lu, IEEE Access 7, 26102 (2019).
  * Zeni  _et al._ (2022) C. Zeni, A. Anelli, A. Glielmo, and K. Rossi, Physical Review B 105, 165141 (2022).
  * Chen  _et al._ (2010) Y. Chen, A. Wiesel, Y. C. Eldar, and A. O. Hero, IEEE transactions on signal processing 58, 5016 (2010).
  * Swinburne and Perez (2018) T. D. Swinburne and D. Perez, Phys. Rev. Materials 2, 053802 (2018).
  * Baletto  _et al._ (2002) F. Baletto, R. Ferrando, A. Fortunelli, F. Montalenti, and C. Mottet, The Journal of chemical physics 116, 3856 (2002).
  * Liu  _et al._ (1991) C. Liu, J. Cohen, J. Adams, and A. Voter, Surface science 253, 334 (1991).
  * Facco  _et al._ (2017) E. Facco, M. d’Errico, A. Rodriguez, and A. Laio, Scientific reports 7, 12140 (2017).
  * Ceruti  _et al._ (2014) C. Ceruti, S. Bassis, A. Rozza, G. Lombardi, E. Casiraghi, and P. Campadelli, Pattern recognition 47, 2569 (2014).
  * Lavenstein  _et al._ (2020) S. Lavenstein, Y. Gu, D. Madisetti, and J. A. El-Awady, Science 370, eabb2690 (2020).
  * Bertin  _et al._ (2022) N. Bertin, W. Cai, S. Aubry, A. Arsenlis, and V. V. Bulatov, arXiv preprint arXiv:2210.14343 (2022).
  * Baggio  _et al._ (2019) R. Baggio, E. Arbib, P. Biscari, S. Conti, L. Truskinovsky, G. Zanzotto, and O. Salman, Physical Review Letters 123, 205501 (2019).
  * Onat and Durukanoğlu (2013) B. Onat and S. Durukanoğlu, Journal of Physics: Condensed Matter 26, 035404 (2013).
  * Marian  _et al._ (2004) J. Marian, W. Cai, and V. V. Bulatov, Nature materials 3, 158 (2004).
  * Sills  _et al._ (2018) R. B. Sills, N. Bertin, A. Aghaei, and W. Cai, Physical review letters 121, 085501 (2018).
  * Hirth and Lothe (1991) J. P. Hirth and J. Lothe,  _Theory Of Dislocations_ (Malabar, FL Krieger, 1991).
  * Bac  _et al._ (2021) J. Bac, E. M. Mirkes, A. N. Gorban, I. Tyukin, and A. Zinovyev, Entropy 23, 1368 (2021).
  * Note (2) Offsets are absorbed into the first entry of the per-atom descriptor ${\bf D}_{i}$ (i.e. ${\rm D}_{i1}=1$).
  * Williams  _et al._ (2015) M. O. Williams, I. G. Kevrekidis, and C. W. Rowley, Journal of Nonlinear Science 25, 1307 (2015).
  * Rasmussen (2004) C. E. Rasmussen,  _Gaussian Processes in Machine Learning_ (Springer, Berlin, Heidelberg, 2004).