# Beyond potential energy surface benchmarking: a complete application of
machine learning to chemical reactivity

Xingyi Guan1,3, Joseph Heindel1,3, Taehee Ko4, Chao Yang5, Teresa Head-
Gordon1-3

###### Abstract

We train an equivariant machine learning model to predict energies and forces
for a real-world study of hydrogen combustion under conditions of finite
temperature and pressure. This challenging case for reactive chemistry
illustrates that ML-learned potential energy surfaces (PESs) are always
incomplete as they are overly reliant on chemical intuition of what data is
important for training, i.e. stable or metastable energy states. Instead we
show here that a “negative design” data acquisition strategy is necessary to
create a more complete ML model of the PES, since it must also learn avoidance
of unforeseen high energy intermediates or even unphysical energy
configurations. Because this type of data is unintuitive to create, we
introduce an active learning workflow based on metadynamics that samples a
lower dimensional manifold within collective variables that efficiently
creates highly variable energy configurations for further ML training. This
strategy more rapidly completes the ML PES such that deviations among query by
committee ML models helps to now signal occasional calls to the external ab
initio data source to further molecular dynamics in time without need for
retraining the ML model. With the hybrid ML-physics model we predict the
change in transition state and/or reaction mechanism at finite temperature and
pressure for hydrogen combustion, thereby delivering on the promise of real
application work using ML trained models of an ab initio PES with two orders
of magnitude reduction in cost.

1Kenneth S. Pitzer Theory Center and Department of Chemistry,

2Departments of Bioengineering and Chemical and Biomolecular Engineering

University of California, Berkeley, CA, 94720 USA

3Chemical Sciences Division, Lawrence Berkeley National Laboratory, Berkeley,
CA, 94720 USA

4Department of Mathematics, Penn State University, University Park, PA 16802

5Applied Mathematics and Computational Research Division, Lawrence Berkeley
National Laboratory

corresponding author<EMAIL_ADDRESS>
##  1 Introduction

Machine learning (ML) methods are emerging as an alternative to ab initio
molecular dynamics (AIMD) and physical-based potential energy functions (or
force fields), once trained on high quality ab initio energies and forces
associated with a given conformation of nuclei. Starting with the generalized
neural network representation of high dimensional potential energy surfaces
(PESs) proposed by Behler and Parrinello1, their work inspired additional
state-of-the-art approaches for non-reactive systems, such as ANI-1 for
predicting energies of organic molecules2, 3, 4, kernel methods such as
Gradient Domain Machine Learning (GDML) that predict atomic forces directly5,
6, and the introduction of other invariant models based on message-passing
architectures developed as generic molecular property predictors7 A key
development more recently are ML models that are equivariant to translations
and rotations through their architectures, showing testing superiority in
accuracy benchmarks with greatly reduced quantities of reference data for
training.8, 9, 10, 11, 12, 13, 14, 15, 16 In what follows we use NewtonNet 16,
a physics inspired message passing equivariant neural network, as the
underlying deep learning model for energy and force prediction for the
surprisingly difficult case of hydrogen combustion, for which we developed the
HCombustion dataset 17 of energies and forces generated using the
$\omega$B97X-V 18 DFT functional with the cc-pVTZ basis set. In the case of
hydrogen combustion there are at least 19 reaction channels, multiple stable
and unstable intermediates that are dependent on a given reactive channel, and
complications can arise due to creation of radical species and alternative
spin states during the combustion process.

As such, and regardless of the architecture, the reliability of the ML model
still heavily relies on the diversity of the training data, especially for
chemically reactive systems that must visit high energy states when undergoing
chemical transformations. ML models by their nature interpolate between known
training data, but its extrapolation capability is limited, i.e. predictions
can be unreliable when molecules or their configurations are dissimilar to
those in the training set. Thus, in order to achieve meaningful chemically
reactive simulations with ML potentials, the training sets should cover a wide
range of structural space with variable energy stability.19 However, it is
challenging to formulate a priori a dataset that is balanced and diverse for a
given reactive system, and it is not unusual for the ML model to still suffer
from an overfitting problem that can lead to models with good accuracy on its
original test set but are error-prone when applied to MD simulations,
especially in the case for gas phase chemical reactivity in which energy
configurations are highly diverse.

Traditionally, active learning (AL) is a powerful strategy for reducing the
amount of labeled data required to develop an ML potential, and selecting
informative configurations for labeling.20 Typically AL uses a query by
committee strategy, in which variance among a set of identical architectures
but stochastically initiated ML models select the most informative data points
for labeling, reducing the data generation effort and improving the accuracy
of the ML model. However AL informed ML is still not a panacea without having
more information. In this study, we propose an active learning workflow for
chemical reactivity that utilizes a different information source - namely the
sampling efficiencies that are inherent in statistical mechanics methods for
rare events. In particular, we have formulated an AL workflow that expands on
the originally formulated HCombustion dataset 17 by formulating collective
variables (CVs) to first systematically sample a lower manifold of all the
intrinsic reaction coordinates of the 19 reactive channels, and then to
stochastically sample with metadynamics to take advantage of its known better
ergodicity.21, 22 It is mostly irrelevant whether the CVs are "good", in the
sense that the first purpose of the CVs is to trap both the intrinsic reaction
coordinate and to sample in the orthogonal directions for much higher energy
configurations. Within the smaller subspace metadynamics fills the free energy
wells to more rapidly find high energy states such that the ML model learns
about configurations whose energies are simply inaccessible, or to find new
reactive intermediates which become accessible with temperature and pressure
but are of the wrong energy ordering.

This active learning strategy allows us to reach a final hydrogen combustion
ML model that is more diverse and balanced, and more importantly, an energy
surface that is relatively smooth. We then utilize the variance of the ML
committee models to report back on additional but much fewer PES rough spots,
at which we substitute a direct call to an ab initio force to complete a local
in time molecular dynamics update, until the ML models recover accurate forces
to continue the trajectory, all without further retraining. We illustrate the
completeness of the ML hydrogen combustion model with a metadynamics
application using well formulated CVs inspired by diffusion maps23 to discover
the entropic contribution of free energy transition states for hydrogen
combustion reaction channels.

##  2 Results

###  2.1 Active Learning workflow

NewtonNet16 has been previously trained on the HCombustion dataset17(which we
refer to as the original dataset) comprised of ab initio molecular dynamics
(AIMD) sampling around the intrinsic reaction coordinate (IRC) at 0K near the
transition state, as well as systematic normal modes with respect to the IRC
to produce PES curvature data, all specific to 19 elemetary reactions
(Supplementary Table S1).24, 25 This yielded a total of $\sim$290,000
potential energies and $\sim$1,260,000 nuclear force vectors. Using $\sim$5k
data per reaction, the trained NewtonNet model reaches very good accuracy for
energies (mean absolute error (MAE) of 0.14 kcal/mol) and forces (MAE of 0.33
kcal/mol/Å).

However, when the ML model is applied to molecular simulation, we find it to
be highly error-prone, requiring for example new data corresponding to
atomization energies to understand the inherent stability of molecules.
Furthermore, while the model is trained on relatively stable geometries as
well as metastable transition states, it lacks the knowledge of any highly
unstable state, and thus predicts configurations with the wrong energy
ordering as well as unphysical geometries that were predicted by the L model
to be energetically stable. In Figure 1(a), we show representative structures
for these two cases; on the left is an unphysical geometry that is close to
dissociation, with DFT energy of -129.16 kcal/mol, but the ML model predicts
its energy to be as stable as -324.05 kcal/mol. On the right is the hydronium
ion that appears in one of the trajectory with an artificially low energy of
-295.56 kcal/mol predicted by model, while its actual DFT energy is -189.98
kcal/mol. Figure 1(b) shows such a trajectory, starting from the transition
state of rxn16, proceeding through a series of unphysical states, and
eventually ending in configurations where all the bonds are broken.

Hence an MD simulation driven with the model will quickly generate
trajectories that go to unphysical regions because the quantum mechanical data
will always be small-scale and insufficient for generating a complete
potential energy surface. It also emphasizes how unintuitive data acquisition
is in regards creating a robust ML model of the PES, especially for chemical
reactivity, which requires negative design principles, i.e. to add
unphysically high energy species to the training of the ML model. To make this
systematic and tractable we formulate an IRC dilation dataset within a lower
manifold of collective coordinates that encapsulates the IRC for each reaction
channel. Using the example of Rxn16 whose AIMD and normal mode data shown in
Figure 1c, each geometry along the IRC curve is proportionally scaled with
multiple ratios of CVs to generate new high energy species as shown in Figure
1d. This addition to the dataset helps the models to learn bond dissociation
and contraction, as well as spanning a more diverse chemical space than the
original dataset.

Figure 1: (a)Two representative structures that the original ML model predicts
with large error. (b) A trajectory for Rxn16 driven with the original ML model
visualized on two reaction coordinates CN1 (H4-H5) and CN2(O2-H5). (c,d) Rxn16
DFT potential energy surfaces visualized on two reaction coordinates CN1
(H4-H5) and CN2(O2-H5). (c) the original HCombustion dataset with AIMD and
normal mode data (d) the dialation dataset. Each CN represents the formation
or breaking of relevant bonds in the reaction process. Energy of each point is
color coded in (c) and (d), with red color meaning points with energy higher
than the Boltzmann weighting threshold.

Although the dilation data is helpful, the PES for hydrogen combustion remains
incomplete. We thus introduce an active learning workflow in which we trained
four models with 1000 structures per reaction selected from the HCombustion
dataset and 200 structures per reaction from the dilation dataset using the
same architecture but different initial parameters. These four models serve as
a starting point of an iterative process, outlined in Figure 2, to
systematically improve the ML model through active learning. To allow for a
fast turn around time, we use relatively small epoch size and large learning
rate in the training through the active learning rounds. More relevant here is
that we chose metadynamics, a rare events sampling method, to more efficiently
sample previously unseen and unstable structures through an external biasing
potential that forces the system to explore regions of high (free) energy.
While in the usual context poor selection of the low-dimensional descriptors
affects the rate at which transitions are enhanced, in this particular context
we are using metadynamics as a tool to fill in the holes in the ML PES in
which the goodness of the CVs is less important than the ability to sample
diverse conformations of high energy variance.

Figure 2: Schematic illustration of active learning workflow using query by
committee. The four NewtonNet models serve as a starting point of an iterative
process, where each round of active learning consists of the following steps:
(1) Perform several short metadynamic simulations to explore the configuration
space in a lower dimension. (2) When the four models disagree outside standard
deviation, collect a representative subset of structures to be included in the
training set through downsampling. (3) Perform DFT calculation of energies and
atomic forces. (4) Retrain the ensemble of ML models with the updated training
set. The details of each of the steps are described in the method section.

To illustrate the active learning approach, we consider Rxn18 as an example in
which the potential energy surface is projected onto two reaction coordinates
CN(O2-O5) and CN(O5-H4), with results shown in Figure 3. The ML model
performance was tracked by analyzing both the original data points derived
from AIMD and normal modes calculations,

Figure 3: Potential energy surface in collective coordinates (left) and change
in energy and force mean absolute error(MAE) as active learning round proceeds
(right) for Reaction 18. The potential energy surface projected onto two
reaction coordinate CN(O2-O5) and CN(O5-H4). Energy of each point is color
coded, with red color meaning points with energy higher than the Boltzmann
weighting threshold. (a,b) Test set of rxn18 in the original dataset composed
of data generated from AIMD and normal modes calculations. (c,d) Data added
between active learning round 1-20, sampled with 2000fs metadynamics. (e,f)
Data added between active learning round 21-33, sampled with 5000fs
metadynamics. (g,h) Data added between active learning round 34-48, sampled
with 10000fs metadynamics.

and the newly added data points accumulated during the active learning
procedure. Because we used longer metadynamics simulations for sampling as the
active learning rounds proceeded and as errors decreased, we show this by
dividing the active learning data into four batches: the original data derived
from normal mode and AIMD near the transition state of the IRC (Figure
3(a,b)), data sampled with 2000fs metadynamics between active learning round
1-20 (Figure 3(c,d)), data sampled with 5000fs metadynamics between active
learning round 21-33 (Figure 3(e,f)), data sampled with 10000fs metadynamics
between active learning round 34-48 (Figure 3(g,h)).

As shown on the left side of Figure 3, active learning spans more of the
reaction coordinate space and in favor of learning high energy states, which
indicates that our sampling method is effective in learning points to avoid in
the MD trajectory. The lack of points in the upper right quadrant in early
stages of active learning indicates that the dilation dataset was sufficient
(Figure 3c), and that the model became more accurate near the reaction pathway
basins since there is little sampling in this region during later rounds of
the active learning process using metadynamics (Figure 3(e,g)). The right side
of Figure 3 shows that there is a very significant improvement in model error
when predicting on the newly added data. Using the original ML model these new
sampled points from active learning created large energy and force errors as
high as $\sim$8 kcal/mol and $\sim$12 kcal/mol/Å respectively, which of course
would create an untenable MD trajectory if those geometries are visited during
an actual application. But after $\sim$50 rounds of active learning, the model
prediction on new geometries are hugely improved, with 0.97 kcal/mol/Å error
in forces on the data points collected between round 21-33, and 1.24
kcal/mol/Å force error on the data points collected between round 34-48.

The somewhat larger errors are still excellent given the much larger energy
range required to generate physically meaningful trajectories, which derives
from the more significant improvements of the new data points added with the
active learning workflow. It also reflects another important strategy in
active learning, which is that high energy configurations must be learned but
they do not have to be as accurate, and this can be reflected in the loss
function through Boltzmann weighting (see Methods). This is supported by the
fact that the original data points are still predicted with relatively small
errors in energy and forces with the active learning model ( 0.4 kcal/mol
energy MAE and 0.8 kcal/mol/Å force MAE).

###  2.2 Committer Analysis and the Free Energy Surface

To investigate the outcome of the active learning to create an ML potential
for hydrogen combustion, we performed committer analysis on the elementary
reaction channels for hydrogen combustion, starting the trajectories from the
IRC transition state of the reaction with temperature set to 500 K. Diatomic
reactions (rxn05, 06, 07 and 08) and barrierless reaction 15 were not
considered. Reaction 12 has two spin states doublet and quartet in the
original dataset, and only the transition state with lower energy (12b
quartet) was considered because the model was trained with energy and forces
from the lower energy spin state.

Table 1 presents the committer results for both the original model without
active learning and the final model after all active learning rounds. The
final ML model gives quite different committer statistics than the original
model, in which the AL model shows a more even chance to commit to reactant
and product for a majority of the reactions. Hence although the AIMD and NM
sampling provided data representation in the tube around the IRC, it is
insufficient such that the commitor analysis is qualitatively different in a
majority of the reactions. But once the model is more complete using our AL
strategy, there is a shift in the forward/backward committer distribution in
most of the reaction channels in which the 0 K IRC is still a good estimator.
Even so, at 500 K there is some evidence of a non-negligible entropy component
to the free energy transition state.

Table 1: Committer statistics using (a) the original model (b) the final model after active learning rounds | Original Model | Final Model
---|---|---
rxn | reactant(%) | product(%) | reactant(%) | product(%)
01 | 12 | 88 | 25 | 75
02 | 0 | 100 | 44 | 56
03 | 48 | 49 | 50 | 50
04 | 100 | 0 | 44 | 56
09 | 86 | 14 | 62 | 38
10 | 58 | 41 | 55 | 45
11 | 91 | 9 | 55 | 45
12 | 28 | 72 | 45 | 55
13 | 55 | 44 | 65 | 35
14 | 7 | 93 | 40 | 60
16 | 0 | 100 | 48 | 52
17 | 96 | 4 | 49 | 51
18 | 22 | 78 | 43 | 57
19 | 95 | 5 | 52 | 48

With a more complete ML PES obtained through the active learning procedure, we
ran long metadynamics trajectories and reconstructed their free energy surface
to determine the transition state free energy for reactions 9, 13, 14, and 18.
Unlike in the exploration phase where any CV can help discover new structures,
in the real application to reconstruct the free energy surface, a good choice
of CVs is important. Due to the diffusive nature of gas phase reactions, the
distance between two fragments can be arbitrary and lead to convergence
problem if we directly use distance as the collective variable. Here we use
CVs inspired by our recent study of diffusion maps in order to assess CVs
suitable for each reaction as evaluated through their correlations with the
diffusion coordinates. We found that the coordination number(CN) between atom
i and j:

$CN(ij)=\frac{2}{1+exp(3(r_{ij}-r_{eq}))}$ (1)

is often a good CV, where $r_{eq}$ is the equilibrium distance between i and
j, usually selected as the equilibrium bond length unless otherwise specified.
Restraints on certain bond distance and bond angles are also necessary to keep
the metadynamics trajectory inside the relevant chemical space defined by the
CVs.

We also exploit the strength of having an ensemble of ML models in which the
standard deviation between the models measures the reliability in the
prediction of ML forces on a given configuration. We therefore devise a hybrid
mode of running a simulation in which we use ML forces when the standard
deviation among models is relatively small, and with ab initio forces at the
same level of theory used for ML training when the model predictions are not
as reliable (Figure 4).

Figure 4: Schematic illustration of the new workflow for rebuilding the free
energy surface. The metadynamics proceeds with a model that utilizes query by
committee ML trajectories to signal when forces have degraded. In lieu of
retraining, the AIMD forces are directly substituted to complete the MD time
step.

As we support below, this hybrid mode can be (1) much less expensive than
fully retraining the ML models with the new DFT data, and (2) will be more
accurate and maintain stability when uncertainties arise, although with
appropriate parallelization one could imagine a constant retraining while the
MD trajectory proceeds with the hybrid model. Ultimately this involves a
tradeoff between how complete is the active learning process vs. The relative
cost of the ML and DFT forces; in our case on average it takes $\sim$0.11
seconds for a ML force step and $\sim$12 seconds to generate the ab initio
forces.

Figure 5 shows the reconstructed free energy surfaces using the hybrid model
for the more interesting reaction channels where the reactant-to-product
ratios are uneven. In Rxn 9, when we use CN(O1-H3) and CN(O2-H3) as the
collective variables we find that the IRC transition state leans towards the
reactant side on the free energy surface, consistent with the 62:38 committer
ratio, and thus shifts quite significantly on the free energy surface toward
the product well (Figure 5a). For Rxn 10 we found CN(O2-O3) and CN(O2-H4) to
be good collective variables, where CN(O2-O3) monitors the reaction progress
and CN(O2-H4) helps separates the transition state from the stable wells.
Consistent with the 55:45 reactant-to-product committer ratio the IRC
transition state leans towards the reactant side, but on the free surface two
possible transition pathways are evident (Figure 5b). One pathway is very
similar to the IRC but with the transition state moving closer to the product
well, while the other pathway demarcated by a smaller CN(O2-H4) value involves
breaking the O2-H3 and forming the H3-H4 molecule in an asynchronous manner.

Figure 5: Free energy surface reconstructed from metadynamics for hydrogen
combustion. (a) Reaction 9, (b) Reaction 10, (c) Reaction 14, and (d) Reaction
1. The original IRC pathway is labeled with a red dot for the transition state
and the free energy transition state is labeled as a purple star.

The IRC transition state of both Rxn14 and Rxn01 lean very heavily toward the
product side on the free energy surface (Figure 5c,d), as is consistent with
the much more skewed commitor reactant-product ratios in Table 1. Accordingly,
the free energy transition state shifts dramatically in both reactions, with
saddle point regions that exhibit very large free energy stabilization of 8-12
kcal/mole relative to the IRC transition state arising from entropic effects.
As a control we tested Rxn 17 and Rxn 18 in which the reactant-to-product
ratio in the committer analysis is close to 50-50. In both cases the IRC and
free energy transition state are consistent indicating that the coordination
numbers tend to be good collective variables and that the entropic factors
governing these reactions are small.

Through this process, we find that the number of unreliable region of PES are
relatively infrequent and the deviations are small enough such that energy
conservation is not badly effected. For the Rxn 18 trajectory the full DFT
force was called 1364 times over a total of 1 million steps, which means
99.86% of the trajectory is driven by the ML model, and only 0.14% of steps is
calculated through the ab initio update. For the Rxn10 trajectory the DFT
force was called only 136 times, which means that 99.99% of the total MD steps
were generated by the ML forces. The small chance to go through ab initio
updates reflects that the learned PES after our negative design and
metadynamics active learning procedures is largely complete, especially for
Rxn 10. But it is important to emphasize that without the ab initio calls the
MD trajectories would become corrupted beyond repair and this hybrid ensemble
approach allows for highly efficient simulation with a time scale similar to
the pure ML run but with more stability, offering a preventative and
alternative to occasional model hallucinations.

##  3 DISCUSSION AND CONCLUSIONS

Machine learning (ML) methods have traditionally aimed to train a faithful
surrogate model of ab initio energy and forces for different configurational
arrangements of a system in order to explore reaction chemistry more
efficiently than ab initio molecular dynamics. For real applications it is
likely that the interpolative nature of machine learning will always suffer
from data insufficiencies that will limit its total replacement of physics
based methods. However, this work emphasizes that we can exploit multiple
hybrid model strategies to accelerate chemical predictions that will be
reliable, especially in regards chemical reactivity. The first recognition is
that we need a negative design aspect to data acquisition. Although there has
been some focus on collecting data on possible intermediates or in the
metastable regions of the potential energy surface, in fact very unphysical
states also should be part of the data. But nearly all reference datasets are
insufficient as they lack unusual and/or high energy configurations because
they are not considered to be chemically relevant.

To uncover this more unintuitive data, we devise an active learning workflow
where we use enhanced sampling within a lower manifold of collective
coordinates along with a query by committee stategy to systematically improve
the ML PES in its uncertain region. By conscious negative design through the
dialation approach and more diverse sampling with metadynamics, we achieved a
more complete and robust ML representation of the hydrogen combustion reactive
PES. Metadynamics, as a rare event sampling method, turns out to be an
efficient sampling tool for unstable structures, which helps the AL workflow
to identify holes in the PES landscape to inform the ML model through
retraining with such data. Using metadynamics only as a sampling tool, the
tricky CV selection step can be avoided by starting with any reasonable or
intuitive CVs.

Ultimately to address real world application studies using ML potentials, any
a priori formulated reference dataset, even with active learning approaches,
will be too small-scale for true PES completeness. In this case any large
variance among ML models provides a nice strategy to invoke a hybrid model,
which smooths over rough spots on the ML PES with calls to the ab initio data
source. This hybrid ML-AIMD model allowed us to reconstruct and perform
analyses on free energy surfaces for many of the reaction channels for
hydrogen combustion with much greater accuracies, but at hoped for ML
computational efficiencies reaching two orders of magnitude improvement over
AIMD.

One alternative approach could be delta learning that learns the difference
between a lower level theory and a higher level theory. This type of
simulation would require calling both the lower theory method and the ML model
at each and every step, and combine them to mimic the PES of the higher level
theory. This approach can avoid the detrimental effect from ML rough spots by
making ML correction zero whenever the ML models disagree. However, the
quality and cost of this type of simulation would highly depend on the
underlying lower level theory method, and switching between the two surface
can affect the energy conservation in the trajectory. We have explored the
idea of using a physical model (C-GeM 26) to represent the long range behavior
and use ML to delta-learn the rest of the interaction, but found that the
intrinsic error from the C-GeM model adds challenge for ML models to learn the
more complex pattern. At such, we give up on this delta learning idea and
decide to stick with the hybrid model we devise.

##  4 METHODS

###  4.1 Dialation Data Preparation

To address the problem that the initial Dataset from Ref. 17 is missing high
energy states, we prepared some additional data by proportionally scaling each
geometry in the IRC is with multiple ratio (0.6, 0.7, 0.8, 0.9, 1.1, 1.2, 1.3,
1.4 , 1.6, 1.8, 2.0, 2.4, 2.8, 3.2). The energy and forces of these new
geometries were obtained with QChem27 force calculation using range separated
hybrid meta-GGA functional $\omega$B97X-V 18 with the cc-pVTZ basis set.

###  4.2 Active Learning: Query by Committee

In the field of molecular modeling, active learning has found another valuable
application for reducing ML errors through improving data selectivity and
efficiency. Here we exploit the query by committee idea to find a label the
most informative data points through this iterative process:

  1. 1.

Perform short metadynamic simulations to explore the configuration space in a
lower dimension.

  2. 2.

When the four models disagree outside standard deviation, collect a
representative subset of structures to be included in the training set through
downsampling.

  3. 3.

Perform DFT calculation of energies and atomic forces to label the new data.

  4. 4.

Retrain the ensemble of ML models with the updated training set.

The details of each of the steps are described in the following subsections.

###  4.3 Training of the NN PES

The NewtonNet model was used to train the NN potential in this work. Details
of this method can be found in Ref. 16. The initial four models are trained
for 2000 epochs with 1000 data points per reaction from the original dataset
with addtional 200 data points from dialation data, initialized with different
weights. The following models are trained with exactly the same setting but
with additional data sampled from short metadynamics simulations. The cutoff
radius was set to 5.0 Å. The initial learning rate was set to 0.001 and with
0.7 learning rate decay. The loss is defined as

$\begin{split}\mathcal{L}&=\frac{\lambda_{E}}{M}\sum_{m}^{M}w_{m}\left(\widetilde{E}_{m}-E_{m}\right)^{2}\\\
&+\frac{\lambda_{F}}{M}\sum_{m}^{M}\frac{w_{m}}{3N_{m}}\sum_{i}^{N_{m}}\left\|\widetilde{\bm{F}}_{mi}-\bm{F}_{mi}\right\|^{2}\\\
&+\frac{\lambda_{D}}{M\times
N_{m}}\sum_{m}^{M}w_{m}\sum_{i}^{N_{m}}\left(1-\frac{{\bm{\mathcal{F}}}_{mi}^{L}\cdot\bm{F}_{mi}}{\left\|{\bm{\mathcal{F}}}_{mi}^{L}\right\|\left\|\bm{F}_{mi}\right\|}\right)\end{split}$

where the prefactor of the energy error $\lambda_{E}$ is set to 1, the
prefactor of force error $\lambda_{F}$ is set to 20, and the prefactor
$\lambda_{D}$ for latent force direction is set to 1.

We also use the following Boltzmann weighting factor $w_{m}$, defined as

$w_{m}=\begin{cases}1&\text{ if }E_{m}<=E_{thresh}\\\
exp({\frac{(-(E_{m}-E_{thresh})}{k_{B}T})})&\text{ if
}E_{m}>E_{thresh},\par\end{cases}$

to bias the training towards data point within a relevant energy scale.
$E_{thresh}$ is a per-atom quantity that put less weighting on all data points
with energy higher than this number. It is chosen to be 16.744 kcal/mol/atom
for the Hydrogen combustion system, which is 10 kcal/mol/atom higher than the
highest per atom energy among all reaction channels in IRCs. To completely
converge the ML model, we added one final training step with all previously
added data, using a larger epoch size (5000 steps), and a more patient
learning rate decay to give the final model that we later use for determining
predictions on the free energy surface.

###  4.4 Metadynamics

New structures are sampled through short metadynamics trajectories for 6
reaction channels: rxn09, rxn10, rxn13, rxn16, rxn17 and rxn18. For each step
in the metadynamics simulation, the atomic forces are evaluated by four NN PES
models simultaneously. Outlier among the 4 predictions is removed if the
absolute difference between the outlier in question and the closest number is
larger than the 95% confidence limit value of the Dixon Q’s test. Then the
mean of model predicted forces is modified by plumed 28 to allow for enhanced
sampling. In this work, all enhanced sampling simulations are performed with
the well-tempered metadynamics, in which a Gaussian centered at the visited
point is periodically added to the potential. The simulation is driven through
the atomic simulation environment(ASE)29 with a specifically tailored
calculator that provide energy and forces for a given structure through the
above protocol. The simulation is conducted at 300K with increasing length as
we have more active learning rounds: 2ps between active learning round 1-20,
5ps round 21-33 and 10ps between round 34-48.

During the simulation, standard deviation on atomic forces over an ensemble of
NN potentials is monitored. Whenever the maximum over the atoms exceeds a
predefined threshold (2 kcal/mol/Å), the configuration is selected for further
downsampling.

###  4.5 Clustering and Down Selection

Each molecule is first represented as a Coulomb matrix30 that includes the
nuclear charges (Zi and Zj) of atom i and j along with their Cartesian
coordinates (Ri and Rj).

$\bm{C}_{ij}=\begin{cases}0.5Z_{i}^{2.4},i=j\\\
\frac{Z_{i}Z_{j}}{\lvert\bm{R_{i}}-\bm{R_{j}}\rvert},i\neq j\par\end{cases}$
(2)

To reduce the dimension of the dataset while retaining the majority of
structural information, the Coulomb matrix was transformed into the eigen-
spectrum by solving the eigenvalue problem $\bm{Cv}=\lambda\bm{v}$ subject to
the constraint $\lambda_{i}\geq\lambda_{i+1}$. The Mini Batch KMeans
clustering algorithm was then applied to categorize the sub-datasets into
smaller clusters based on the eigen-spectrum. The value of K (the number of
clusters) is chosen automatically with a scaled inertia approach 31. The
scaled inertia is formulated as

$\text{Scaled Inertia}=\frac{I(K)}{I(K=1)}+\alpha K$ (3)

where the inertia(I) is the sum of squared distance of samples to their
closest cluster center:

$I(K)=\sum_{i=1}^{N}{(x_{i}-C_{k})^{2}}$ (4)

where N is the number of samples and $C_{k}$ is the centroid of a cluster.
$\alpha$ is a manually tuned factor that gives penalty to the number of
clusters, here we chose $\alpha=0.0002$. We chose the K value that gives the
minimum scaled inertia among all K<300 to do the mini batch K-means clustering
on all molecules from a given reaction channel. Afterward, we randomly picked
a structure from each cluster, whose energy and forces will be calculated and
then be included into the new training set. The final data set had 48,582 data
points sampled in the AL procedure.

###  4.6 DFT single point force calculations

The structures that passes the clustering and downselection process are
gathered for labeling and retraining. They are labeled through DFT force
calculations through QChem using $\omega$B97X-V 18 functional with the cc-pVTZ
basis set. Some of the structures are tricky to deal with due to the diffusive
nature of gas phase simulation. Therefore, for all structures with atomic
oxygen that is separated from the rest of the structure for more than 2
Åapart, we introduced and additional run with FRAGMO initial guess for SCF
calculation. 32 We compared the final energy obtained from a calculation with
fragmo initial guess and the one with the usual superposition of atomic
densities(SAD) inital guess, and take whichever one with lower energy to give
the final energy and forces for labeling the data.

###  4.7 Reconstructing free energy surface with metadynamics

The free energy surfaces were calculated from longer metadynamic simulations.
The Langevin thermostat was used to maintain temperature at 300K, with a
friction coefficient of 0.002 a.u. In the metadynamics simulations, the
Gaussians adopted have an initial height of 5 kJ/mol and width of 0.05 for the
CVs. A Gaussian was deposited every 100 step with a bias factor equal to 10.
The simulations were 200 ps long with step size of 0.2 fs. The free energy
were calculated using the sum_hills utility in Plumed28 and its surface
plotted with python matplotlib.

##  5 ACKNOWLEDGMENTS

We thank the CPIMS program, Office of Science, Office of Basic Energy
Sciences, Chemical Sciences Division of the U.S. Department of Energy under
Contract DE-AC02-05CH11231 for support of the machine learning approach to
hydrogen combustion. This work was supported by the U.S. Department of Energy,
Office of Science, Office of Advanced Scientific Computing, and Office of
Basic Energy Sciences, via the Scientific Discovery through Advanced Computing
(SciDAC) program for the collective variables. This work used computational
resources provided by the National Energy Research Scientific Computing Center
(NERSC), a U.S. Department of Energy Office of Science User Facility operated
under Contract DE-AC02-05CH11231.

##  6 DATA AND CODE AVAILABILITY

Coordinates of geometries used in the training:
https://doi.org/10.6084/m9.figshare.23290115.v1. The full workflow code can be
found in https://github.com/THGLab/ H2Combustion_AL.

##  7 SUPPORTING INFORMATION

We provide a table describing the hydrogen combustion reaction channels 1-19,
and the amount of data accumulated over active learning rounds.

##  8 DECLARATION OF INTERESTS

The authors declare no competing interests.

##  9 AUTHOR CONTRIBUTIONS

X.G. and T.H.G. designed the project. X.G. carried out the AIMD simulations,
metadynamic calculations, and active learning. X.G. and T.H.G designed the
collective coordinates with the help of J.H.,T.K.,C.Y. All authors discussed
the results and made comments and edits to the manuscript.

## References

  * 1 Behler, J. & Parrinello, M.  Generalized neural-network representation of high-dimensional potential-energy surfaces.  _Physical Review Letters_ 98, 146401 (2007).  URL https://link.aps.org/doi/10.1103/PhysRevLett.98.146401.
  * 2 Smith, J. S., Isayev, O. & Roitberg, A. E.  Ani-1, a data set of 20 million calculated off-equilibrium conformations for organic molecules.  _Scientific Data_ 4, 170193 (2017).  URL https://doi.org/10.1038/sdata.2017.193.
  * 3 Smith, J. S., Nebgen, B., Lubbers, N., Isayev, O. & Roitberg, A. E.  Less is more: Sampling chemical space with active learning.  _The Journal of Chemical Physics_ 148, 241733 (2018).  URL https://doi.org/10.1063/1.5023802.
  * 4 Smith, J. S. _et al._ Approaching coupled cluster accuracy with a general-purpose neural network potential through transfer learning.  _Nature Communications_ 10, 2903 (2019).  URL https://doi.org/10.1038/s41467-019-10827-4.
  * 5 Chmiela, S., Sauceda, H. E., Müller, K.-R. & Tkatchenko, A.  Towards exact molecular dynamics simulations with machine-learned force fields.  _Nature Communications_ 9, 3887 (2018).  URL https://doi.org/10.1038/s41467-018-06169-2.
  * 6 Chmiela, S. _et al._ Machine learning of accurate energy-conserving molecular force fields.  _Science Advances_ 3, e1603015 (2017).  URL http://advances.sciencemag.org/content/3/5/e1603015.abstract.
  * 7 Schutt, K. T., Sauceda, H. E., Kindermans, P. J., Tkatchenko, A. & Müller, K. R.  Schnet - a deep learning architecture for molecules and materials.  _Journal of Chemical Physics_ 148 (2018).
  * 8 Thomas, N. _et al._ Tensor field networks: Rotation- and translation-equivariant neural networks for 3d point clouds.  _arXiv preprint arXiv:1802.08219_ (2018).
  * 9 Drautz, R.  Atomic cluster expansion for accurate and transferable interatomic potentials.  _Phys. Rev. B_ 99, 014104 (2019).  URL https://link.aps.org/doi/10.1103/PhysRevB.99.014104.
  * 10 Anderson, B., Hy, T. S. & Kondor, R.  Cormorant: Covariant molecular neural networks.  _Adv. Neural Inf. Process. Syst._ 32 (2019).
  * 11 Qiao, Z., Welborn, M., Anandkumar, A., Manby, F. R. & Miller, I., Thomas F.  Orbnet: Deep learning for quantum chemistry using symmetry-adapted atomic-orbital features.  _The Journal of Chemical Physics_ 153, 124111 (2020).  URL https://doi.org/10.1063/5.0021955.
  * 12 Batzner, S. _et al._ Se(3)-equivariant graph neural networks for data-efficient and accurate interatomic potentials.  _arXiv preprint arXiv:2101.03164_ (2021).
  * 13 Glick, Z. L., Koutsoukas, A., Cheney, D. L. & Sherrill, C. D.  Cartesian message passing neural networks for directional properties: Fast and transferable atomic multipoles.  _The Journal of Chemical Physics_ 154, 224103 (2021).
  * 14 Schütt, K., Unke, O. & Gastegger, M.  Equivariant message passing for the prediction of tensorial properties and molecular spectra.  In Meila, M. & Zhang, T. (eds.) _Proceedings of the 38th International Conference on Machine Learning_ , vol. 139 of _Proceedings of Machine Learning Research_ , 9377–9388 (PMLR, 2021).  URL https://proceedings.mlr.press/v139/schutt21a.html.
  * 15 Unke, O. T. _et al._ Spookynet: Learning force fields with electronic degrees of freedom and nonlocal effects.  _Nature Communications_ 12, 7273 (2021).  URL https://doi.org/10.1038/s41467-021-27504-0.
  * 16 Haghighatlari, M. _et al._ Newtonnet: a newtonian message passing network for deep learning of interatomic potentials and forces.  _Digital Discovery_ 1, 333–343 (2022).
  * 17 Guan, X. _et al._ A benchmark dataset for hydrogen combustion.  _Scientific Data_ 9, 215 (2022).
  * 18 Mardirossian, N. & Head-Gordon, M.  $\omega$B97X-V: A 10-parameter, range-separated hybrid, generalized gradient approximation density functional with nonlocal correlation, designed by a survival-of-the-fittest strategy.  _Phys. Chem. Chem. Phys._ 16, 9904–9924 (2014).
  * 19 Kulichenko, M. _et al._ Uncertainty-driven dynamics for active learning of interatomic potentials.  _Nature Computational Science_ 3, 230–239 (2023).
  * 20 Shapeev, A., Gubaev, K., Tsymbalov, E. & Podryabinkin, E.  Active learning and uncertainty estimation (1970).  URL https://link.springer.com/chapter/10.1007/978-3-030-40245-7_15#Sec8.
  * 21 Laio, A. & Parrinello, M.  Escaping free-energy minima.  _Proceedings of the National Academy of Sciences_ 99, 12562–12566 (2002).  URL https://doi.org/10.1073/pnas.202427399.
  * 22 Barducci, A., Bussi, G. & Parrinello, M.  Well-tempered metadynamics: A smoothly converging and tunable free-energy method.  _Phys. Rev. Lett._ 100, 020603 (2008).  URL https://link.aps.org/doi/10.1103/PhysRevLett.100.020603.
  * 23 Ko, T. _et al._ Using diffusion maps to analyze reaction dynamics for a hydrogen combustion benchmark dataset (2023).  2304.09296.
  * 24 Bertels, L. W., Newcomb, L. B., Alaghemandi, M., Green, J. R. & Head-Gordon, M.  Benchmarking the performance of the reaxff reactive force field on hydrogen combustion systems.  _The Journal of Physical Chemistry A_ 124, 5631–5645 (2020).
  * 25 Li, J., Zhao, Z., Kazakov, A. & Dryer, F. L.  An updated comprehensive kinetic model of hydrogen combustion.  _Int. J. Chem. Kinet._ 36, 566–575 (2004).
  * 26 Guan, X., Leven, I., Heidar-Zadeh, F. & Head-Gordon, T.  Protein c-gem: A coarse-grained electron model for fast and accurate protein electrostatics prediction.  _Journal of Chemical Information and Modeling_ 61, 4357–4369 (2021).
  * 27 Shao, Y. _et al._ Advances in molecular quantum chemistry contained in the q-chem 4 program package.  _Molecular Physics_ 113, 184–215 (2014).
  * 28 Tribello, G. A., Bonomi, M., Branduardi, D., Camilloni, C. & Bussi, G.  Plumed 2: New feathers for an old bird.  _Computer Physics Communications_ 185, 604–613 (2014).
  * 29 Larsen, A. H. _et al._ The atomic simulation environment—a python library for working with atoms.  _Journal of Physics: Condensed Matter_ 29, 273002 (2017).
  * 30 Rupp, M., Tkatchenko, A., Müller, K.-R. & von Lilienfeld, O. A.  Fast and accurate modeling of molecular atomization energies with machine learning.  _Physical Review Letters_ 108 (2012).
  * 31 Herman-Saffar, O.  An approach for choosing number of clusters for k-means (2021).  URL https://towardsdatascience.com/an-approach-for-choosing-number-of-clusters-for-k-means-c28e614ecb2c.
  * 32 Khaliullin, R. Z., Cobar, E. A., Lochan, R. C., Bell, A. T. & Head-Gordon, M.  Unravelling the origin of intermolecular interactions using absolutely localized molecular orbitals.  _The Journal of Physical Chemistry A_ 111, 8753–8765 (2007).


# Non-classical nucleation of zinc oxide from a physically-motivated machine-
learning approach.

Gaétan Laurens  Institut Lumière Matière, UMR5306 Université Lyon 1-CNRS,
Université de Lyon, 69622 Villeurbanne Cedex, France  Jacek Goniakowski  CNRS,
Sorbonne Université, Institut des NanoSciences de Paris, UMR 7588, 4 Place
Jussieu, F-75005 Paris, France  Julien Lam<EMAIL_ADDRESS>Centre
d’élaboration des Matériaux et d’Etudes Structurales, CNRS (UPR 8011), 29 rue
Jeanne Marvig, 31055 Toulouse Cedex 4, France

###### Abstract

Observing non-classical nucleation pathways remains challenging in simulations
of complex materials with technological interests. This is because it requires
very accurate force fields that can capture the whole complexity of their
underlying interatomic interactions and an advanced structural analysis. Here,
we first report the construction of a machine-learning force field for zinc
oxide interactions using the Physical LassoLars Interaction Potentials
approach which allows us to be predictive even for untrained structures. Then,
we carried out freezing simulations from a liquid and observed the crystal
formation with atomistic precision. Our results, which are analyzed using a
data-driven approach based on bond order parameters, demonstrate the presence
of both prenucleation clusters and two-step nucleation scenarios thus
retrieving seminal predictions of non-classical nucleation pathways made on
much simpler models.

##  I Introduction

While crystals in material science are ubiquitous, the mechanisms of their
formation which spans from nucleation to crystal growth remain one of the most
intriguing process in nature. Better understanding crystallization would allow
for a rational control of material engineering and possibly the development of
novel functional materials and technological applications. From the
fundamental point of view, numerous works have been dedicated to elucidating
the emergence of the nucleation core Lutsko and Lam (2020); Lutsko (2019);
Wolde and Frenkel (1997); Auer and Frenkel (2001); Ten Wolde  _et al._ (1995)
and its role in controlling the final crystal structure Sleutel and Van
Driessche (2014); Van Driessche  _et al._ (2018); Menon  _et al._ (2020); Lam
and Lutsko (2018); Amodeo  _et al._ (2020); Desgranges and Delhommelle (2007a,
2006). For instance, it is now possible to observe the crystal birth with
electron microscopy Sleutel and Van Driessche (2014); Van Driessche  _et al._
(2018); Yamazaki  _et al._ (2017); Sleutel  _et al._ (2014); Ogata  _et al._
(2020), and colloidal science has also provided numerous experimental results
on nucleation van Blaaderen  _et al._ (1997); Leunissen  _et al._ (2005);
Velikov  _et al._ (2002); de Nijs  _et al._ (2015). Yet, numerical simulations
remain the principle instrument to investigate crystallization at the
atomistic level Sosso  _et al._ (2016). In this context, most of these works
require large scale simulations in order to observe the phase transition and
have therefore only focused on materials for which the interactions are simple
including hard-spheres Auer and Frenkel (2001); Coli and Dijkstra (2021);
Russo and Tanaka (2012); Leoni and Russo (2021), Lennard-Jones Lutsko (2019);
Ten Wolde  _et al._ (1995); Lam and Lutsko (2018); Desgranges and Delhommelle
(2007a), water Leoni and Russo (2021); Lin  _et al._ (2018); Metya and
Molinero (2021); Qiu and Molinero (2018); Qiu  _et al._ (2019); Lupi  _et al._
(2014); Moore and Molinero (2011); Lupi  _et al._ (2017); Russo  _et al._
(2014), as well as metallic potentials like embedded-atom model (EAM) Menon
_et al._ (2020); Amodeo  _et al._ (2020); Desgranges and Delhommelle (2007b).

Prompted by this large body of fundamental achievements, it becomes timely to
reach the same level of understanding for crystallization with more complex
materials in order to target more diverse technological applications. So far,
the need for large scale simulations have prevented from using quantum
accurate modeling including density functional theory (DFT), and the research
field dedicated to constructing novel interaction potentials to bridge this
computational gap has been ever expanding Combettes  _et al._ (2020); Laurens
_et al._ (2020); Calvo  _et al._ (2017); Zipoli and Curioni (2013);
Rajasekaran  _et al._ (2016); Senftle  _et al._ (2016); Brommer  _et al._
(2015); Vashishta  _et al._ (2008). In particular, the past decade has seen
the emergence of innovative types of interaction potentials based on machine-
learning algorithms Behler and Parrinello (2007). Various approaches have been
proposed such as Artificial Neural Networks Behler and Parrinello (2007),
Gaussian approximation potentials Bartók  _et al._ (2010), Linearized
potentials Seko  _et al._ (2015, 2014); Takahashi  _et al._ (2017); Seko  _et
al._ (2019); Goryaeva  _et al._ (2019); Benoit  _et al._ (2020), Spectral
Neighbor Analysis Potential Thompson  _et al._ (2015); Wood and Thompson
(2018), Symmetric Gradient Domain Machine learning Chmiela  _et al._ (2017,
2018), and Moment Tensor Potentials Shapeev (2016); Novikov  _et al._ (2020).
The success of those machine-learning interaction potentials (MLIPs) is seen
through the large variety of studied materials, namely pure metals Novoselov
_et al._ (2019); Seko  _et al._ (2015); Takahashi  _et al._ (2018); Zeni  _et
al._ (2018); Botu  _et al._ (2017), organic molecules Bereau  _et al._ (2018);
Sauceda  _et al._ (2019); Bartók  _et al._ (2017); Veit  _et al._ (2019),
water Nguyen  _et al._ (2018); Bartók  _et al._ (2013); Morawietz  _et al._
(2012); Natarajan and Behler (2016); Morawietz and Behler (2013), amorphous
materials Deringer and Csányi (2017); Bartók  _et al._ (2018); Caro  _et al._
(2018); Deringer  _et al._ (2018a, b); Sosso  _et al._ (2018), and hybrid
perovskites Jinnouchi  _et al._ (2019).

The oxide materials has been much less investigated using MLIPs Sivaraman  _et
al._ (2020); Artrith  _et al._ (2013); Eckhoff  _et al._ (2020a); Artrith  _et
al._ (2011); Sundararaman  _et al._ (2018); Eckhoff  _et al._ (2020b); Artrith
_et al._ (2017); Artrith and Urban (2016). With those materials, the
complexity lies in the emergence of long-range electrostatic effects difficult
to capture with traditional machine-learning approaches Ko  _et al._ (2021a,
b). In addition, oxide materials often exhibit a rich structural landscape
composed of numerous polymorphs Navrotsky  _et al._ (2008); Machala  _et al._
(2011); Sponza  _et al._ (2015); Tsybulya and Kryukova (2008). Being able to
model those polymorphs with a unique set of atomic interaction remains
challenging even for advanced machine-learning methods.

In this article, we focus on the modeling of zinc oxide (ZnO) using our newly
developed approach named Physical LassoLars Interaction Potential (PLIP) which
is based on a physically-motivated mathematical formulation of the potential
and on a constrained linear regression scheme to adjust its parameters Benoit
_et al._ (2020). Beyond the artificial neural network potentials taking into
account electrostatic effects Artrith  _et al._ (2011, 2013), since we focus
on modeling ZnO crystallization, we have constructed potential able to account
for six most stable bulk ZnO polymorphs and their low index surfaces, as well
as for ZnO liquid and amorphous phases.

In particular, we will first describe our numerical method and demonstrate the
ability of our potential to model zinc oxide in various situations. Then, we
will present the results of molecular dynamics (MD) simulations using our PLIP
to study the freezing of a bulk liquid mixture of zinc and oxygen atoms. Our
results show that while traditional interaction potentials of zinc oxide may
allow for the formation of meta-stable crystals, the use of a more precise
machine-learning interaction potential leads to the emergence of only the most
stable crystal structure. Finally, we also employ data-driven approaches to
distinguish between those different polymorphs and study at the atomistic
scale the birth of the zinc oxide crystal. Our results indicate that numerous
nucleation mechanisms are in competition including a two-step nucleation
process and the emergence of core-shell crystals within the out-of-equilibrium
liquid. From a material point of view, since zinc oxide crystals are found in
numerous applications, $e.g.$ photocatalysisMirzaeifard  _et al._ (2020),
piezzoelectricityMahmood  _et al._ (2021), drug deliveryXiong (2013),
antibacterialAbebe  _et al._ (2020), and gas detection Yuan  _et al._ (2019),
our results provide new insights into the material processing.

##  II Methods

###  II.1 Machine-learning interaction potential using the PLIP approach

As in most MLIP approaches, the total potential energy $E$ is decomposed as
the sum of independent atomic energies $E_{i}$: $E=\sum^{N}_{i=0}E_{i}$, where
$N$ corresponds to the number of atoms of the considered configuration. Then,
a linear model consists in approximating $E_{i}$ as a linear combination of
descriptors $X_{n}^{i}$: $E_{i}=\sum_{n}\omega_{n}X_{n}^{i}$ where
$\omega_{n}$ are the linear coefficients that must be determined. More
specifically, our PLIP model employs three types of descriptors which
explicitly follow a many-body order expansion:

$\displaystyle[2B]_{n}^{i}=\sum_{j}f_{n}(R_{ij})\times f_{c}(R_{ij}),$ (1)
$\displaystyle[3B]_{n,l}^{i}=\sum_{j}\sum_{k}f_{n}(R_{ij})f_{c}(R_{ij})f_{n}(R_{ik})f_{c}(R_{ik})cos^{l}(\theta_{ijk}),$
(2) $\displaystyle[NB]_{n,m}^{i}=\left(\sum_{j}f_{n}(R_{ij})\times
f_{c}(R_{ij})\times f_{s}(R_{ij})\right)^{m},$ (3)

where $R_{ij}$ is the distance between atoms $i$ and $j$, $\theta_{ijk}$ is
the angle centered around the atom $i$, and $l$ and $m$ are two positive
integers. $f_{c}(R_{ij})$ is the same cut-off function as in the seminal work
of Behler and Parinello Behler and Parrinello (2007), namely
$f_{c}(R_{ij})=0.5\left(1+cos(\pi(R_{ij}/R_{cut}))\right)$, where the cut-off
distance $R_{cut}$ is set at 6 Å. The shift function is chosen in the from:
$f_{s}(u)=6u^{5}-15u^{4}+10u^{3}$, where $u=(R_{ij}-r_{1})/(r_{2}-r_{1})$, and
$r_{1}$ (resp. $r_{2}$) is defined as 95 % (resp. 105 %) of a short distance
equal to $1.1$ Å. Regarding the basis functions $f_{n}(R_{ij})$, we previously
demonstrated the advantage of coupling different classical functions at the
same time Benoit  _et al._ (2020). Yet, in this particular study on zinc oxide
where numerical efficiency is crucial, we work only with Gaussian functions
for which the width and the central positions are respectively listed as
follows: [0.5, 1.0, 1.5] and [0.5, 1.0, 1.5 … 5.5, 6.0]. In addition, we vary
the integers $l$ and $m$ respectively from $0$ to $5$ Å-1, and from $4$ to $7$
Å. Altogether, when taking into account the binary nature of the zinc oxide
system, our model is made of $1981$ available descriptors.

In order to match the first-principle database, the model is fitted using the
LassoLars approach which allows for a well-informed selection of the most
preponderant descriptors and a reduction in the complexity of the obtained
potential Benoit  _et al._ (2020); Zeni  _et al._ (2021). In this study, we
built three PLIPs using for each approximately $150$ descriptors selected
among an order of magnitude more descriptors.

###  II.2 First principle calculations and training database

Reference GGA-DFT calculations are performed with VASP Kresse and Hafner
(1993); Kresse and Furthmuller (1996), using PW91 exchange-correlation
functional Perdew  _et al._ (1996), and the projector augmented wave method
Kresse and Joubert (1999). Standard zinc and soft oxygen (energy cutoff of 270
eV) pseudopotentials provided by VASP are used in all calculations, enabling
an efficient structural relaxation of systems composed of several hundreds of
atoms. GGA results obtained with the soft and the full (energy cutoff of 400
eV) oxygen pseudopotentials show satisfactory agreement (differences of bulk
parameters $a$ and $c$ smaller than 0.01 Å, and cohesion energy differences
$\Delta E_{\mathrm{coh}}$ below 0.01 eV/ZnO). The present results on six ZnO
polymorphs [wurtzite (WRZ), zinc blend (ZBL), body centered tetragonal (BCT),
sodalite (SOD), $h$-BN (HBN), and cubane (CUB) crystallographic structures]
coincide very well with the existing hybrid HSE03 reference Sponza  _et al._
(2015) [See Tab. 1], in particular concerning the relative stability of the
considered polymorphs. Such agreement between the two types of modeling
indicates that the chosen GGA-DFT approach leads to sufficiently accurate
results while enabling for larger scale calculations.

Table 1: Calculated (present PW91 and reference HSE03 Sponza  _et al._ (2015)) characteristics of the six bulk ZnO polymorphs: lattice parameters $a$ and $c$ (Å), and cohesion energy differences with respect to ground state (WRZ) structure $\Delta E_{\mathrm{coh}}$ (eV/ZnO). | PW91 | HSE03
---|---|---
|  a, c |  $\Delta E_{\mathrm{coh}}$ |  a, c |  $\Delta E_{\mathrm{coh}}$
WRZ |  3.28, 5.30 |  -0.00 |  3.25, 5.24 |  -0.00
ZBL |  3.27, 3.27 |  -0.01 |  3.24, 3.24 |  -0.03
BCT |  5.61, 3.29 |  -0.05 |  5.55, 3.25 |  -0.05
SOD |  5.67, 5.67 |  -0.14 |  5.61, 5.61 |  -0.14
HBN |  3.46, 4.59 |  -0.14 |  3.43, 4.50 |  -0.14
CUB |  6.28, 6.28 |  -0.23 |  6.22, 6.22 |  -0.24

Since the training database includes both ordered and disordered structures,
in DFT calculations we systematically use relatively large supercells (16-19 Å
large cuboids), containing 320-480 atoms, and sample the Brillouin zone with a
single $\Gamma$ point. In ordered structures, the atomic coordinates of all
ions are relaxed until residual forces dropped below 0.01 eV/Å and, in the
bulk structures, all components of the stress tensor are smaller than 0.01
eV/Å3. The disordered structures are collected along short high-temperature MD
runs with no further relaxation as to produce non-vanishing forces on all
ions. For those first MD runs, we use a classical Buckingham ZnO
potentialBinks and Grimes (1993) which was previously employed in numerous
studies of ZnOWang  _et al._ (2010); Lin  _et al._ (2013); Kulkarni  _et al._
(2005, 2006). We would like to emphasize that using a classical potential
instead of first-principle calculations at this stage favors a more rapid
sampling of a larger variety of configurations.

In order to include information on the characteristics of low coordinated
ions, for each of the six ZnO polymorphs we also consider at least one low
index non-polar surface, whereas three surface orientations are considered for
the most stable (WRZ) polymorph. In this latter case, beyond the most stable
(10-10) and (11-20) non-polar orientations, we also include the polar
(0001)/(000-1) surfaces. As to impede the emergence of a macroscopic dipole
moment, we use an asymmetric slab with the conventional (2$\times$2)
reconstructed surfaces at which one surface oxygen (zinc) ion of each four is
removed at the oxygen (zinc) termination Goniakowski  _et al._ (2007a). As a
consequence, in this case the calculated surface energy corresponds to an
average value of the zinc- and oxygen-terminated surfaces. In all
calculations, slabs of 6-12 atomic layers are separated by about 15 Å of
vacuum, and the atomic coordinates of all ions are fully relaxed. The
calculated surface energies, Tab. 2, are fully consistent with the existing
computational evidence Claeyssens  _et al._ (2005); Freeman  _et al._ (2006);
Goniakowski  _et al._ (2007b); Morgan (2009); Demiroglu and Bromley (2013);
Sponza  _et al._ (2015).

Table 2: Calculated (PW91) surface energies of low index surfaces of the six ZnO polymorphs $E_{\mathrm{surf}}$ (J/m2). |  orientation |  $E_{\mathrm{surf}}$
---|---|---
WRZ |  (10-10) |  0.89
WRZ |  (11-20) |  0.93
WRZ |  (0001) |  1.30
ZBL |  (110) |  0.92
BCT |  (010) |  0.49
BCT |  (001) |  1.14
SOD |  (001) |  1.00
HBN |  (0001) |  0.36
CUB |  (001) |  0.48

Three PLIPs are developed with different training datasets. Indeed, we train a
first PLIP on a database composed only of the 6 bulk polymorphs melted up to
5000 K, a second based on ZnO surfaces heated up to 2000 K along with the
first PLIP database, and a third one containing additional amorphous bulk
obtained after rapidly cooling liquid structures obtained at 2000 K. Those
three PLIPs are respectively denoted V1, V2, and V3 in the following. To
validate our PLIPs, we perform a series of calculations where bulk and surface
properties are computed for the six ZnO polymorphs, and compared to DFT
results. We employ the conjugate gradient algorithm, with a force convergence
criterion of 10-9 eV/Å. Moreover, we investigate the efficiency of our PLIPs
in the context of disordered ZnO structures, such as liquid and amorphous
phases. For this purpose, simulations are carried out by starting with 500
atoms randomly disposed and heated at 3000 K, and then cooled to 300 K with a
rate of 1800 K/ns. Final structures are analyzed by computing radial
distribution functions (RDF).

###  II.3 Molecular dynamics simulations of freezing

After developing our PLIP models, we perform freezing molecular dynamics
simulations. In particular, the system is initialized with 1000 randomly
positioned atoms of zinc and the same number of oxygen atoms. Then, atomic
positions are first optimized to remove large forces due to overlapping atoms.
We further proceed to a first equilibration at 5000 K during 1 ps in the NVT
ensemble with a Nose-Hoover thermostat. From preliminary freezing simulations
(not shown) we found out that crystallization only occurs in the temperature
range of 1500 K to 1000 K. Therefore, for our reported simulations, the
freezing is operated in two steps. In the first step, the system is quenched
from 2500 K to 1500 K during 100 ps which allows us to quickly obtain a first
liquid structure. In the second step, the system is further slowly cooled down
from 1500 K to 1000 K during 10 ns. These two cooling steps are both performed
in the NPT ensemble at P = 1 bar. Finally, the atomic positions and box size
of the obtained structures are optimized using a conjugate gradient algorithm.

The Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) package
Plimpton (1995) with a custom-made pair style for PLIP modeling is used in all
the simulations.

###  II.4 Data-driven approach for the structural analysis

Analysis of the obtained structures and the crystallization dynamics requires
the usage of novel order parameters which are able to distinguish between each
of the six crystal polymorphs of zinc oxide. Such a challenging task is
clearly out of the reach of classical order parameter including Common
neighbor analysis, Voronoi cell topology, or Polyhedral template matching
which are mostly designed for monodisperse systems and for classical crystal
structures like face-centered cubic and body centered cubic. In the case of
zinc oxide crystallization, the very rich structural landscape requires the
use of a structurally agnostic order parameter, namely one that is
transferable for all kinds of crystal structures.

In this context, we develop a data-driven approach for the structural
analysis. In practice, averaged Steinhardt parameters denoted
$\overline{q_{l}}$ with $l\in[3:6]$ are computed for zinc atoms when
considering the entire system (zinc and oxygen atoms) and only inter-species
distances (Zn-Zn and O-O). We use the pyscal Python package Menon  _et al._
(2019) with the adaptive methods to compute the neighbor lists. This
calculation leads to a Steinhardt vector denoted $\bm{Q}$ of $12$ components
(6 for the entire system and 6 for inter-species distances) for each atom.
Then, we compute the same vector in the six bulk polymorphs to construct a
reference database averaged over all the zinc atoms in those crystal
structures. If a studied atom has the same local structure as one in the bulk
polymorph, its Steinhardt vector must be close to that of the database which
can be measured using Euclidean distance in vector space that we denote
$D_{Y-X}$. In Fig. 1, we test those Euclidean distances on nanoparticles of
known structures and show that the method is able to identify the correct
structure. In addition, the method retrieves that the center of the
nanoparticle is more similar to the bulk than its edges.

Finally, in order to quantify these Euclidean distances, we compute a
likelihood parameter denoted $S$ measured along a path between two reference
structures named A and B using:

$S=\frac{e^{-\lambda D_{AX}}+2e^{-\lambda D_{BX}}}{e^{-\lambda
D_{AX}}+e^{-\lambda D_{BX}}}$ (4)

where $\lambda$ is adjusted so that $D_{AB}=2.3/\lambda$. As such, we obtain a
score that quantifies the local structure: an atom which is similar to
structure A (resp. B) gives $S\sim 1$ (resp. $S\sim 2$). Therefore, we
consider that the value of $1.5$ discriminates between the two types of
structure.

In retrospective, our approach combines the use of (1) a vector of Steinhardt
parameters which was already suggested for colloidal systems and amorphous
models (Boattini  _et al._ , 2020, 2019), (2) a similarity measurement with
the identity function as Kernel, and (3) a path collective variable which was
previously employed for the construction of order parameters based on
permutation invariant vectors Pipolo  _et al._ (2017).

Figure 1: Structural analysis of three different nanoparticles obtained from
Viñes et al. Viñes  _et al._ (2017) that are respectively in the WRZ (a), BCT
(b) and SOD (c) crystal phases. We measured the distances with respect to the
same three different crystal phases in bulk.

##  III Results

###  III.1 Validation of the PLIP model

####  III.1.1 Bulk crystals

As a first step to validate the developed PLIPs, the crystal structures of the
six polymorphs were optimized using the three PLIPs, V1, V2, and V3. The
resulting lattice parameter $a$ and cohesion energy $E_{coh}$ relative to the
lowest-energy WRZ phase were compared to those obtained using DFT in order to
evaluate the reproduction accuracy of the PLIPs. Errors with respect to DFT
values for these two quantities are displayed in Fig. 2, along with results
obtained with the Buckingham potential.

Regarding the three PLIPs, they give approximately the same results, with the
V2 and V3 PLIPs being slightly more accurate than the V1 one. The already
accurate results obtained with V1 indicate that training with only bulk ZnO
structures is sufficient to reach good model for the bulk crystals.

Accuracy on ZnO surfaces was also investigated and the errors on the surface
energy are computed both using PLIPs and the Buckingham potential [See Fig.
2.(c)]. The PLIP errors are generally much lower than what is obtained with
the Buckingham potential with an improvement of a factor between 3 and 8.5.
When comparing each PLIP version, we note that while V2 and V3 are accurate by
construction, V1 in most cases leads to very good results in terms of surface
energy although it was not trained with any surface structures.

Figure 2: Errors (%) with respect to the DFT values of bulk lattice parameter
$a$, cohesion energy $E_{coh}$ relative to the wurtzite polymorph, and the
relative surface energy $\gamma$ determined from the relaxed structures
obtained with the three PLIPs as well as with the Buckingham potential.

####  III.1.2 Liquid

After the crystal properties of ZnO, we were interested in measuring the
accuracy of the model for amorphous structures. For this purpose, a hot liquid
was cooled until reaching the ambient temperature with a cooling rate allowing
for the formation of the amorphous form of ZnO. The final structure is
analyzed by calculating the partial RDFs, as displayed in Fig. 3, when using
the PLIPs, the Buckingham potential, as well as DFT calculations from
Ref.Mora-Fonz and Shluger (2019) for comparison. The RDFs coming from the PLIP
structures fit perfectly those of DFT structures (solid black line), while a
mismatch is clearly visible between the DFT and the Buckingham (dashed black
line) curves. Moreover, even if all the PLIP RDFs are similar, a slight
advantage is given to the V3 curve (green line), especially at the first peak.
However, it appears that the training of V1 and V2 models on hot liquid
structures seems sufficient to describe the amorphous structures stabilized at
room temperature.

Figure 3: Partial radial distribution functions $g(r)$ of final structures
obtained after cooling a hot liquid. Structures calculated with the three
PLIPs are compared with those computed using DFT (solid black lines) and the
Buckingham potential (dashed black lines).

####  III.1.3 Nanostructures

Figure 4: (a) Relative error on the clusters’s energy per unit of ZnO relative
to the Wurtzite bulk energy. Corresponding DFT energies are extracted from
Viñes et al. Viñes  _et al._ (2017). (b) Snapshots representing some of the
studied clusters.

For this last test of our PLIP models, we employed clusters of ZnO obtained by
Viñes et al. Viñes  _et al._ (2017). They focused on three polymorphs, namely
BCT, WRZ and SOD, and also worked with single-caged (SC) structures. They
managed to find the most probable low energy structures for systems made of
tens to hundreds of atoms and we computed the single-point energy with our
PLIP models. Fig. 4 shows that the three versions of PLIP give good results
when comparing to DFT calculations Viñes  _et al._ (2017). On the contrary,
the Buckingham model leads to relative errors that can be larger than 50%.
Surprisingly, while the largest clusters exhibit bulk-like structures, the
smallest ones are qualitatively different. Yet, our PLIP models are able to
accurately reproduce the energy even for those smallest clusters which is an
evidence for the transferability towards untrained structures. In particular,
the single cage SC-$108$ and SC-$90$ do not show any bulk-like structure since
these are shell cluster, and the PLIP models still give very small errors.

Before closing on this validation section, it is worth noting that although
our PLIP models were only trained with DFT forces, they are still able to
reach very accurate results for energetic properties (surface, bulk, and
nanoparticles). From now on, because V3 is built with a more complete
database, calculations presented in the remaining article will only employ the
V3 potential.

###  III.2 Freezing simulations and crystallization mechanisms

Following this first section dedicated to the validation of the PLIP model, we
will now describe results obtained with freezing simulations. As a starting
point, Fig. 5.(a,b) shows the structures obtained at the end of the freezing
simulation. When comparing results obtained with PLIP and the Buckhingham
model, it is difficult to differentiate between the two images. In order to
provide a more quantitative picture, the Euclidean distances in Steinhardt
vector space between the obtained structures and those in the bulk database
are computed and shown in Fig 5.c. It appears that the smallest distance for
the PLIP results is obtained with the WRZ structure while for the Buckingham
potential, it corresponds to the BCT structure. This quantitative result
suggests that the classical potential leads to a qualitatively different
result, $i.e.$ the growth of a BCT structure against the Wurtzite structure
observed in both experiments and DFT calculations. Being able to crystallize
the most stable structure is a strong justification for using the PLIP model
over a more classical potential.

Figure 5: Structures obtained after freezing simulations using (a) the PLIP
model and (b) the Buckingham potential, after a slow cooling during $10$ ns.
(c) Euclidean distances between the obtained structures and the bulk perfect
crystals computed in the Steinhardt vector space.

Thus, one of the key advantages for using the PLIP model over classical
interactions such as the Buckingham potential is that, along with the
quantitatively better results for bulk and surface properties, it also leads
to qualitatively better outcome of freezing simulations. In addition, when
comparing to DFT, the main advantage is the possibility for larger scale
simulations. This is essential in studies of crystallization where freezing
simulations have to cover more than a few nanoseconds because shorter ones
result in amorphous structures. Furthermore, much larger unit cells allowed by
our PLIP help avoiding any nonphysical behavior due to interactions between
nucleation cores in neighbouring cells.

As to get an insight into the nucleation scenarios, three individual runs were
performed using the same thermodynamic conditions and varying only the initial
atom positions and velocities. We investigate in more detail the behavior of
the number of atoms in the predominant WRZ and BCT environments which are
measured using the path order parameter $S$. We observed two types of scenario
regrouped in Fig. 6 and Fig. 7.

Figure 6: (a,c) Number of atoms detected in WRZ or BCT phases along the
cooling path $i.e.$ with respect to time and temperature, for two MD runs. The
beginnings of the curve increase are highlighted in grey and zoomed in (b,d).
Colored dashed lines point out where the structures (e)-(i) are picked along
the freezing simulations. Atoms identified as WRZ and BCT according to the
likelihood parameter $S$ are colored in green and blue, respectively.

On the one hand, three stages can be identified in Fig. 6.(a,c): (1) an abrupt
and quick increase around 7 ns and 1150 K marks the nucleation and growth of
WRZ and BCT structures from liquid ZnO, (2) a plateau more or less broad
indicates the stabilization of coexisting BCT and WRZ phases [See Fig. 6.i],
and (3) a phase transformation of BCT into WRZ visible by the sudden drop of
the curve of the former, compensated by the increase of the latter. In the
following, we will focus particularly on what happens during the first stage
(1) [See Fig. 6.(b,d)]. First, as it can be expected, growth of the WRZ phase
is initiated by a nucleation of the first WRZ grain. The behavior of $N_{WRZ}$
at the very beginning of the MD simulation, Fig. 6.(b,d), shows that several
small WRZ ZnO clusters form and dissolve in the ZnO liquid before the actual
growth starts. From a structural point of view, its stability seems to be
enhanced by the formation of three-membered ZnO rings [See Fig. 6.(e,f)],
characteristic of the Wurtzite (0001) surface. Indeed, in absence of such
rings, similarly-sized ZnO clusters tend to dissolve. Furthermore, the
emerging WRZ clusters are surrounded by a BCT phase taking the form of
relatively thin and less crystallized structures [See Fig. 6.(g,h)]. Then,
during the growth, the grain becomes significantly elongated, with virtually
parallel (0001) facets at its opposite sides. This suggests that
crystallization along the (0001) atomic planes is favored at this growth
stage. Altogether, our observations suggest the following first scenario of
nucleation: (1) Stochastic appearance of unstable WRZ clusters that do not
exhibit three-membered rings, (2) Emergence of a stable WRZ cluster made of
three-membered rings and surrounded by a shell of BCT phase, (3) Continuous
growth of WRZ with a reminiscent BCT shell, (4) After the main part of the
growth processes, unstable BCT structures occupying the edges of the periodic
WRZ are consumed by the WRZ crystal. From a fundamental perspective, the
appearance of such pre-nucleation cluster made of a core of the stable phase
surrounded by a shell of an unstable phase contradicts the classical
nucleation scenario in which the first emerging crystal is made of the most
stable phase only. Such non-classical nucleation scenario has been also
observed in simpler systems described by hard spheres and Lennard-Jones
interactions. Because we managed to observe this in a system as complex as
zinc oxide modeled with machine-learning interaction potential, we confirm the
universality of this phenomenon.

Figure 7: (a) Number of atoms detected in WRZ or BCT phases along the cooling
path for the third MD run. Colored dashed lines point out where the structures
(b)-(d) are picked along the freezing simulation. Atoms identified as WRZ and
BCT according to the likelihood parameter $S$ are colored in green and blue,
respectively.

On the other hand, Fig. 7.(a) shows a completely different scenario of
nucleation. Indeed, while $N_{WRZ}$ and $N_{BCT}$ were previously increasing
at the same time and temperature, here, the unstable BCT emerges sooner than
its stable WRZ counterpart. Fig. 7.(b,c) confirm that the emerging crystal is
indeed made of mostly BCT phase. Then, once the BCT is large enough, within
its structure, one can observe a solid-solid transition where BCT turns into
WRZ. In this case, the BCT structure becomes porous, with multiple embedded
WRZ clusters and is not stable enough to continue persisting. It is remarkable
that while the initial BCT structure is much larger than in the first
scenario, it can not survive as long because of this porous growth of WRZ. As
in the first scenario, we observe an other example of non-classical nucleation
with a two-step mechanisms: (1) Metastable BCT crystals constitute the
nucleation core and (2) The stable WRZ crystal only appears within large BCT
crystals. An interesting feature of this second scenario is that we no longer
observe growth of elongated WRZ grains since it originates from multiple
embedded smaller clusters. The elongated grains that were previously observed
should therefore be favored by the presence of BCT shells.

##  IV Conclusions

Despite the key technological importance of ZnO in many applications, its
crystallization features are at present relatively poorly understood. This is
due to overwhelming computational effort necessary for dedicated simulations,
which excludes a direct use of DFT-based methods and requires fine interatomic
potentials, able to correctly account for the structural diversity of ZnO.

With the goal to access the structure of critical nucleus in a homogeneous
liquid ZnO phase and to follow the subsequent growth of the ZnO crystal, we
have constructed and validated a new robust Machine Learning Interatomic
Potential (MLIP) suitable for large-scale simulations (nanosecond-long
molecular dynamics of several hundreds of atoms). The potential is based on
physically-motivated mathematical formulation and a constrained LassoLars
method is used for identification and adjustment of its pertinent parameters.
The training database was composed of DFT-GGA results on a variety of ordered
crystalline structures, but includes also disordered and amorphous
configurations, such as obtained in high-temperature molecular dynamics
simulations. We show that the new interaction potential successfully
reproduces the fine structural and energetic DFT-GGA characteristics of six of
the most stable ZnO bulk polymorphs as well as of their low-index surfaces.
Moreover, it also improves the simulated structural characteristics of liquid
ZnO such that they become quasi-indistinguishable from the DFT ones. Finally,
the obtained potential was successfully tested on completely untrained
structures made of nanoparticles in different crystal phases. Altogether,
these extensive testings demonstrate the reliability and the transferability
of our PLIP model which was then employed to study the crystallization of zinc
oxide.

In particular, we have performed molecular dynamics simulations on freezing of
bulk ZnO liquid. First, our results show that while traditional interaction
potentials predict the formation and growth of meta-stable crystal nuclei, the
PLIP model leads to a qualitatively different behavior in which the most
stable crystal structure emerges. Analysis of the crystalline phase in such
complex structural landscape was made posible by employing a data-driven
approach based on Steinhardt’s parameters. Then, we were able to access at the
atomistic level the birth of the initial zinc oxide crystal in the surrounding
liquid and to follow its behavior upon further growth. While performing three
independent runs at the same thermodynamic conditions, we observed two
different scenarios which are both in contradiction to classical nucleation
theory: (1) Prenucleation clusters made of two coexisting phases, WRZ in the
center and BCT in the shell and (2) Two-step nucleation process with
metastable BCT emerging first and turning into stable WRZ afterwards. Similar
scenarios of nucleation were previously observed in much simpler systems
including Lennard-Jones, hard-sphere, and pure metals. In our case, the
combination of our Physical LassoLars interaction potential that allows for a
quantum-accurate modeling of zinc oxide with our data-driven approach for a
refined structural analysis of crystal phases leads us to confirm the
universality of seminal predictions of non-classical nucleation processes.

## Acknowledgement

JL acknowledges financial support of the Fonds de la Recherche Scientifique -
FNRS. Computational resources have been provided by the Consortium des
Equipements de Calcul Intensif (CECI), by the Fédération Lyonnaise de
Modélisation et Sciences Numériques (FLMSN) and by he Regional Computer Center
CALMIP in Toulouse. JL thanks David Mora Fonz and Anup Pandey for fruitful
discussions.

## References

  * Lutsko and Lam (2020) J. F. Lutsko and J. Lam, Phys. Rev. E 101, 052122 (2020).
  * Lutsko (2019) J. F. Lutsko, Sci. Adv. 5, 7399 (2019).
  * Wolde and Frenkel (1997) P. R. t. Wolde and D. Frenkel, Science 277, 1975 (1997).
  * Auer and Frenkel (2001) S. Auer and D. Frenkel, Nature 409, 1020 (2001).
  * Ten Wolde  _et al._ (1995) P. R. Ten Wolde, M. J. Ruiz-Montero, and D. Frenkel, Phys. Rev. Lett. 75, 2714 (1995).
  * Sleutel and Van Driessche (2014) M. Sleutel and A. E. S. Van Driessche, Proc. Natl. Acad. Sci. U.S.A. 111, E546 (2014).
  * Van Driessche  _et al._ (2018) A. E. S. Van Driessche, N. Van Gerven, P. H. H. Bomans, R. R. M. Joosten, H. Friedrich, D. Gil-Carton, N. A. J. M. Sommerdijk, and M. Sleutel, Nature 556, 89 (2018).
  * Menon  _et al._ (2020) S. Menon, G. Díaz Leines, R. Drautz, and J. Rogal, J. Chem. Phys. 153, 104508 (2020).
  * Lam and Lutsko (2018) J. Lam and J. F. Lutsko, Nanoscale 10, 4921 (2018).
  * Amodeo  _et al._ (2020) J. Amodeo, F. Pietrucci, and J. Lam, J. Phys. Chem. Lett. 11, 8060 (2020).
  * Desgranges and Delhommelle (2007a) C. Desgranges and J. Delhommelle, Phys. Rev. Lett. 98, 235502 (2007a).
  * Desgranges and Delhommelle (2006) C. Desgranges and J. Delhommelle, J. Am. Chem. Soc. 128, 15104 (2006).
  * Yamazaki  _et al._ (2017) T. Yamazaki, Y. Kimura, P. G. Vekilov, E. Furukawa, M. Shirai, H. Matsumoto, A. E. S. Van Driessche, and K. Tsukamoto, Proc. Natl. Acad. Sci. U.S.A. 114, 2154 (2017).
  * Sleutel  _et al._ (2014) M. Sleutel, J. Lutsko, A. E. S. Van Driessche, M. A. Durán-Olivencia, and D. Maes, Nat. Commun. 5, 1 (2014).
  * Ogata  _et al._ (2020) A. F. Ogata, A. M. Rakowski, B. P. Carpenter, D. A. Fishman, J. G. Merham, P. J. Hurst, and J. P. Patterson, J. Am. Chem. Soc. 142, 1433 (2020).
  * van Blaaderen  _et al._ (1997) A. van Blaaderen, R. Ruel, and P. Wiltzius, Nature 385, 321 (1997).
  * Leunissen  _et al._ (2005) M. E. Leunissen, C. G. Christova, A.-P. Hynninen, C. P. Royall, A. I. Campbell, A. Imhof, M. Dijkstra, R. van Roij, and A. van Blaaderen, Nature 437, 235 (2005).
  * Velikov  _et al._ (2002) K. P. Velikov, C. G. Christova, R. P. A. Dullens, and A. van Blaaderen, Science 296, 106 (2002).
  * de Nijs  _et al._ (2015) B. de Nijs, S. Dussi, F. Smallenburg, J. D. Meeldijk, D. J. Groenendijk, L. Filion, A. Imhof, A. van Blaaderen, and M. Dijkstra, Nat. Mater. 14, 56 (2015).
  * Sosso  _et al._ (2016) G. C. Sosso, J. Chen, S. J. Cox, M. Fitzner, P. Pedevilla, A. Zen, and A. Michaelides, Chem. Rev. 116, 7078 (2016).
  * Coli and Dijkstra (2021) G. M. Coli and M. Dijkstra, ACS Nano 15, 4335 (2021).
  * Russo and Tanaka (2012) J. Russo and H. Tanaka, Sci. Rep. 2, 1 (2012).
  * Leoni and Russo (2021) F. Leoni and J. Russo, Phys. Rev. X 11, 031006 (2021).
  * Lin  _et al._ (2018) C. Lin, G. Corem, O. Godsi, G. Alexandrowicz, G. R. Darling, and A. Hodgson, J. Am. Chem. Soc. 140, 15804 (2018).
  * Metya and Molinero (2021) A. K. Metya and V. Molinero, J. Am. Chem. Soc. 143, 4607 (2021).
  * Qiu and Molinero (2018) Y. Qiu and V. Molinero, J. Phys. Chem. Lett. 9, 5179 (2018).
  * Qiu  _et al._ (2019) Y. Qiu, A. Hudait, and V. Molinero, J. Am. Chem. Soc. 141, 7439 (2019).
  * Lupi  _et al._ (2014) L. Lupi, A. Hudait, and V. Molinero, J. Am. Chem. Soc. 136, 3156 (2014).
  * Moore and Molinero (2011) E. B. Moore and V. Molinero, Nature 479, 506 (2011).
  * Lupi  _et al._ (2017) L. Lupi, A. Hudait, B. Peters, M. Grünwald, R. Gotchy Mullen, A. H. Nguyen, and V. Molinero, Nature 551, 218 (2017).
  * Russo  _et al._ (2014) J. Russo, F. Romano, and H. Tanaka, Nat. Mater. 13, 733 (2014).
  * Desgranges and Delhommelle (2007b) C. Desgranges and J. Delhommelle, J. Am. Chem. Soc. 129, 7012 (2007b).
  * Combettes  _et al._ (2020) S. Combettes, J. Lam, P. Benzo, A. Ponchet, M.-J. Casanove, F. Calvo, and M. Benoit, Nanoscale 12, 18079 (2020).
  * Laurens  _et al._ (2020) G. Laurens, D. Amans, J. Lam, and A.-R. Allouche, Phys. Rev. B 101, 045427 (2020).
  * Calvo  _et al._ (2017) F. Calvo, N. Combe, J. Morillo, and M. Benoit, J. Phys. Chem. C 121, 4680 (2017).
  * Zipoli and Curioni (2013) F. Zipoli and A. Curioni, New J. Phys. 15, 123006 (2013).
  * Rajasekaran  _et al._ (2016) G. Rajasekaran, R. Kumar, and A. Parashar, Mater. Res. Express 3, 035011 (2016).
  * Senftle  _et al._ (2016) T. P. Senftle, S. Hong, M. M. Islam, S. B. Kylasa, Y. Zheng, Y. K. Shin, C. Junkermeier, R. Engel-Herbert, M. J. Janik, H. M. Aktulga, T. Verstraelen, A. Grama, and A. C. T. van Duin, npj Comput. Mater. 2, 1 (2016).
  * Brommer  _et al._ (2015) P. Brommer, A. Kiselev, D. Schopf, P. Beck, J. Roth, and H.-R. Trebin, Model. Simul. Mater. Sci. Eng. 23, 074002 (2015).
  * Vashishta  _et al._ (2008) P. Vashishta, R. K. Kalia, A. Nakano, and J. P. Rino, J. Appl. Phys. 103, 083504 (2008).
  * Behler and Parrinello (2007) J. Behler and M. Parrinello, Phys. Rev. Lett. 98, 146401 (2007).
  * Bartók  _et al._ (2010) A. P. Bartók, M. C. Payne, R. Kondor, and G. Csányi, Phys. Rev. Lett. 104, 136403 (2010).
  * Seko  _et al._ (2015) A. Seko, A. Takahashi, and I. Tanaka, Phys. Rev. B 92, 054113 (2015).
  * Seko  _et al._ (2014) A. Seko, A. Takahashi, and I. Tanaka, Phys. Rev. B 90, 024101 (2014).
  * Takahashi  _et al._ (2017) A. Takahashi, A. Seko, and I. Tanaka, Phys. Rev. Materials 1, 063801 (2017).
  * Seko  _et al._ (2019) A. Seko, A. Togo, and I. Tanaka, Phys. Rev. B 99, 214108 (2019).
  * Goryaeva  _et al._ (2019) A. M. Goryaeva, J.-B. Maillet, and M.-C. Marinica, Comput. Mater. Sci. 166, 200 (2019).
  * Benoit  _et al._ (2020) M. Benoit, J. Amodeo, S. Combettes, I. Khaled, A. Roux, and J. Lam, Mach. Learn.: Sci. Technol. 2, 025003 (2020).
  * Thompson  _et al._ (2015) A. P. Thompson, L. P. Swiler, C. R. Trott, S. M. Foiles, and G. J. Tucker, J. Comput. Phys. 285, 316 (2015).
  * Wood and Thompson (2018) M. A. Wood and A. P. Thompson, J. Chem. Phys. , 241721 (2018).
  * Chmiela  _et al._ (2017) S. Chmiela, A. Tkatchenko, H. E. Sauceda, I. Poltavsky, K. T. Schütt, and K.-R. Müller, Sci. Adv. 3, e1603015 (2017).
  * Chmiela  _et al._ (2018) S. Chmiela, H. E. Sauceda, K.-R. Müller, and A. Tkatchenko, Nat. Commun. 9, 1 (2018).
  * Shapeev (2016) A. V. Shapeev, Multiscale Model. Simul. (2016).
  * Novikov  _et al._ (2020) I. S. Novikov, K. Gubaev, E. V. Podryabinkin, and A. V. Shapeev, Mach. Learn.: Sci. Technol. 2, 025002 (2020).
  * Novoselov  _et al._ (2019) I. I. Novoselov, A. V. Yanilkin, A. V. Shapeev, and E. V. Podryabinkin, Comput. Mater. Sci. 164, 46 (2019).
  * Takahashi  _et al._ (2018) A. Takahashi, A. Seko, and I. Tanaka, J. Chem. Phys. 148, 234106 (2018).
  * Zeni  _et al._ (2018) C. Zeni, K. Rossi, A. Glielmo, Á. Fekete, N. Gaston, F. Baletto, and A. De Vita, J. Chem. Phys. 148, 241739 (2018).
  * Botu  _et al._ (2017) V. Botu, R. Batra, J. Chapman, and R. Ramprasad, J. Phys. Chem. C 121, 511 (2017).
  * Bereau  _et al._ (2018) T. Bereau, R. A. DiStasio, A. Tkatchenko, and O. A. von Lilienfeld, J. Chem. Phys. , 241706 (2018).
  * Sauceda  _et al._ (2019) H. E. Sauceda, S. Chmiela, I. Poltavsky, K.-R. Müller, and A. Tkatchenko, J. Chem. Phys. 150, 114102 (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).
  * Veit  _et al._ (2019) M. Veit, S. K. Jain, S. Bonakala, I. Rudra, D. Hohl, and G. Csányi, J. Chem. Theory Comput. 15, 2574 (2019).
  * Nguyen  _et al._ (2018) T. T. Nguyen, E. Székely, G. Imbalzano, J. Behler, G. Csányi, M. Ceriotti, A. W. Götz, and F. Paesani, J. Chem. Phys. , 241725 (2018).
  * Bartók  _et al._ (2013) A. P. Bartók, M. J. Gillan, F. R. Manby, and G. Csányi, Phys. Rev. B 88, 054104 (2013).
  * Morawietz  _et al._ (2012) T. Morawietz, V. Sharma, and J. Behler, J. Chem. Phys. 136, 064103 (2012).
  * Natarajan and Behler (2016) S. K. Natarajan and J. Behler, Phys. Chem. Chem. Phys. 18, 28704 (2016).
  * Morawietz and Behler (2013) T. Morawietz and J. Behler, J. Phys. Chem. A 117, 7356 (2013).
  * Deringer and Csányi (2017) V. L. Deringer and G. Csányi, Phys. Rev. B 95, 094203 (2017).
  * Bartók  _et al._ (2018) A. P. Bartók, J. Kermode, N. Bernstein, and G. Csányi, Phys. Rev. X 8, 041048 (2018).
  * Caro  _et al._ (2018) M. A. Caro, V. L. Deringer, J. Koskinen, T. Laurila, and G. Csányi, Phys. Rev. Lett. 120, 166101 (2018).
  * Deringer  _et al._ (2018a) V. L. Deringer, N. Bernstein, A. P. Bartók, M. J. Cliffe, R. N. Kerber, L. E. Marbella, C. P. Grey, S. R. Elliott, and G. Csányi, J. Phys. Chem. Lett. 9, 2879 (2018a).
  * Deringer  _et al._ (2018b) V. L. Deringer, M. A. Caro, R. Jana, A. Aarva, S. R. Elliott, T. Laurila, G. Csányi, and L. Pastewka, Chem. Mater. 30, 7438 (2018b).
  * Sosso  _et al._ (2018) G. C. Sosso, V. L. Deringer, S. R. Elliott, and G. Csányi, Mol. Simul. 44, 866 (2018).
  * Jinnouchi  _et al._ (2019) R. Jinnouchi, J. Lahnsteiner, F. Karsai, G. Kresse, and M. Bokdam, Phys. Rev. Lett. 122, 225701 (2019).
  * Sivaraman  _et al._ (2020) G. Sivaraman, A. N. Krishnamoorthy, M. Baur, C. Holm, M. Stan, G. Csányi, C. Benmore, and Á. Vázquez-Mayagoitia, npj Comput. Mater. 6, 1 (2020).
  * Artrith  _et al._ (2013) N. Artrith, B. Hiller, and J. Behler, Phys. Status Solidi B 250, 1191 (2013).
  * Eckhoff  _et al._ (2020a) M. Eckhoff, K. N. Lausch, P. E. Blöchl, and J. Behler, J. Chem. Phys. 153, 164107 (2020a).
  * Artrith  _et al._ (2011) N. Artrith, T. Morawietz, and J. Behler, Phys. Rev. B 83, 153101 (2011).
  * Sundararaman  _et al._ (2018) S. Sundararaman, L. Huang, S. Ispas, and W. Kob, J. Chem. Phys. 148, 194504 (2018).
  * Eckhoff  _et al._ (2020b) M. Eckhoff, F. Schönewald, M. Risch, C. A. Volkert, P. E. Blöchl, and J. Behler, Phys. Rev. B 102, 174102 (2020b).
  * Artrith  _et al._ (2017) N. Artrith, A. Urban, and G. Ceder, Phys. Rev. B 96, 014112 (2017).
  * Artrith and Urban (2016) N. Artrith and A. Urban, Comput. Mater. Sci. 114, 135 (2016).
  * Ko  _et al._ (2021a) T. W. Ko, J. A. Finkler, S. Goedecker, and J. Behler, Nat. Commun. 12, 1 (2021a).
  * Ko  _et al._ (2021b) T. W. Ko, J. A. Finkler, S. Goedecker, and J. Behler, Acc. Chem. Res. 54, 808 (2021b).
  * Navrotsky  _et al._ (2008) A. Navrotsky, L. Mazeina, and J. Majzlan, Science 319, 1635 (2008).
  * Machala  _et al._ (2011) L. Machala, J. Tuček, and R. Zbořil, Chem. Mater. 23, 3255 (2011).
  * Sponza  _et al._ (2015) L. Sponza, J. Goniakowski, and C. Noguera, Phys. Rev. B 91, 075126 (2015).
  * Tsybulya and Kryukova (2008) S. V. Tsybulya and G. N. Kryukova, Phys. Rev. B 77, 024112 (2008).
  * Mirzaeifard  _et al._ (2020) Z. Mirzaeifard, Z. Shariatinia, M. Jourshabani, and S. M. Rezaei Darvishi, Ind. Eng. Chem. Res. 59, 15894 (2020).
  * Mahmood  _et al._ (2021) N. Mahmood, H. Khan, K. Tran, P. Kuppe, A. Zavabeti, P. Atkin, M. B. Ghasemian, J. Yang, C. Xu, S. A. Tawfik, M. J. S. Spencer, J. Z. Ou, K. Khoshmanesh, C. F. McConville, Y. Li, and K. Kalantar-Zadeh, Mater. Today 44, 69 (2021).
  * Xiong (2013) H.-M. Xiong, Adv. Mater. 25, 5329 (2013).
  * Abebe  _et al._ (2020) B. Abebe, E. A. Zereffa, A. Tadesse, and H. C. A. Murthy, Nanoscale Res. Lett. 15, 1 (2020).
  * Yuan  _et al._ (2019) H. Yuan, S. A. A. A. Aljneibi, J. Yuan, Y. Wang, H. Liu, J. Fang, C. Tang, X. Yan, H. Cai, Y. Gu, S. J. Pennycook, J. Tao, and D. Zhao, Adv. Mater. 31, 1807161 (2019).
  * Zeni  _et al._ (2021) C. Zeni, K. Rossi, A. Glielmo, and S. de Gironcoli, J. Chem. Phys. 154, 224112 (2021).
  * Kresse and Hafner (1993) G. Kresse and J. Hafner, Phys. Rev. B 47, 558 (1993).
  * Kresse and Furthmuller (1996) G. Kresse and J. Furthmuller, Phys. Rev. B 54, 11169 (1996).
  * Perdew  _et al._ (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Physical review letters 77, 3865 (1996).
  * Kresse and Joubert (1999) G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
  * Binks and Grimes (1993) D. J. Binks and R. W. Grimes, J. Am. Ceram. Soc. 76, 2370 (1993).
  * Wang  _et al._ (2010) J. Wang, P. Xiao, M. Zhou, Z. R. Wang, and F. J. Ke, J. Appl. Phys. 107, 023512 (2010).
  * Lin  _et al._ (2013) K.-H. Lin, S.-J. Sun, S.-P. Ju, J.-Y. Tsai, H.-T. Chen, and J.-Y. Hsieh, J. Appl. Phys. 113, 073512 (2013).
  * Kulkarni  _et al._ (2005) A. J. Kulkarni, M. Zhou, and F. J. Ke, Nanotechnology 16, 2749 (2005).
  * Kulkarni  _et al._ (2006) A. J. Kulkarni, M. Zhou, K. Sarasamak, and S. Limpijumnong, Phys. Rev. Lett. 97, 105502 (2006).
  * Goniakowski  _et al._ (2007a) J. Goniakowski, F. Finocchi, and C. Noguera, Reports on Progress in Physics 71, 016501 (2007a).
  * Claeyssens  _et al._ (2005) F. Claeyssens, C. L. Freeman, N. L. Allan, Y. Sun, M. N. R. Ashfold, and J. H. Harding, J. Mater. Chem. 15, 139 (2005).
  * Freeman  _et al._ (2006) C. L. Freeman, F. Claeyssens, N. L. Allan, and J. H. Harding, Phys. Rev. Lett. 96, 066102 (2006).
  * Goniakowski  _et al._ (2007b) J. Goniakowski, C. Noguera, and L. Giordano, Phys. Rev. Lett. 98, 205701 (2007b).
  * Morgan (2009) B. J. Morgan, Phys. Rev. B 80, 174105 (2009).
  * Demiroglu and Bromley (2013) I. Demiroglu and S. T. Bromley, Phys. Rev. Lett. 110, 245501 (2013).
  * Plimpton (1995) S. Plimpton, J. Comput. Phys. 117, 1 (1995).
  * Menon  _et al._ (2019) S. Menon, G. D. Leines, and J. Rogal, Journal of Open Source Software 4, 1824 (2019).
  * Boattini  _et al._ (2020) E. Boattini, S. Marín-Aguilar, S. Mitra, G. Foffi, F. Smallenburg, and L. Filion, Nat. Commun. 11, 1 (2020).
  * Boattini  _et al._ (2019) E. Boattini, M. Dijkstra, and L. Filion, J. Chem. Phys. 151, 154901 (2019).
  * Pipolo  _et al._ (2017) S. Pipolo, M. Salanne, G. Ferlat, S. Klotz, A. M. Saitta, and F. Pietrucci, Phys. Rev. Lett. 119, 245701 (2017).
  * Viñes  _et al._ (2017) F. Viñes, O. Lamiel-Garcia, F. Illas, and S. T. Bromley, Nanoscale 9, 10067 (2017).
  * Mora-Fonz and Shluger (2019) D. Mora-Fonz and A. L. Shluger, Phys. Rev. B 99, 014202 (2019).