# Evaluating approaches for on-the-fly machine learning interatomic potentials
for activated mechanisms sampling with the activation-relaxation technique
nouveau

Eugène Sanscartier<EMAIL_ADDRESS>Département de physique
and Regroupement québécois sur les matériaux de pointe, Université de
Montréal, Case Postale 6128, Succursale Centre-ville, Montréal, Québec H3C
3J7, Canada.  Félix Saint-Denis  Département de physique and Regroupement
québécois sur les matériaux de pointe, Université de Montréal, Case Postale
6128, Succursale Centre-ville, Montréal, Québec H3C 3J7, Canada.  Karl-Étienne
Bolduc  Département de physique and Regroupement québécois sur les matériaux
de pointe, Université de Montréal, Case Postale 6128, Succursale Centre-ville,
Montréal, Québec H3C 3J7, Canada.  Normand Mousseau
<EMAIL_ADDRESS>https://normandmousseau.com Département de
physique and Regroupement québécois sur les matériaux de pointe, Université de
Montréal, Case Postale 6128, Succursale Centre-ville, Montréal, Québec H3C
3J7, Canada.

###### Abstract

In the last few years, much efforts have gone into developing general machine-
learning potentials able to describe interactions for a wide range of
structures and phases. Yet, as attention turns to more complex materials
including alloys, disordered and heterogeneous systems, the challenge of
providing reliable description for all possible environments becomes ever more
costly. In this work, we evaluate the benefits of using specific versus
general potentials for the study of activated mechanisms in solid-state
materials. More specifically, we test three machine-learning fitting
approaches using the moment-tensor potential to reproduce a reference
potential when exploring the energy landscape around a vacancy in Stillinger-
Weber silicon crystal and silicon-germanium zincblende structure using the
activation-relaxation technique nouveau (ARTn). We find that a targeted on-
the-fly approach specific and integrated to ARTn generates the highest
precision on the energetic and geometry of activated barriers, while remaining
cost-effective. This approach expands the type of problems that can be
addressed with high-accuracy ML potentials.

††preprint: AIP/123-QED

##  I Introduction

As computational materials scientists turn to attention to ever more complex
systems, they are faced with two major challenges : (i) how to describe
correctly their physics and (ii) how to reach the appropriate size and time
scale to capture the properties of interest. The first challenge is generally
solved by turning to ab initio methods, Kohn and Sham (1965) that allow the
solution Schrödinger’s equation with reasonably controlled approximations.
Theses approaches, however, suffer from $N^{4}$ scaling which limits their
application to small system sizes and short time scales. The second challenge
is met by a variety of methods that cover different scales. Molecular dynamics
Lindahl (2008), for example, which directly solves Newton’s equation, accesses
typical time scales between picoseconds and microseconds, at the very best.
Other approaches, such as lattice Voter and Doll (1984); Voter (2007) and off-
lattice kinetic Monte-Carlo Henkelman and Jónsson (2001); El-Mellouhi,
Mousseau, and Lewis (2008), by focusing on physically relevant mechanisms, can
extend this time scale to seconds and more, as long the diffusion takes place
through activated processes. Even though these methods are efficient, each
trajectory can require hundreds of thousands to millions of forces
evaluations, becoming too costly with ab initio approaches, forcing modellers
to use empirical potentials in spite of their incapacity at describing
correctly complex environments.

Building on ab initio energy and forces, machine-learned potentials (MLP)
Behler and Parrinello (2007); Bartók, Kondor, and Csányi (2013); Thompson  _et
al._ (2015); Shapeev (2016) open the door to lifting some of this
difficulties, by offering much more reliable physics as a small fraction of
the cost of ab initio evaluations.

Since their introduction, ML potentials have been largely coupled with MD,
focusing on the search for universal potentials able to describe a full range
of structures and phases for a given material Sivaraman  _et al._ (2021);
Kang, Shang, and Liu (2020); Sivaraman  _et al._ (2020). As we turn to more
complex systems such as alloys and disordered and heterogeneous systems, it
becomes more and more difficult to generate such universal potentials, since
the number of possible environments grows rapidly with this complexity. In
this context, the development of specific potentials, with on-the-fly learning
that makes it possible to adapt to new environments, becomes a strategy worth
exploring.

In this work, we focus on the construction of machine-learned potentials
adapted to the sampling of energy landscape dominated by activated mechanisms,
i.e., solid-state systems with local activated diffusion and evolution. This
kinetics is associated with aging and relaxation of disordered materials
Béland  _et al._ (2013), sluggish diffusion in concentrated alloys Osetsky
_et al._ (2018) and defect diffusion in complex environments Restrepo  _et
al._ (2018). A correct computational sampling, using open-ended methods such
as the activation-relaxation technique (ART) Barkema and Mousseau (1996) and
its revised version (ART nouveau or ARTn) Malek and Mousseau (2000a); Jay  _et
al._ (2022a), requires a precise description of local minima and of the
landscape surrounding the first-order saddle points that characterize
diffusion according to the transition-state theory (TST) Truhlar, Garrett, and
Klippenstein (1996). These barriers can be high — reaching many electron-volts
— and involve strained configurations that can be visited only very rarely
with standard molecular dynamics. Yet, because of this need for long-time
evolution, these problems have been largely out of reach of ab initio
approaches and have been studied, until now, mostly with empirical potentials.

More specifically, we compare three machine learning approaches to train a
Moment Tensor Potential (MTP) Shapeev (2016); Novikov  _et al._ (2020) for the
diffusion of a vacancy in silicon and silicon-germanium alloy as sampled with
ARTn. The first approach consists in a pure MD learning, fitted at various
temperatures, following steps that echo the work of Novoselov et al. Novoselov
_et al._ (2019); the second approach adds an on-the-fly training during ARTn
runs and the third one focuses on a purely on-the-fly training during ARTn
runs. While on-the-fly learning is not a new approach and have been use before
in various context, they have never been compared for the specific application
of activated kinetics in solid state systems.

To generate the statistics necessary to offer solid conclusions, we select to
use the Stilliger-Weber empirical potential both for Si Stillinger and Weber
(1985) and SiGe Ethier and Lewis (1992). Clearly, the machine-learned
potentials generated here have therefore limited physical relevance, these are
sufficient to allow us to assess the relative quality of our three approaches.

Results underline the efficiency gain in developing targeted ML potentials for
specific applications, comparing the cost of fitting Si with SiGe, it also
shows the rapid increase in computation complexity associated with moving from
element to alloy systems, which emphasizes the usefulness of a specific
approach such as the one applied here to activated processes.

##  II Methodology

###  II.1 ML Potential

The Moment Tensor Potential (MTP) Shapeev (2016); Novikov  _et al._ (2020) is
a linear model of functions $B_{\alpha}(\mathfrak{r}_{i})$ built from
contractions of moment tensor descriptors defined by the local neighborhood
relative position $\mathfrak{r}_{i}$ of atom $i$ within a sphere of influence
of radius $r_{c}$ respecting a set invariances. This model has been shown to
be fast while giving accuracy on the order of $\sim$meV/atom and requiring few
hundreds to thousands of reference potential calls Zuo  _et al._ (2020) on-
the-fly.

MTP have been used on a wide variety of problems including on-the-fly MD
simulation Novoselov  _et al._ (2019); Novikov  _et al._ (2020); Podryabinkin
and Shapeev (2017), search and minimization of new alloys Podryabinkin  _et
al._ (2019); Gubaev  _et al._ (2019) and diffusion processes Novoselov  _et
al._ (2019) on systems counting one or multiple species. In the following we
offer a summary of the method; more information on MTP is available in Ref.
Novikov  _et al._ , 2020.

MTP approximates atomic configuration energy as sum of local contributions. A
local contribution is obtained through a sum over the included basis
$\\{B_{\alpha}(\mathfrak{r}_{i})\\}$ as a linear combination of
$B(\mathfrak{r}_{i})$ and $\xi_{\alpha}$,

$V(\mathfrak{r}_{i})=\sum_{\alpha=1}^{m}\xi_{\alpha}B_{\alpha}(\mathfrak{r}_{i}).$
(1)

The “level” of a potential gives the number of different possible tensor
$M_{\mu,\nu}\left(\mathfrak{r}_{i}\right)$ descriptors. The
$B_{\alpha}(\mathfrak{r}_{i})$ functions of Eq. 1 are constructed by a
tensorial contraction of different $M_{\mu,\nu}\left(\mathfrak{r}_{i}\right)$
and the number of different tensorial contraction sets $m$ in Eq. 1. More
information on MTP is available in Ref. Novikov  _et al._ , 2020.

The total energy of a N-atom configuration ($\mathfrak{R}$) is then given by
the sum of N local contributions,

$E(\mathfrak{R})=\sum_{i=1}^{N}V(\mathfrak{r}_{i})=\sum_{i=1}^{N}\sum_{\alpha=1}^{m}\xi_{\alpha}B_{\alpha}(\mathfrak{r}_{i}),$
(2)

and the forces are obtained by taking the gradient of this quantity,

$\textbf{F}(\mathfrak{R})=-\nabla\sum_{i=1}^{N}\sum_{\alpha=1}^{m}\xi_{\alpha}B_{\alpha}(\mathfrak{r}_{i}).$
(3)

The parameters $\xi_{\alpha}$ are obtained by minimizing the loss function:

$\sum_{\mathfrak{R}\in\mathrm{A}}\left[w_{e}\left(E(\mathfrak{R})-\hat{E}(\mathfrak{R})\right)^{2}+w_{f}\sum_{i}^{N}\left|\textbf{f}_{i}(\mathfrak{R})-\hat{\textbf{f}}_{i}(\mathfrak{R})\right|^{2}\right]\rightarrow\min_{\xi}.$
(4)

Here $\mathrm{A}$ is the training set made of configurations with known energy
and forces. The goal is to minimize the difference between $E(\mathfrak{R})$,
$\textbf{f}_{i}(\mathfrak{R})$(real value) and $\hat{E}(\mathfrak{R})$,
$\hat{\textbf{f}}_{i}(\mathfrak{R})$(predicted by model), respectively, for
all element in $\mathrm{A}$. Weights on contribution from energy and forces
($w_{e}$ and $w_{f}$) are set to one.

###  II.2 Learning On-The-Fly Tools

On-the-fly atomic machine learning potential (OTF) involves the repeated
training of the model potential as new atomic environments are generated
through various procedures.

Following the work of Shapeev and collaborators Novikov  _et al._ (2020), the
reliability of the potential to evaluate a given configuration is done using
the D-optimality criterion which state that the best training set is the one
with maximal volume. The volume can be seen as a measure of the domain of the
training set where the model allows reliable predictions. Assuming that the
training set is of maximal volume, any configuration within or beyond this
volume is inside or outside the known domain of the model, respectively, which
we interpret as interpolation and extrapolation. We grade with the
D-optimality criterion of a given configuration by assessing the model
reliability. For this, Shapaeev et al. introduce a selection algorithm
(MaxVol) that tests whether or not this configuration should be added to the
training set or replace a configuration already in it While a detailed
description can be found in Ref. Podryabinkin and Shapeev (2017), we provide
here a brief summary of the retained approach.

The selection and extrapolation-grade algorithm can be applied using either a
local-energy or a global-energy descriptor.

The local-energy descriptor is presented as a rectangular matrix
$\mathrm{G}_{N\times\mathrm{m}}$ formed by the basis elements
$B_{\alpha}(\mathfrak{r}_{i})$ associated with the neighborhood
$\mathfrak{r}_{i}$ of all $N$ atoms:

$\mathrm{G}=\left(\begin{array}[]{ccc}B_{1}(\mathfrak{r}_{1})&\ldots&B_{m}(\mathfrak{r}_{1})\\\
\vdots&\ddots&\vdots\\\
B_{1}(\mathfrak{r}_{N})&\ldots&B_{m}(\mathfrak{r}_{N})\end{array}\right).$ (5)

For a given configuration, the global-energy description reduces this
information to a vector $\mathbf{g}$,

$\mathbf{g}=\left(\begin{array}[]{ccc}b_{1}(\mathfrak{R})&\ldots&b_{m}(\mathfrak{R})\end{array}\right),$
(6)

where each term, $\\{b_{\alpha}(\mathfrak{R})\\}$ is a sum over all
neighborhoods for a specific basis element
$\\{B_{\alpha}(\mathfrak{r}_{i})\\}$:

$\\{b_{\alpha}(\mathfrak{R})\\}=\sum_{i=0}^{N}\\{B_{\alpha}(\mathfrak{r}_{i})\\}.$

For the global-energy descriptor of a given configuration with a training set
$\mathrm{A}$, solving for $c_{ij}$, in

$\left(\begin{array}[]{ccc}c_{11}&\ldots&c_{1m}\end{array}\right)\mathrm{A}=\mathbf{g}.$
(7)

Assuming that $\mathrm{A}$ is of maximal volume, $0\leq c_{ij}\leq 1$.
Conversely, if one find $c_{ij}>1$ then $\mathrm{A}$ is not of maximal volume.
The later case mean that the model is extrapolating on this configuration and
that updating $\mathrm{A}$ with this configuration will increase its volume.
The extrapolation grade, $\gamma$, is then defined as the largest component of
$c_{ij}$,

$\gamma(\mathfrak{R})=\max\left|c_{ij}\right|.$ (8)

The same approach is used for the local-energy description, applying Eq. 7
with the rows of matrix $\mathrm{G}$ rather than the vector $\mathbf{g}$. For
non-linear and other ML potentials, in Eq. 5, one adapt this algorithm
performing the ML model’s gradient with respect to its parameters.

In practice, for $\gamma(\mathfrak{R})$ below a certain threshold
$\gamma_{0}$, the model interpolation is reliable, while for
$\gamma_{0}<\gamma(\mathfrak{R})<\gamma_{max}$, the model cannot be applied
with confidence, but can be adapted by adding this configuration to the
training set. When $\gamma(\mathfrak{R})>\gamma_{max}$, the configuration is
too far from the training set and it is rejected as the model cannot be
adapted with confidence. In this work, we set $\gamma_{0}=1.1$ and
$\gamma_{max}=2.2$, taking the lower values of the parameters studied in Ref.
Podryabinkin and Shapeev (2017), unless specified otherwise.

###  II.3 On-The-Fly Learning Cycle Workflow

Our workflow is similar to that of Ref. Novikov  _et al._ , 2020, with main
differences discussed in Section II.6. We follow the same general machine-
learning on-the-fly workflow for all sampling approaches tested here.

We split each simulation in one or multiple sequences of atomic configurations
generated using either MD or ARTn. Each run unfolds as follows (see Fig. 1):

  1. 1.

Launch a sequence during which configurations are generated according to a
sampling algorithm (MD or ARTn).

At each iteration step the extrapolation-grade $\gamma$ is evaluated.

    1. (a)

If $0<\gamma<\gamma_{max}$, the energy and forces of the configuration are
evaluated with MTP;

    2. (b)

if $\gamma_{0}<\gamma<\gamma_{max}$, the configuration is set aside for an
update of MTP parameters;

    3. (c)

else if $\gamma>\gamma_{max}$, energy and forces of the configuration are not
evaluated with MTP and the configuration is not kept for update. The sequence
is stopped and we go directly to the update step (step 3).

  2. 2.

Move on next to the iteration in the sequence (step 1).

  3. 3.

The model is updated, if at at least one configuration as been set aside for
an update of MTP (i) at the end of a sequence or (ii) at any moment during the
sequence if $\gamma>\gamma_{max}$.

  4. 4.

If there is an update, restart a new sequence (go to step 1), else stop if no
configuration with $\gamma>\gamma_{0}$ has been set aside during the
predefined maximum length of the sequence.

The moment tensor potential model update is defined as follows (see Fig. 1,
right-hand side):

  1. 1.

A selection is made from the set aside configurations (with
$\gamma>\gamma_{0}$) using MaxVol Podryabinkin and Shapeev (2017).

  2. 2.

Each selected configuration is evaluated by the reference model

  3. 3.

The training set is updated with the new evaluated configurations

  4. 4.

The moment tensor potential is fitted on the new training set accordingly to
Eq. 4

More details of this procedure can be found in Ref. Podryabinkin and Shapeev,
2017.

Figure 1:  On-the-fly machine learning workflow used with MD and ARTn (on the
left). A potential update can take place at two points: when the sequence ends
or when $\gamma>\gamma_{max}$. The updating procedures are given in the box on
the right.

###  II.4 MD and ARTn

Two sampling approaches are used to generate a sequence of configurations: (1)
molecular dynamics (MD) as implemented within LAMMPS Thompson  _et al._ (2022)
and (2) the activation-relaxation technique nouveau (ARTn) algorithm developed
by Mousseau and collaborators Barkema and Mousseau (1996); Malek and Mousseau
(2000a); Jay  _et al._ (2022b). Since MD is well known, we only give below a
brief summary of ARTn.

ARTn is designed to explore the potential energy landscape of atomic systems
through the identification of local transition states connecting nearby local
minima. Its workflow can be summarized in three main steps (see, for a recent
in depth discussion of the ARTn version used in this work, see Ref. Jay  _et
al._ , 2022b):

  1. 1.

Leaving the harmonic well: starting from an energy minimum, an atom and its
neighbours are moved iteratively in a direction selected at random until a
direction of negative curvature on the potential energy surfaces,
$\mathbf{d}(\lambda_{\mathrm{min}})$ with $\lambda_{\mathrm{min}}$, the lowest
eigenvalue of the Hessian matrix, smaller than zero, emerges; this indicates
the presence of a nearby first-order saddle point;

  2. 2.

Converging to a first-order saddle point: the system is then pushed in the
direction of negative curvature $\mathbf{d}(\lambda_{\mathrm{min}})$ while the
force is minimized in the perpendicular plane, until the total force $F$
passes below a threshold near $F_{0}$, which indicates the saddle point have
been reached;

  3. 3.

Relaxing into a new minimum: the system is then pushed over the saddle point
and relaxed into a connected new minimum.

At each step $\lambda_{\mathrm{min}}$ and $\mathbf{d}(\lambda_{\mathrm{min}})$
are found using an iterative Lanczos method Lanczos (1950); Malek and Mousseau
(2000b); Jay  _et al._ (2022a). Perpendicular relaxation during activation and
global minimization are done using the Fast Inertial Relaxation Engine (FIRE)
algorithm Bitzek  _et al._ (2006).

Generated events are accepted or rejected according to the Metropolis
algorithm, where the acceptation probability $p$ is given by

$p=\min\left(1,e^{-\beta\Delta E}\right),$ (9)

with $\Delta E=E_{\mathrm{saddle}}-E_{\mathrm{minimum}}$, the energy
difference between the saddle and a connected minima and $\beta=1/k_{B}T$
where $k_{B}$ is the Boltzmann factor and $T$ is a fictitious temperature,
since thermal deformations are not taken into account. Potential energy
landscape exploration consist of generating a number of event.

###  II.5 Systems studied

The fitting approaches are tested on two physical systems: (i) a Si diamond
structure with Stillinger-Weber as a reference potential Stillinger and Weber
(1985); and (ii) a SiGe zincblende structure using the Stillinger-Weber
potential with parameters from Ref. Ethier and Lewis (1992). Both models count
215 atoms and a vacancy. These potentials as selected as they are
computationally light and facilitate the accumulation of statistics at the
level needed here; to obtain physically-relevant results, one should train on
DFT and not on empirical potentials.

The Si system is fitted with a ML potential set at level 16, with 92 moment
tensor functions ($B(\mathfrak{R})$, Eq. 1). For SiGe, a potential at this
level (16) generates errors on the barrier of the order of 0.5 eV, which
indicates that a richer set of parameters is needed to describe the chemical
diversity and a level 20 is chosen for this system, with 288 moment tensor
functions. The relation between the number of moment tensor functions for Si
and energy error is presented in Supplemental Fig. 1.

###  II.6 Fitting approaches

To evaluate the reliability of the various on-the-fly approaches to reproduce
the reference potential on configurations of interest for complex materials,
the training set is limited to structures visited during MD or ARTn
simulations within the conditions described below. No additional information
regarding alternative crystalline structures, defects, surfaces, pressure,
etc. is provided.

For each of these two systems, we compare the following approaches:

  1. 1.

ML-MD: The MTP potential is train OTF on MD simulations. The potential is then
evaluated, without further update, in ARTn simulation.

  2. 2.

OTF-MDART: Starting from the ML-MD generated potential, the MTP is re-trained
following the OTF procedure during ARTn simulations.

  3. 3.

OTF-ART: Training of the potential is done uniquely during ARTn runs with OTF.

The ML-MD approach is in line with Ref. Novoselov  _et al._ , 2019 where a
potential is trained OTF during MD. However, while the potential is trained
with MD, its accuracy is evaluated during ARTn activated process search.

####  II.6.1 ML-MD: simulations details

Nine sets of MTP ML-MD potentials are developed and trained independently
during NVT MD simulations. Each set is trained at one specific simulation
temperature ranging from 300 K to 2700 K by step of 300 K and starting from
the same 215 atom crystalline structure with a vacancy. Each set consists of
ten independently constructed MTP potentials for statistical purpose.

Training takes place on a series of sequences, each run for a maximum of 100
ps, with steps of 1 fs, with an average of 75 ps per cycle. MTP potentials
require about $34\pm 14$ and $93\pm 43$ learning cycles for Si and SiGe to be
converged: the MTP potential is considered having learned the potential when
no configuration generated during a 100 ps second is found in the
extrapolating zone of the potential (with $\gamma>\gamma_{0}$).

As long as this is not the case, the sequence is restarted from the same
initial structure with different initial velocities. To facilitate
convergence, ML-MD potentials are fitted over three sets of progressively more
restricted reliability extrapolation parameter $\gamma_{0}$. Moreover because
MD leads to global deformation, the extrapolation is computed using global
descriptors (see tab. 1).

The final potential is then evaluated, in a fixed form, in ARTn simulations.

Table 1: Extrapolation and selection hyper-parameter values used for the three on-the-fly approaches used in this work. approach: |  | $\gamma_{0}$ | $\gamma_{max}$ |  grade- mode
---|---|---|---|---
ML-MD |  | 5.5/3.3/1.1 | 60/10/2.2 | global
OTF-MDART |  | 1.1 | 2.2 | local
OTF-ART |  | 1.1 | 2.2 | local

####  II.6.2 OTF ARTn simulations details

Each ARTn simulation is launched for 1500 events, with 24 parallel independent
searches, for a total of 36 000 generated events. For ARTn, a sequence is
either a search for a saddle point (successful or failed) or a minimization
from the saddle to minimum.

At each point, 24 sequences are generated in parallel, and the configuration
selected for an update of the potential is made on the combined set of
configurations to generate one training set. Sequence are restarted from the
last accepted position or, in the case of the vacancy in Si, the ground state.
When an activation step generates a configuration with
$\gamma(\mathfrak{R})>\gamma_{max}$, it is relaunched with the same initial
deformation. As with MD, ten independent ARTn runs are launched for
statistics.

In the bulk, diffusion of the vacancy in Si takes place through a symmetric
mechanism bringing the vacancy from one state to an identical one so all ARTn
event searches are effectively started from the same state. Starting from a
zincblende structure, SiGe evolves according to an accept-reject Metropolis
with a fictitious temperature of $0.5$ eV Mousseau and Barkema (1999). Since
the configurations explored by ARTn are locally deformed; the extrapolation
grade for ARTn generated configurations used for the OTF-MDART and OTF-ART
approaches are evaluated with the local descriptors.

###  II.7 Analysis

Following the standard approach, the error is computed on the energy and force
differences between the MLP and reference potentials computed on the same
structures. Here, however, this error is only measured on configurations
generated during the ARTn procedure.

For the energy:

$\Delta E=|E_{MLP}(X_{MLP})-E_{ref}(X_{MLP})|,$ (10)

and, for the forces:

$\Delta
F=\frac{1}{N}\sum_{i=0}^{N}\sqrt{\|\mathbf{f}_{MLP}^{(i)}(X_{MLP})-\mathbf{f}_{ref}^{(i)}(X_{ref})\|^{2}},$
(11)

where the positions $X_{MLP}$ are obtained from a simulation run with the
machine-learned potential and the energy on this exact configuration is
computed with the reference and the machine-learned potentials. The same is
done for the error on forces.

Since this work is focused on the correct description of first-order
transition states, we also compute the minimum and saddle barrier positions
and energy convergence errors ($\Delta X_{\text{conv}}$, $\Delta
E_{\text{conv}}$) as

$\displaystyle\Delta X_{\text{conv}}$
$\displaystyle=\sqrt{\sum_{i=0}^{N}\|\mathbf{x}^{(i)}_{MLP}-\mathbf{x}^{(i)}_{ref}\|^{2}},$
(12) $\displaystyle\Delta E_{\text{conv}}$
$\displaystyle=|E_{MLP}(X_{MLP})-E_{ref}(X_{ref})|,$ (13)

where $X_{MLP}$ and $X_{ref}$ are the positions corresponding to minimum or
saddle point as defined by the MLP and the reference potentials respectively,
with $E_{MLP}(X_{MLP})$ and $E_{ref}(X_{ref})$ the corresponding energies; by
definition, forces are zero at these points defined by the respective
potentials.

While $X_{MLP}$ and $E_{MLP}(X_{MLP})$ are obtained on the ARTn trajectories,
$X_{ref}$ and $E_{ref}(X_{ref})$ are obtained after reconverging the minima or
the saddle point using the reference potential starting from $X_{MLP}$ and
following the ARTn procedure.

From an energy barrier $\delta E(X)$, the energy barrier error $\Delta\delta
E_{barrier}$ is given by

$\Delta\delta E_{barrier}=|\delta E_{MLP}(X_{MLP})-\delta E_{ref}(X_{ref})|.$
(14)

If no trend is observed between the different temperatures where potentials
are trained, we calculate their average and deviation in order to to
effectively compare them with other approach.

##  III Results

Figure 2:  Number of calls to the reference potential for each of the machine-
learned potentials developed for Si as function of the temperature referring
to the one used during MD training. Since configurations are relaxed to zero K
in ARTn simulations, there is no associated temperature for this procedure.
Vertical bars represent the standard deviation computed on ten independent
realisations. Figure 3:  Average energy (top) and mean absolute forces
(bottom) errors per atom for Si measured over all configurations generated
along pathways in ARTn for the three approaches. Temperature refers to the one
used during MD training. Vertical bars represent the standard deviation
computed on ten independent realisations. Figure 4:  Average energy barrier
error for Si as defined by Eq. 14 for all events generated in ARTn for the
three approaches. Temperature refers to the one used during MD training.
Vertical bars represent the standard deviation computed on ten independent
realisations. Figure 5:  Mean position error on all saddle point for Si.
Temperature refers to the one used during MD training. Vertical bars represent
the standard deviation computed on ten independent realisations.

In this section, we first examine results for a vacancy in c-Si to establish
the methods then consider the same approaches on the more complex SiGe alloy.

###  III.1 ML-MD

The ML-MD approach serves as a benchmark to assess the efficiency of the
various approaches in sampling energy barriers and diffusion mechanisms. Here,
ten independent ML potentials are generated through on-the-fly MD simulations
at 9 different target temperatures ranging from 300 to 2700 K by step of 300 K
and require between $253\pm 60$, at 300 K, and $369\pm 85$ evaluations of the
reference potential, at 2700 K, to complete learning cycles (see Fig. 2).

For the purpose of this work, the quality of the ML-MD potential is evaluated
on configurations generated with ARTn as local activated events associated
with vacancy in a crystalline environment are generated. To avoid non-physical
results, when a ARTn-generated configuration shows a $\gamma>200$, the
configuration is rejected, the event search is stopped and a new event search
is launched from the same initial minimum.

Fig. 3 shows the standard validation error on energy and forces calculated
over all configurations generated along pathways for the 36 000 successful
events and 10 080 failed saddle searches (a success rate of 78 %). The error
on energy increases almost exponentially with the sampling temperature,
ranging from $0.44\pm 0.36$ meV/atom at 300 K to $5.1\pm 1.7$ meV/atom at
2700K. The error on forces is essentially constant at 0.0123 eV/Å, on average,
between 300 and 1800 K, and increases rapidly at high temperature, to reach
0.0256 eV/Å at 2700 K.

Since, the focus of this work is on transition states, Fig. 4 displays the
error on the energy barriers as a function of MD-fitting temperature, computed
with Eq. 12 and averaged over all generated barriers. This error is relatively
uncorrelated of the MD temperature simulation with an average of $0.056\pm
0.022$ eV, with minimum error of $0.024\pm 0.01$ eV at 2400 K and maximum of
$0.08\pm 0.03$ eV at 1200 K. This error is lower than that for a general point
on the energy landscape (Fig. 3) in part because it is computed as a
difference between saddle and initial minimum.

Errors on the position of the saddle point, associated with the capacity to
reproduce correctly their geometry, are given in Fig. 5. The top panel
indicates the average distance between saddle points converged with the
reference and the ML potentials: it decreases from $0.16\pm 0.05$ Å at 300 K
to a minimum of $0.09\pm 0.02$ Å between 1500 and 2100 K, going up at the two
highest temperatures (2400 and 2700 K).

Overall, this straightforward fitting approach based on constant-temperature
MD runs provides accurate diffusion barriers, ranging from 0.51 to more than 4
eV, for a vacancy in crystalline silicon at a low computational costs (263 to
369 evaluations of the reference potential).

###  III.2 Revisiting ML-MD potential in ARTn: the OTF-MDART adjusting
approach

To evaluate the possibility of improving on ML-MD potentials for activated
events, potentials are on-the-fly re-trained during ARTn learning cycles (OTF-
MDART). Fig. 2 gives the number of calls to the reference potential for this
procedure during the ARTn runs (dashed orange line) as well as the total
number of calls, including those made during ML-MD fitting (solid orange
line). The number of calls during ARTn learning cycles ranges from $979\pm
153$ at 300 K to to $136\pm 38$ at 2700 K for a total of $1232\pm 177$ to
$505\pm 109$ respectively, when including ML-MD calls.

The error on energy and forces remains correlated with the ML-MD temperature:
it is higher when the error is higher at ML-MD trained temperature. This
correlation is particularly strong when retraining MD potentials fitted
between 1500 and 2700 K (Fig. 3, solid orange line). Error on energy for OTF-
MDART is almost constant between 300 and 2400 K, at 0.22 meV/atom, rising to
1.9 meV/atom at 2700 K, lower by 50 to 63 % than ML-MD. As similar improvement
is observed on the forces, which range from 0.0103 eV/Å, on average, between
300 and 1800 K, increasing to 0.0173 eV/Å at 2700 K, representing a 16 % to 32
% decrease in error.

Table 2:  Average energy barrier error and mean position error on all saddle point for Si. The average error for ML-MD and OTF-MDART training is taken over all temperature sets. The standard deviation computed on all temperature sets. Errors |  | ML-MD | OTF-MDART | OTF-ART
---|---|---|---|---
$\Delta\delta E_{barrier}$ (eV) |  |  0.056$\pm$0.022 |  0.040$\pm$0.012 |  0.039$\pm$0.008
$\Delta X_{\text{conv}}$ (Å) |  |  0.114$\pm$0.029 |  0.072$\pm$0.010 |  0.072$\pm$0.006

Between 300 and 1500 K, retrained potentials with OTF-MDART show more constant
energy barrier errors than pure ML-MD models (Fig. 4), with an error of about
$0.036$ eV (OTF-MDART) vs average of $0.064$ eV (ML-MD) a 44 % improvement. At
the highest temperature — 1800 to 2700 K, however, as OTF-MDART calls for less
learning cycles, errors and fluctuations are not reduced with respect to ML-
MD. Interestingly, though, improvements on the saddle position is observed at
all temperatures for OTF-MDART (Fig. 5) with an average error of $0.072\pm
0.010$ Å.

Overall, by retraining ML-MD potential in ARTn, errors are reduced and results
are more consistent, i.e., error distributions are narrower, irrespective of
the temperature used in the initial MD training. This additional retraining
leads to a 50 % to 96 % decrease in energy error (Fig. 3), a 29 % improvement
for average energy barrier errors (Tab. 2) and a 37 % reduction on mean saddle
positions errors but with an additional number of calls to the reference
potential increasing between 37 to 490 %.

Figure 6:  Fraction of original MD configurations (left scale) and total
number of MD configurations (right scale) remaining in the final training set
(TS) for Si. Temperature refers to the one used during MD training. Vertical
bars represent the standard deviation computed on ten independent
realisations.

These results can be understood by looking at the fraction of MD-generated
configurations that remain in the training set at the end of the simulation
(Fig. 6): at temperatures between 300 and 1200 K, none of the ML-MD
configurations remain in the final training set; this proportions goes from
from 1.3 to 38 % between 1500 and 2700 K (left-hand axis, blue line). At these
temperatures, the system melts and generates a wider range of configurations.
Since these configurations are far from ARTn-generated configurations, the
selection algorithm keeps them in the set even though they do not help reduce
errors for the configurational space of interest with ARTn.

###  III.3 The OTF-ART adjusting approach

Given the results for OTF-MDART, we now turn to an OTF approach entirely
integrated in ARTn, in an attempt to increase accuracy, and reduce the cost
and waste of evaluations of the reference potential.

Ten independents on-the-fly ML potential are generated entirely in ARTn for a
total of 36 000 events starting from the same initial minimum. Each potential
is trained initially from the same one configuration (the initial minimum), in
the training set. Each parallel event search goes trow a learning cycle if
needed and as the simulation progresses learning cycle become rarer. The
values are averaged over the ten simulation and as the simulation go through
learning.

With an average total of $628\pm 283$ reference potential evaluations, the
cost of the OTF-ART is between that of ML-MD and OTF-MDART. Along pathways,
the average energy error for these potentials is of $0.22\pm 0.03$ meV/atom,
on par with OTF-MDART potential based on low-temperature ML-MD fitting, and 49
% lower than the 300 K ML-MD potential. Errors on forces, at 0.011$\pm$0.001
eV/Å, are in between ML-MD (0.012 eV/Å) and OTF-MDART (0.010 eV/Å) at low
training temperature. Comparing with the 2700 K potential fitting in MD, OFT-
ART error is 57 % lower than ML-MD (0.026 eV/Å) and 36 % lower than OTF-MDART
(0.017 eV/Å).

Focusing on barrier energy, the average error is $0.039\pm 0.008$ eV (see Fig.
4), about 2.5 % lower than OTF-MDART and 30.3 % better than ML-MD. The error
of $0.072\pm 0.006$ Å on the converged saddle position is similar to the
$0.072\pm 0.010$ Å obtained with OTF-MDART and 37 % lower than with ML-MD
($0.114$ Å).

###  III.4 Reproducing the dominant diffusion mechanism

Figure 7:  ARTn-generated energy barrier distributions for vacancy-diffusion
events in Si, including direct barrier (from ground state) and inverse
barriers (from excited states), as generated with the MTP model (orange) and
re-converged using the reference model (blue) from saddles and minima position
originally found with the MTP model.

The exploration of the energy landscape around the vacancy leads to a
generation of wide range of activated mechanisms and associated barriers (both
forward, associated with the diffusion of the vacancy, and backward, from the
final minima back to the saddle point). Fig. 7 presents the complete
distribution of generated direct and inverse barriers connected to the ground
state. The peak near 0 eV (around $10^{-2}$ to $10^{-1}$ eV) is associated
with the inverse barrier to to direct saddle at 2.38, 2.70 eV and higher (up
to 5.5 eV), except for the inverse 0.45 eV barrier which is linked to the 2.87
eV direct barrier. Direct barriers at 0.51 eV represent symmetric first
neighbor vacancy diffusion while barriers at 2.38 and 2.70 eV are associated
with more complex vacancy diffusion mechanism El-Mellouhi, Mousseau, and
Ordejón (2004). Events with barriers at 2.38, 2.70 eV, for example, involve a
vacancy diffusion through complex bond-exchanges. Spectator events Kumeda,
Wales, and Munro (2001) where the diamond network around the vacancy is
transformed by a bond switching are also generated. This mechanism was
proposed by Wooten, Winer, and Weaire (WWW) to describe the amorphization of
silicon Wooten, Winer, and Weaire (1985). The main spectator event occurs as
two neighbors of the vacancy are pushed together allowing the creation of a
bound associated with the 2.87 eV barrier. Other mechanisms involve strong
lattice distortion and bond formation not involving direct neighbors of the
vacancy with very high energy barriers El-Mellouhi, Mousseau, and Ordejón
(2004) of in between 3.2 and 4.0 eV.

Table 3: Average energy barrier errors and mean saddle position error on the 0.51 eV vacancy diffusion for Si. The average error for ML-MD and OTF-MDART training is taken over all temperature sets. The standard deviation computed on all temperature sets. Errors |  | ML-MD | OTF-MDART | OTF-ART
---|---|---|---|---
$\Delta\delta E_{barrier}$ (eV) |  |  0.026$\pm$0.015 |  0.022$\pm$0.011 |  0.019$\pm$0.005
$\Delta X_{\text{conv}}$ (Å) |  |  0.088$\pm$0.036 |  0.040$\pm$0.017 |  0.047$\pm$0.018

Since vacancy diffusion for this system is dominated by a 0.51 eV single
barrier mechanism, with the next barrier at 2.35 eV, an accurate description
of the dominant mechanism is essential to correctly capture defect kinetics in
Si. Tab. 3 presents the error on this barrier for the three approaches
described above. With an error of 0.019$\pm$0.005 eV, a relative error of 3.7
%, OTF-ART offers the closest reproduction of the reference barrier, followed
by OTF-MDART and ML-MD, with a respective error of 0.022$\pm$0.011 (relative
error of 4.3 %) and 0.026$\pm$0.015 (5.1 %). Overall, the error on energy
barrier is lower than that on the total energy presented above
(0.046$\pm$0.006 eV for OTF-ART, for example), due to a partial error
cancellation associated with energy difference taken to measure the barrier.

The validity of the barrier is also measured by the precision on the saddle
geometry. For the 0.51 eV barrier, ML-MD converges with an error on the
position of 0.088$\pm$0.036Å, with OTF-MDART and OTF-ART giving error almost
50 % lower, at 0.040$\pm$0.017Å and 0.047$\pm$0.018Å respectively.

###  III.5 SiGe system

Figure 8:  SiGe barrier histogram, including direct barrier (from a states accepted) and inverse barriers (from excited states), as found on-the-fly by the MTP model(orange) and re-converge by the reference model(blue) from saddles and minimums position originally given by MTP. Figure 9:  Number of calls to the reference potential for each of the OTF machine-learned potentials developed for SiGe as a function of the temperature referring to the one used during MD training. Since configurations are relaxed to zero K in ARTn simulations, there is no associated temperature for this procedure. Vertical bars represent the standard deviation computed on ten independent realisations. Figure 10:  Average energy (top) and mean absolute forces (bottom) errors for SiGe measured over all configurations generated along pathways in ARTn for the three approaches. Temperature refers to the one used during MD training. Vertical bars represent the standard deviation computed on ten independent realisations. Table 4: Average energy barrier errors and mean saddle position error on all barriers for SiGe. The average error for ML-MD and OTF-MDART training is taken over all temperature sets. The standard deviation computed on all temperature sets. Errors |  | ML-MD | OTF-MDART | OTF-ART
---|---|---|---|---
$\Delta\delta E_{barrier}$ (eV) |  |  0.082$\pm$0.024 |  0.072$\pm$0.014 |  0.066$\pm$0.015
$\Delta X_{\text{conv}}$ (Å) |  |  0.091$\pm$0.020 |  0.076$\pm$0.013 |  0.070$\pm$0.014

Having shown the interest of developing a specific potential by applying on-
the-fly learning directly to activated events on a simple system such as c-Si
with a vacancy, we test this approach with a more complex alloy with the same
overall reference potential to facilitate comparison. Starting from a ordered
zincblende structure, the diffusion of a vacancy creates chemical disorder
that complexifies the landscape visited as shown by the continuous
distribution of activated barriers, including both direct and inverse
barriers, found as the vacancy diffuses (Fig. 8); we note that the lowest
barrier for a vacancy diffusing is around 0.6 eV, with lower barriers
associated, as for Si, with reverse jumps from metastable states. The energy
barrier distribution for a vacancy diffusing in SiGe (Fig. 8) is much more
complex than for Si due to the chemical disorder that builds as the vacancy
diffuses.

As stated in the methodology, the additional complexity of the system imposes
a richer machined-learning potential, with a larger set of parameters to
encompass the greater diversity in the components and the configurations, due
to chemical disorder. Combined, these two levels of complexity (set of
parameters and configurational) result in an overall higher numbers of calls
to the reference potential as compared to Si, irrespective of the approach
used (see Fig. 9 (SiGe) vs. Fig. 2(Si)): while ML-MD requires between 380
evaluations of the reference potential at 300 K and 1549 at 2700 K, OTF-MDART
needs a total of around 3465 calculations of the reference potential,
irrespective of the temperature as original ML-MD configurations are
progressively removed from the training set. This efforts results in a number
of calls to the reference potential for OTF-MDART 4 % higher than with OTF-ART
(3329 on average). To reduce computational costs, we omit the 1500 K run, as
statistical behavior is smooth in this temperature region.

Figure 11:  Percentage of search interruptions during ML-MD potential
evaluation in ARTn ($\gamma>200$) for Si and SiGe as function of ML-MD
training temperature. Vertical bars represent the standard deviation computed
on ten independent realisations.

To disentangle the two contributions, we compare with the cost of fitting a Si
potential with the same level 20 potential as used for SiGe. Following the
full OTF-ART procedure, creating a Si MLP requires 2926 calls to the reference
potential. The intrinsic complexity of the landscape contributes therefore to
about a 14 % increase of the Si baseline calls count. In terms of accuracy,
the Si MLP level 20 leads to an average error on energy of 0.1 meV/atom, about
50 % lower than with the level 16 potential described above (0.22 meV/atom).
For SiGe, this error is (0.42 meV/atom), two times higher than for Si MLP
level 16 and four times that of Si MLP level 20.

This can be understood by the number of different configurations visited: as
opposed to the Si system where each initial minimum is identical (as the
vacancy moves in an otherwise perfect elemental crystal), the binary system is
transformed as the vacancy diffuses, as the chemical order is slowly
destroyed: each of the 24 ARTn parallel trajectories used to define the
potential over 1500 events evolves independently according to a probability
given by the Metropolis algorithm with a fictitious temperature (since network
itself is structurally at 0K) of 0.5 eV (Eq. 9), providing a rich range of
local environments.

Fitting a potential is clearly harder: with the parameters used — when a
configuration graded at $\gamma>200$ is encountered, the ARTn event search is
stopped —, no event could be generated using the ML-MD potential at 300 K and
600 K, which explains the absence of data for this temperatures in Fig. 10 and
4. For SiGe, the error on energy (see Fig. 10) with the ML-MD at 900 K and
above ranges from 0.5 meV/atom to 1.4 meV/atom, as a function of temperature.
On average, these errors are between 14 % and 69 % lower with OTF-MDART or
OTF-ART at around 0.43 meV/atom.

The OTF-ART approach gives an error in energy barrier of $0.066\pm 0.015$ eV
which represent a 19.5 % and 8.3 % lower error from the ML-MD ($0.082\pm
0.024$ eV) and OTF-MDART ($0.072\pm 0.014$ eV) respectively (Tab. 4). The
errors on the converged saddle position for OTF-ART and OTF-MDART are similar
at $0.070\pm 0.014$ Å and $0.076\pm 0.013$ Å, respectively, and represent a 23
% lower error than with ML-MD ($0.092$ Å). This accuracy is similar to that
obtained with Si, in contrast to total energy and energy barrier errors.

We note that the advantage of ML-MD for SiGe is overstated as shown by the
proportion of events generated with ML-MD potential that are interrupted due
to a too large extrapolation grade, $\gamma>200$ for both SiGe and Si (Fig.
11): for SiGe between 85 % and 30 % of events are aborted between 300 K and
1200 K, respectively. This proportion falls to zero percent failure at 1800 K.

##  IV Discussion

We compare three approaches aimed at the construction of potentials with
machine learning on-the-fly for the exploration of activated mechanism of the
potential energy landscape. We evaluate these by computing their efficiency at
reproducing the energy landscape around a vacancy in two systems, a relatively
simple Si diamond system (Fig. 7) and a more complex SiGe zincblende system
that disorders under vacancy diffusion (Fig. 8), both described with a
reference empirical potential to allow for better statistical analysis.

The first approach, which sets the comparison level, constructs a more general
machine learning potential with molecular dynamics (ML-MD), the second on-the-
fly adjusts this this generated potential, during the search for activated
events using ARTn, while the third approaches constructs a specifically on-
the-fly trained a potential during search of activated events (OTF-ART). The
efficiency of these three procedures is measured on the quality of the
reproduction of the reference potential during the search for activated event.

The baseline, defined by the ML-MD, is competitive with previously published
work. Energy errors for the more standard ML MD approach with a level 16
potential range from $0.44\pm 0.36$ meV/atom at 300 K to $5.1\pm 1.7$ meV/atom
at 2700 K (Fig. 3), an order of magnitude lower or similar than the 4 meV/atom
on an MTP potential of level 24 for Si obtained by Zuo et al. Zuo  _et al._
(2020), with the difference explained by the fact that activated events
involve local deformations from a zero-temperature crystal with a vacancy and
that DFT potentials are more difficult to fit than empirical ones Novikov  _et
al._ (2020). That is, using of DFT, we would expect that these conclusions
hold as we are looking at difference between our three approaches, not
claiming specific accuracy of these.

Similarly, the relative energy error on the dominant 0.51 eV diffusion barrier
for SW Si is of 5.1 % (0.026 eV) with the ML-MD approach and 3.7 % (0.019 eV)
with the OTF-ART. Using the same MTP potential trained using an OTF MD with an
ab initio reference potential, Novoselov et al. find a 0.20 eV barrier for
vacancy diffusion in Si as compared with 0.18 eV with the reference potential,
an error of 0.02 eV or a 10.0 % relative error.

Overall, the ML-MD approach, especially when run at temperatures between 900
and 1800 K, can generate a generic ML potential with reasonable precision for
describing activated mechanisms in Si and SiGe. Developing a more specific OTF
potential, generated directly with ARTn on activated trajectories, however,
offers a more accurate description of both the energy and geometry at the
barriers.

It is possible to recover this precision by adjusting the original MD
potential during ARTn runs, however, this increases the number of calls to the
reference potential, raising the total costs beyond that of OTF-ART while
largely erasing work made during ML-MD training phase: for Si, between 300 and
1200 K, none of the ML-MD configurations are retained while around 1.3 to 12.5
% are retained for the potential trained in range of 1500 and 2700 K (Fig. 6,
right-hand axis, orange line), but at the cost of lowering the precision on
barriers.

Moving to a more complex system, such an evolving binary alloy, increases the
overall cost of the procedure in terms of calls to the reference potential, as
more parameters need to be fit. Here also, the gain on using a specific
potential constructed from ARTn trajectories is notable, both in the average
errors and their fluctuations. Indeed, the ML-MD potential presents
considerable instabilities while generating activated trajectories as it can
be see by the number of configurations considered out-side of the potential’s
scope ($\gamma>200$), see Fig. 11. This can also be see by comparing the
number of evaluations of the reference potential as function of the
temperature for OTF-MDART; for Si this number decreases while, for SiGe, the
total number of evaluations of the reference potential remains constant. This
mean that the potential requires significant re-adjustment, even from its
high-temperature training, to adapt to the created disorder of the evolving
SiGe system.

## Conclusion

We compare the advantage of using a more general vs specific machine-learned
potential (MLP) to describe activated mechanisms in solid. To do so, we
generate first an MLP constructed with the Moment Tensor Potential formalism
Novikov  _et al._ (2020); Shapeev (2016) to replicate Stillinger-Weber
potential for Si and SiGe crystals with a single vacancy using a standard
molecular dynamics procedure (MD-ML).

Comparing the quality of the reproduction of activated mechanisms with a ML
potential further refined during an activation-relaxation technique nouveau
sampling of the energy landscape and a potential unique constructed on-the-fly
within ARTn, we show that while a general potential can deliver high accuracy
for both the barrier geometries and their related energies, error and
fluctuations around the average value are significantly lowered by
constructing a specific potential, with a number of calls to the reference
potential that is lower than a combined approach (MD + ARTn) for a similar
precision.

The advantage of using a specific potential remains when looking at more
complex materials, such the SiGe alloys considered here, even though the
advantage in terms of calls to the reference is strongly reduced. The next
steps will involve applying this strategy with DFT as reference potential to
attack problems that have long been out of reach of computational materials
sciences, allowing a much closer connection between modeling and experience.

##  V Supplementary Material

See the supplementary material for a figure of the relation between accuracy
and the level of the MTP potential for Si. Also present are figures displaying
the accuracy of the saddle position for the 0.51 eV barrier in Si and all
barriers in SiGe.

##  VI Code and data availability

The ARTn packages as well as the data reported here are distributed freely.
Please contact Normand Mousseau
(normand.mousseau@umontreal.ca).

###### Acknowledgements.

This project is supported through a Discovery grant from the Natural Science
and Engineering Research Council of Canada (NSERC). Karl-Étienne Bolduc is
grateful to NSERC and IVADO for summer scholarchips. We are grateful to Calcul
Québec and Compute Canada for generous allocation of computational resources.

## References

  * Kohn and Sham (1965) W. Kohn and L. J. Sham, “Self-consistent equations including exchange and correlation effects,” Phys. Rev. 140, A1133–A1138 (1965).
  * Lindahl (2008) E. R. Lindahl, “Molecular dynamics simulations,” in  _Molecular modeling of proteins_ (Springer, 2008) pp. 3–23.
  * Voter and Doll (1984) A. F. Voter and J. D. Doll, “Transition state theory description of surface self-diffusion: Comparison with classical trajectory results,” The Journal of chemical physics 80, 5832–5838 (1984).
  * Voter (2007) A. F. Voter, “Introduction to the kinetic monte carlo method,” in  _Radiation effects in solids_ (Springer, 2007) pp. 1–23.
  * Henkelman and Jónsson (2001) G. Henkelman and H. Jónsson, “Long time scale kinetic monte carlo simulations without lattice approximation and predefined event table,” The Journal of Chemical Physics 115, 9657–9666 (2001), https://doi.org/10.1063/1.1415500 .
  * El-Mellouhi, Mousseau, and Lewis (2008) F. El-Mellouhi, N. Mousseau, and L. J. Lewis, “Kinetic activation-relaxation technique: An off-lattice self-learning kinetic Monte Carlo algorithm,” Physical Review B 78, 153202 (2008).
  * Behler and Parrinello (2007) J. Behler and M. Parrinello, “Generalized Neural-Network Representation of High-Dimensional Potential-Energy Surfaces,” Phys. Rev. Lett. 98, 146401 (2007).
  * Bartók, Kondor, and Csányi (2013) A. P. Bartók, R. Kondor, and G. Csányi, “On representing chemical environments,” Physical Review B 87, 184115 (2013).
  * Thompson  _et al._ (2015) A. P. Thompson, L. P. Swiler, C. R. Trott, S. M. Foiles, and G. J. Tucker, “Spectral neighbor analysis method for automated generation of quantum-accurate interatomic potentials,” Journal of Computational Physics 285, 316–330 (2015).
  * Shapeev (2016) A. V. Shapeev, “Moment tensor potentials: A class of systematically improvable interatomic potentials,” Multiscale Modeling & Simulation 14, 1153–1173 (2016).
  * Sivaraman  _et al._ (2021) G. Sivaraman, J. Guo, L. Ward, N. Hoyt, M. Williamson, I. Foster, C. Benmore, and N. Jackson, “Automated development of molten salt machine learning potentials: application to licl,” The Journal of Physical Chemistry Letters 12, 4278–4285 (2021).
  * Kang, Shang, and Liu (2020) P.-L. Kang, C. Shang, and Z.-P. Liu, “Large-scale atomic simulation via machine learning potentials constructed by global potential energy surface exploration,” Accounts of Chemical Research 53, 2119–2129 (2020).
  * Sivaraman  _et al._ (2020) G. Sivaraman, A. N. Krishnamoorthy, M. Baur, C. Holm, M. Stan, G. Csányi, C. Benmore, and Á. Vázquez-Mayagoitia, “Machine-learned interatomic potentials by active learning: amorphous and liquid hafnium dioxide,” npj Computational Materials 6, 1–8 (2020).
  * Béland  _et al._ (2013) L. K. Béland, Y. Anahory, D. Smeets, M. Guihard, P. Brommer, J.-F. Joly, J.-C. Pothier, L. J. Lewis, N. Mousseau, and F. Schiettekatte, “Replenish and Relax: Explaining Logarithmic Annealing in Ion-Implanted $c$-Si,” Phys. Rev. Let. 111, 105502 (2013).
  * Osetsky  _et al._ (2018) Y. N. Osetsky, L. K. Beland, A. V. Barashev, and Y. Zhang, “On the existence and origin of sluggish diffusion in chemically disordered concentrated alloys,” Current Opinion in Solid State and Materials Science 22, 65–74 (2018).
  * Restrepo  _et al._ (2018) O. A. Restrepo, N. Mousseau, M. Trochet, F. El-Mellouhi, O. Bouhali, and C. S. Becquart, “Carbon diffusion paths and segregation at high-angle tilt grain boundaries in $\alpha$-fe studied by using a kinetic activation-relation technique,” Physical Review B 97, 054309 (2018).
  * Barkema and Mousseau (1996) G. T. Barkema and N. Mousseau, “Event-Based Relaxation of Continuous Disordered Systems,” Physical Review Letters 77, 4358–4361 (1996).
  * Malek and Mousseau (2000a) R. Malek and N. Mousseau, “Dynamics of Lennard-Jones clusters: A characterization of the activation-relaxation technique,” Physical Review E 62, 7723–7728 (2000a).
  * Jay  _et al._ (2022a) A. Jay, M. Gunde, N. Salles, M. Poberžnik, L. Martin-Samos, N. Richard, S. d. Gironcoli, N. Mousseau, and A. Hémeryck, “Activation–Relaxation Technique: An efficient way to find minima and saddle points of potential energy surfaces,” Computational Materials Science 209, 111363 (2022a).
  * Truhlar, Garrett, and Klippenstein (1996) D. G. Truhlar, B. C. Garrett, and S. J. Klippenstein, “Current status of transition-state theory,” The Journal of physical chemistry 100, 12771–12800 (1996).
  * Novikov  _et al._ (2020) I. S. Novikov, K. Gubaev, E. V. Podryabinkin, and A. V. Shapeev, “The mlip package: moment tensor potentials with mpi and active learning,” Machine Learning: Science and Technology 2, 025002 (2020).
  * Novoselov  _et al._ (2019) I. Novoselov, A. Yanilkin, A. Shapeev, and E. Podryabinkin, “Moment tensor potentials as a promising tool to study diffusion processes,” Computational Materials Science 164, 46–56 (2019).
  * Stillinger and Weber (1985) F. H. Stillinger and T. A. Weber, “Computer simulation of local order in condensed phases of silicon,” Physical review B 31, 5262 (1985).
  * Ethier and Lewis (1992) S. Ethier and L. J. Lewis, “Epitaxial growth of si1- xgex on si (100) 2$\times$ 1: A molecular-dynamics study,” Journal of materials research 7, 2817–2827 (1992).
  * Zuo  _et al._ (2020) Y. Zuo, C. Chen, X. Li, Z. Deng, Y. Chen, J. Behler, G. Csányi, A. V. Shapeev, A. P. Thompson, M. A. Wood, _et al._ , “Performance and cost assessment of machine learning interatomic potentials,” The Journal of Physical Chemistry A 124, 731–745 (2020).
  * Podryabinkin and Shapeev (2017) E. V. Podryabinkin and A. V. Shapeev, “Active learning of linearly parametrized interatomic potentials,” Computational Materials Science 140, 171–180 (2017).
  * Podryabinkin  _et al._ (2019) E. V. Podryabinkin, E. V. Tikhonov, A. V. Shapeev, and A. R. Oganov, “Accelerating crystal structure prediction by machine-learning interatomic potentials with active learning,” Physical Review B 99, 064114 (2019).
  * Gubaev  _et al._ (2019) K. Gubaev, E. V. Podryabinkin, G. L. Hart, and A. V. Shapeev, “Accelerating high-throughput searches for new alloys with active learning of interatomic potentials,” Computational Materials Science 156, 148–156 (2019).
  * Thompson  _et al._ (2022) A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in ’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott, and S. J. Plimpton, “LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales,” Comp. Phys. Comm. 271, 108171 (2022).
  * Jay  _et al._ (2022b) A. Jay, M. Gunde, N. Salles, M. Poberžnik, L. Martin-Samos, N. Richard, S. de Gironcoli, N. Mousseau, and A. Hémeryck, “Activation–relaxation technique: An efficient way to find minima and saddle points of potential energy surfaces,” Computational Materials Science 209, 111363 (2022b).
  * Lanczos (1950) C. Lanczos, “An iteration method for the solution of the eigenvalue problem of linear differential and integral operators,” (1950).
  * Malek and Mousseau (2000b) R. Malek and N. Mousseau, “Dynamics of lennard-jones clusters: A characterization of the activation-relaxation technique,” Physical Review E 62, 7723 (2000b).
  * Bitzek  _et al._ (2006) E. Bitzek, P. Koskinen, F. Gähler, M. Moseler, and P. Gumbsch, “Structural relaxation made simple,” Physical review letters 97, 170201 (2006).
  * Mousseau and Barkema (1999) N. Mousseau and G. T. Barkema, “Exploring High-Dimensional Energy Landscapes,” Computing in Science & Engineering 1, 74–82 (1999).
  * El-Mellouhi, Mousseau, and Ordejón (2004) F. El-Mellouhi, N. Mousseau, and P. Ordejón, “Sampling the diffusion paths of a neutral vacancy in silicon with quantum mechanical calculations,” Phys. Rev. B 70, 205202 (2004).
  * Kumeda, Wales, and Munro (2001) Y. Kumeda, D. J. Wales, and L. J. Munro, “Transition states and rearrangement mechanisms from hybrid eigenvector-following and density functional theory.: application to c10h10 and defect migration in crystalline silicon,” Chemical physics letters 341, 185–194 (2001).
  * Wooten, Winer, and Weaire (1985) F. Wooten, K. Winer, and D. Weaire, “Computer generation of structural models of amorphous si and ge,” Phys. Rev. Lett. 54, 1392–1395 (1985).