---
abstract: 'Atomistic molecular dynamics simulation is an important tool for predicting materials properties. Accuracy depends crucially on the model for the interatomic potential. The gold standard would be quantum mechanics (QM) based force calculations, but such a first-principles approach becomes prohibitively expensive at large system sizes. Efficient machine learning models (ML) have become increasingly popular as surrogates for QM. Neural networks with many thousands of parameters excel in capturing structure within a large dataset, but may struggle to extrapolate beyond the scope of the available data. Here we present a highly automated active learning approach to iteratively collect new QM data that best resolves weaknesses in the existing ML model. We exemplify our approach by developing a general potential for elemental aluminum. At each active learning iteration, the method (1) trains an ANI-style neural network potential from the available data, (2) uses this potential to drive molecular dynamics simulations, and (3) collects new QM data whenever the neural network identifies an atomic configuration for which it cannot make a good prediction. All molecular dynamics simulations are initialized to a disordered configuration, and then driven according to randomized, time-varying temperatures. This nonequilibrium molecular dynamics forms a variety of crystalline and defected configurations. By training on all such automatically collected data, we produce ANI-Al, our new interatomic potential for aluminum. We demonstrate the remarkable transferability of ANI-Al by benchmarking against experimental data, e.g., the radial distribution function in melt, various properties of the stable face-centered cubic (FCC) crystal, and the coexistence curve between melt and FCC.'
author:
- 'Justin S. Smith'
- Benjamin Nebgen
- Nithin Mathew
- Jie Chen
- Nicholas Lubbers
- Leonid Burakovsky
- Sergei Tretiak
- Hai Ah Nam
- Timothy Germann
- Saryu Fensin
- Kipton Barros
bibliography:
- 'references.bib'
title: Automated discovery of a robust interatomic potential for aluminum
---

Introduction
============

Machine learning (ML) methods have become ubiquitous across many scientific disciplines. ML is capable of predicting physical phenomena at a high accuracy compared to reference data and with a relatively low computational cost, while also providing highly automated fitting procedures. ML-based models provide accurate and fast prediction of potential energies, forces, atomic charges and other properties of atomic systems [@Rupp2012; @Bleiziffer2018MachineCalculations; @Butler2018MachineScience; @Sifain2018]. Compared with classical potentials, the flexibility of ML potentials allows them to be fit to large and complex datasets, and suggests the possibility of unprecedented transferability. Applications to materials physics, chemistry, and biology are innumerable. To give some examples, simulations for crystal structure prediction, drug development, materials aging, and high strain/strain-rate deformation would all benefit tremendously from better interatomic potentials.

Machine learning (ML) of interatomic potentials is a rapidly advancing research topic. ML potentials have been fit to a variety of datasets including materials [@Behler2007; @Handley2014; @Szlachta2014; @Li2015; @Thompson2015SpectralPotentials; @Kruglov2017; @Huan2017; @Jindal2017; @Botu2017MachineOutlook; @Schutt2018; @Deringer2018ComputationalTheory; @Suwa2019MachineElectrons; @Pozdnyakov2019FastPotential; @Liu2020StructurePotential] and molecules [@Morawietz2012ACharges; @Smith2017; @Yao2017ThePhysics; @Lubbers2018; @Nguyen2018ComparisonExpansions; @Gastegger2017MachineSpectra; @Smith2019ApproachingLearning; @Schran2019AutomatedGround]. The datasets are calculated from computationally expensive *ab initio* quantum mechanics (QM) methods, most commonly density functional theory (DFT) [@Engel2011DensityTheory]. Trained on this data, the model can predict energy and forces for new atomic configurations. ML potentials should exactly encode translation, rotation, and permutation symmetries. Forces should also be calculated as the exact gradient of the predicted energy. Another important ingredient is an assumption of spatial locality: each atomic force only depends on neighboring atoms within a fixed radius, typically of order 5 to [${10}\,\mathrm{{\textup{\AA}}}$]{}. Long-range Coulomb interactions or dispersion corrections may be added [@Yao2017ThePhysics; @Unke2019PhysNet:Charges].

Compared to ML models, classical potentials have highly restrictive functional forms that might not be appropriate for the full space of atomic configurations. A common challenge, for example, is to design a single classical potential that captures transitions between multiple incompatible crystal phases. Over time, the community has identified a need for more flexible functional forms. For example, the embedded atom method (EAM) [@Daw1983SemiempiricalMetals], has lead to generalizations such as modified EAM (MEAM) [@Lee2010TheSimulations] and multistate MEAM [@Baskes2007MultistateMethod]. Such potentials are generally the tool of choice for molecular dynamics simulations of bulk metals but may struggle to generalize beyond the conditions to which they were trained. This shortcoming can be partly attributed to inflexibility of the assumed functional form. In contrast to classical potentials, neural network-based ML potentials may contain $\sim 10^5$ fitting parameters and make almost no assumptions beyond physical symmetries and spatial locality.

Such flexibility potentially allows ML potentials to achieve the accuracy and generality of the reference QM calculations, but with much greater speed. Because of spatial locality, the computational cost to calculate energy and forces scales linearly with system size. The Neurochem implementation [@2017ASE_ANI; @Smith2017] of our ANAKIN-ME (ANI) model, running on a single modern GPU (Nvidia RTX 2080 Ti), calculates all forces for a box of 1,000 atoms in about [${20}\,\mathrm{ms}$]{}. Classical potentials such as MEAM might be 1 to 2 orders of magnitude faster than well optimized ML potentials, whereas reference QM calculations are typically many orders of magnitude slower. As with classical potentials, one can use ML potentials to perform large-scale molecular dynamics simulations by employing spatial domain decomposition. Scalability is partly limited by the fact that, for typical ML codes, upwards of 1,000 atoms should be assigned to each GPU for greatest efficiency.

Building a transferable ML potential requires a large set of QM calculations [@Yao2017ThePhysics]. Prior work shows that ML potentials can be stunningly accurate over a wide swath of chemical and phase space [@Smith2017; @Yao2017ThePhysics; @Kobayashi2017NeuralAlloys; @Sifain2018; @Zhang2019ActiveSimulation; @Zhang2019DP-GEN:Models], but finding sufficiently diverse and reliable training data, to achieve true transferability, remains a challenge. For example, one could readily capture a limited range of physics (e.g. some hand-selected phases in equilibrium) by training to just that specific type of data. However, such success may be misleading, in that the model may fail to generalize beyond the scope of training data. More complex and challenging datasets will push ML models closer to capturing the true range of both equilibrium and *non-equilibrium* physics for a given material.

Active learning (AL) is an iterative dataset construction approach. New data is added to the training set whenever an estimate for model uncertainty exceeds a threshold. AL has recently become paramount in the construction of ML potentials [@Reker2015Active-learningDiscovery; @Podryabinkin2017ActivePotentials; @Gastegger2017MachineSpectra; @Gubaev2018MachineLearning; @Bernstein2019DeSurfaces]. Previous work has employed AL to drive nonequilibrium sampling of large datasets through organic chemical space, yielding the highly general ANI-1x potential [@Smith2018LessLearning]. Other recent research by Gubaev et al. [@Gubaev2019AcceleratingPotentials] has explored the use of AL with moment tensor potentials to construct atomistic potentials for materials. Zhang et al. also applied AL to materials using the deep potential model [@Zhang2019ActiveSimulation] for MgAl alloys. AL was used by Deringer, Pickard, and Csányi to build an accurate and general model for elemental Boron [@Deringer2018Data-DrivenBoron]. Although AL is rapidly growing into a standard technique, prior studies have commonly utilized substantial human guidance in selecting the physically relevant search space (e.g., by initializing sampling to the correct crystal structures). With such an approach, one may then question the extent to which the learned models generalize beyond the specific conditions targeted by the AL sampling procedure.

In the present study, we demonstrate that [*it is possible to obtain a broadly accurate potential for aluminum with essentially zero human guidance*]{}. Our AL sampling begins with completely disordered atomic configurations and automatically discovers physically relevant crystal phases, including Face Centered Cubic (FCC), Hexagonal Close Packed (HCP), and Body Centered Cubic (BCC), and defects. Our core assumptions are limited to the following: (1) The potential should satisfy physical symmetries (e.g., translation, rotation invariance), (2) The potential is local (no interactions beyond [${7}\,\mathrm{{\textup{\AA}}}$]{}), and (3) The AL sampler is loosely restricted in temperature (up to about [${2000}\,\mathrm{K}$]{}) and pressure (up to about [${50}\,\mathrm{GPa}$]{}). By avoiding any further use of human guidance, we aim to enable true materials discovery without *a priori* knowledge of the phases of interest. Scaled to hundreds of nodes on the Sierra supercomputer, each containing four Tesla V100 GPUs, our AL algorithm required about five days to automatically produce ANI-Al, the general potential for aluminum introduced in this work.

The AL algorithm uses the ML model under construction to drive molecular dynamics sampling for its data collection. All MD trajectories begin from a disordered (high temperature) configuration. By applying time-varying pressures and temperatures, this dynamics generates a variety of interesting atomic configurations. In particular, we can observe nucleation into various poly-crystalline structures, and the formation of crystal defects. Iteratively, new QM data is generated and added to the training dataset whenever the estimated uncertainty of the ML model exceeds a threshold. The ML model is periodically retrained to the new data. About 50 such iterations resulted in our final AL dataset. With this data, we trained models using two distinctly different ML architectures: ANI [@Smith2017] and HIP-NN [@Lubbers2018]. Performance of ANI and HIP-NN was comparable, demonstrating that the key enabler of an ML potential is the quality of the training data. Our ANI-Al potential exhibits state-of-the-art accuracy for elemental Al, as demonstrated by benchmarks both in equilibrium (e.g., properties of the stable FCC and melt phases, and map the coexistence temperature as a function of pressure) and highly out-of-equilibrium (e.g., theoretical energy barriers for the transformation between different crystal structures) conditions.

Building the ANI-A potential {#Methods}
============================

This section presents details of the automated procedure to build ANI-Al, our general purpose machine learning potential for bulk aluminum.

The ANI machine learning model {#meth:animodel}
------------------------------

ANI is a neural network architecture for modeling interatomic potentials. Our prior work with ANI has largely focused on modeling clusters of organic molecules [@Smith2019ApproachingLearning]. A variety of ANI potentials are available online in the ASE\_ANI github repository [@2017ASE_ANI]. Here we present ANI-Al, our ANI model for aluminum in both crystal and melt phases.

Our training data consists of DFT calculations, evaluated on “interesting” atomic configurations, as identified by the active learning procedure (Sec. \[meth:activelearn\]). Section 1.1 of the supplemental information describes our selected DFT/PBE parameters in detail. The input to ANI is an atomic configuration (nuclei positions and species). ANI first transforms the DFT data into a representation amenable to further processing. For this, ANI employs Behler and Parrinello [@Behler2007] type atomic descriptors, but with modified angular symmetry functions [@Smith2017]. Details of all model hyper-parameters are provided in the supplemental information Section 1.1. The output of ANI is a total energy prediction, computed as a sum over local contributions evaluated at each atom. Using backpropagation, one can efficiently calculate all forces as gradients of the predicted energy.

Each DFT calculation outputs the total system energy $E$ and the forces $\mathbf{f}_{j} = \mathrm{d} E / \mathrm{d} \mathbf{r}_j$ for all atoms $j=1\dots M$. The loss function, $$L = {\left(\hat{E}-E\right)^2+\frac{\ell_0}{M} \sum_{j=1}^{M}{\left(\hat{\mathbf f}_{j}-\mathbf f_{j}\right)^2}},$$ is a measure of disagreement between the ANI predictions for energy, $\hat{E}$, and forces, $\hat{\mathbf f}_{j} = \mathrm{d}\hat{E} / \mathrm{d} \mathbf{r}_j$, and the DFT reference data. Training ANI corresponds to optimizing the model parameters to minimize this loss, summed over all DFT calculations in the dataset. We select the characteristic length scale $\ell_0 = 0.01$ to balance energy and force training terms.

Standard optimization techniques require calculating the gradient of $L$ with respect to all model parameters $W$ (modern ML potentials might have $\sim 10^5$ trainable model parameters). Because forces $\hat{\mathbf f}_{j}$ appear in $L$, calculating $\mathrm{d} L / \mathrm{d} W$ involves *second* derivatives of the ANI energy output, i.e., $\mathrm{d}^2 \hat E / \mathrm{d} W \mathrm{d} \mathbf r_j$. Fortunately, direct calculation of all such second partial derivatives is not required. Frameworks such as TensorFlow or PyTorch support iterated backpropagation, thereby enabling efficient calculation of the full gradient $\mathrm{d} L / \mathrm{d} W$ [@Yao2017ThePhysics; @Lubbers2018; @Zhang2018DeepMechanics]. However, we train our models in Neurochem, a C++/CUDA implementation which is many times faster than comparable PyTorch implementations. The source code transformation required to perform backpropagation *twice* would be very difficult to perform manually within a hand-written CUDA/C++ code.

To train to force data, we have devised an efficient finite differencing scheme to approximate $\mathrm{d} L / \mathrm{d} W$. Our scheme involves just two first-order backpropagation calculations, $\mathrm{d} \hat{E}_i[\mathbf r^\pm] / \mathrm{d} W$, evaluated at specially displaced atomic positions, $\mathbf r \mapsto \mathbf r \pm \epsilon \delta \mathbf f$, with $\epsilon$ a small parameter and $\delta \mathbf f = \hat{\mathbf f} - \mathbf f$ the observed force error. Details will be presented in a future work.[^1]

To improve the quality of our predictions, a single ANI model actually employs ensemble-averaging over 8 neural networks. Each neural network in the ensemble is trained to the same data, but using an independent random initialization of the model parameters. We observe that ensemble-averaged energy and force errors can be up to 20% and 40% smaller, respectively, than those of a single neural network prediction.

Active learning {#meth:activelearn}
---------------

![Diagram of the active learning sampling algorithm employed in this work. Multiple such cycles can be run simultaneously, with occasional synchronization points to combine all new data into a single global dataset. The MD sampling, DFT data generation, and ML model training all benefit from GPU-acceleration.[]{data-label="fig:aldatagen"}](images/AL-datageneration.png){width="1.0\columnwidth"}

### Overview

The active learning process employed here is similar to that in previous work, [@Smith2018] adapted for materials problems and efficient parallel execution on hundreds to thousands of nodes on the Sierra supercomputer. We first train an initial model to a dataset of about 400 random disordered atomic configurations, generated as in Sec. \[meth:sampling:randomdisorder\]. Next, we begin the AL procedure. Using the current ML potential, we simulate many MD trajectories, each initialized to a random disordered configuration. During these simulations the temperature is varied to accelerate sampling of new configurations (Sec. \[meth:activelearn:tempschedule\]). As these simulations run, an uncertainty quantification metric [@Seung1992QueryCommittee; @Smith2018] is used to determine whether the model is operating as expected. If ML uncertainty exceeds a threshold, the MD trajectory is terminated and the final atomic configuration is placed on a queue for DFT calculation and addition to the training dataset. Periodically, the ML model is retrained to the updated training dataset. This AL loop is iterated until the cost of the MD simulations becomes prohibitively expensive. Specifically, we terminate the procedure when typical MD trajectories reach about [${250}\,\mathrm{ps}$]{} (about $2.5\times 10^5$ timesteps) without uncovering any weaknesses in the ML model. The final active learned dataset contains 6,352 DFT calculations, each containing 55 to 249 atoms, and having varying levels of disorder.

We emphasize that this active learning procedure is fully automated, and receives essentially no human guidance. For example, each MD simulation is initialized to a random disordered configuration, and the machine must *discover* the existence of crystal phases.

Section 1.2 of the SI describes the active learning procedure in more detail.

### Randomized atomic configurations {#meth:sampling:randomdisorder}

We employ randomized atomic configurations to collect an initial dataset of DFT calculations, and to initialize all MD simulations for AL sampling. The procedure to randomize a supercell is as follows:

1.  Randomly sample each of the three linear dimensions of the orthorhombic supercell uniformly from the range 10.5 to [${17.0}\,\mathrm{{\textup{\AA}}}$]{}

2.  Randomly select a target atomic density $\rho$ uniformly from the range 1.80 to [${4.05}\,\mathrm{g/cm^3}$]{}.

3.  Iteratively place atoms randomly in the supercell. If the proposed new atom lies within a distance $r_\textrm{min} = {{1.8}\,\mathrm{{\textup{\AA}}}}$ of an existing atom (i.e., roughly the van der Waals radius), that placement is rejected as unphysical. Placement of atoms is repeated until the target density $\rho$ has been reached.

### Nonequilibrium temperature schedule {#meth:activelearn:tempschedule}

To maximize the diversity in active learning sampling, we perform the MD simulations with a Langevin thermostat using a temperature that varies in time according to a randomized schedule.

Starting at time $t=0$, and running until $t=t_\textrm{max} = {{250}\,\mathrm{fs}}$, the applied temperature is, $$\label{eq:Tsch}
      T(t) = T_\textrm{start} + \frac{t}{t_\textrm{max}} (T_\textrm{end} - T_\textrm{start}) +  T_\textrm{mod} \sin^2 (\pi t / t_0) \\$$ The first two terms linearly ramp the background temperature. The initial temperature $T_\textrm{start}$ is randomly sampled from the range [${10}\,\mathrm{K}$]{} to [${1000}\,\mathrm{K}$]{}. The final background temperature $T_\textrm{end}$ is randomly sampled from the range [${10}\,\mathrm{K}$]{} to [${600}\,\mathrm{K}$]{}. The third term in Eq.  superimposes temperature oscillations. The modulation scale $T_\textrm{mod}$ is randomly sampled from the range [${0}\,\mathrm{K}$]{} to [${2000}\,\mathrm{K}$]{}. The oscillation period $t_0$ is randomly sampled from the range [${10}\,\mathrm{ps}$]{} to [${50}\,\mathrm{ps}$]{}.

By spawning MD simulations with many different temperature schedules, we hope to observe a wide variety of nonequilibrium processes. Given that each MD simulation begins from a disordered melt configuration, we hope that the nonequilibrium dynamics will automatically produce: (1) nucleation into various crystal structures (in particular, the ground-state FCC crystal), (2) a variety of defect structures and dynamics (dislocation glide, vacancy diffusion, etc.) and (3) rapid quenches into disordered glass phases. Acquiring snapshots from these types of dynamics will be crucial to the diversity of the training dataset and, thus, to the overall generality of the ANI-Al potential.

Accuracy Benchmarks
===================

Here we present a variety of benchmarks for ANI-Al, our machine learned potential for bulk aluminum. As described in Sec. \[Methods\], ANI-Al is trained from over 6,000 DFT calculations that were carefully selected using an iterative “active learning” procedure.

![Crystal energies (per atom) as a function of volume (per atom), relative to the ground state. Solid lines represent ANI-Al predictions and circles represent DFT reference calculations. Vertical bars represent sample variance of the eight neural networks comprising the (ensembled) ANI-Al model. Panel (b) is a magnification of panel (a) near the energy minima. The yellow region, from 11 to [${25}\,\mathrm{{\textup{\AA}}^3/atom}$]{}, indicates the approximate sampling range of the training data. Crystal structures explored: Diamond, simple cubic (SC), body centered cubic (BCC), hexagonal close packed (HCP), and face centered cubic (FCC).[]{data-label="fig:coldcurveAl"}](images/cold_curve_Al_AL_vert.pdf){width="\linewidth"}

Predicting crystal energies
---------------------------

Figure \[fig:coldcurveAl\] shows ANI-Al predicted energies (solid lines) for select crystal structures. ANI-Al correctly predicts that FCC is the lowest energy structure for aluminum. Vertical bars show the sample variance over the eight neural networks that comprise a single ANI-Al model (i.e., the uncertainty measure used within the active learning procedure). DFT reference data is shown in circles.

For both ANI-Al and DFT calculations, energies are measured relative to the FCC ground state. Let $\epsilon_x$ represent the error in the ANI-Al prediction, relative to DFT, for crystal structure $x$ at its energy-minimizing volume. By definition, we employ energy shifts such that $\epsilon_\textrm{fcc} = 0$. After FCC, the second lowest energy structure is HCP, for which the ANI-Al error is $\epsilon_\textrm{hcp} = 0.42\,\textrm{meV}/\textrm{atom}$. Note that FCC and HCP are competing close-packed structures, and both can reasonably be expected to emerge in our active learning dynamics. Specifically, HCP can be viewed, locally, as a stacking fault within FCC, and it is plausible that such configurations are contained within our training dataset. BCC, by contrast, is only physical in aluminum at much higher densities, far beyond the range of our active learning sampling. It is not surprising, therefore, that the ANI-Al error for BCC is an order of magnitude larger, $\epsilon_\textrm{bcc} = 5.3\,\textrm{meV}/\textrm{atom}$. Simple cubic and diamond crystals are less physical still, and we observe $\epsilon_\textrm{sc} = 37\,\textrm{meV}/\textrm{atom}$ and $\epsilon_\textrm{diamond} = - 44\,\textrm{meV}/\textrm{atom}$. Nonetheless, the agreement between ANI-Al and DFT observed in Fig. \[fig:coldcurveAl\] is remarkable. We emphasize that ANI-Al is not explicitly trained to any of these crystal structures. Also note that ANI-Al is capable of extrapolating beyond the range of densities (yellow region) contained in the training data. Similar observations were found in Ref. .

Predicting elastic constants
----------------------------

We can compare ANI-Al predicted elastic constants against experimental data. A particularly important property is the bulk modulus, which corresponds to the curvature of the FCC cold curve at its minimum (Fig. \[fig:coldcurveAl\]b). Experimentally, the FCC bulk modulus is measured to be [${79}\,\mathrm{GPa}$]{} [@Simmons1971SingleHandbook], whereas the ANI-Al prediction is [${77.3}\,\mathrm{GPa}$]{}. The full set of FCC elastic constants is measured experimentally to be, $C_{11}={{114}\,\mathrm{GPa}}$, $C_{12}={{61.9}\,\mathrm{GPa}}$, and $C_{44}={{31.6}\,\mathrm{GPa}}$ [@Simmons1971SingleHandbook]. For ANI-Al, we predict $C_{11}={{117}\,\mathrm{GPa}}$, $C_{12}={{57.2}\,\mathrm{GPa}}$, and $C_{44}={{30.4}\,\mathrm{GPa}}$.

There are several possible sources for the disagreement between ANI-Al and experiment: (1) Incomplete sampling of FCC crystal configurations, (2) failure of the ML model to capture the training data, (3) inconsistency of the DFT data due to lack of convergence in $k$-space basis set, and (4) inaccuracy of the DFT functional itself.

The largest discrepancy in the ANI-Al predicted elastic constants is observed in $C_{12}$, for which the relative error is 8%. This elastic constant measures the linear response of a normal stress to a strain in an orthogonal direction. For the ML model to precisely capture $C_{12}$, the training data should ideally contain bi-axially strained FCC configurations. The mechanisms by which our active learning can generate strained FCC are somewhat limited (e.g., by nucleating lattices that are incommensurate with the supercell, or by formation of frustrated polycrystals), and it is possible that enhancements to the sampling strategy could yield better prediction of elastic constants. For example, future work might employ time-varying applied strains of the supercell, in addition to the time-varying temperatures employed in the present study.

Lack of DFT convergence may also be a limiting factor. Our choice of a $3\times 3\times 3$ $k$-space basis was, in retrospect, too small, given that our supercells have linear dimensions up to [${17}\,\mathrm{{\textup{\AA}}}$]{}. Consequently, we observed an artificial dependence of calculated DFT properties on the size of the supercell. Note, however, that the errors in ANI-Al predicted elastic constants are roughly the same order as good DFT/PBE calculations. For example, prior theoretical work calculated $C_{44} = {{29.5}\,\mathrm{GPa}}$ [@Pham2011Finite-temperatureMeasurements], which is in line with the ANI-Al result, [${30.4}\,\mathrm{GPa}$]{}, both predictions being smaller than experiment, [${31.6}\,\mathrm{GPa}$]{}.

For the task of predicting elastic constants, ANI-Al accuracy is on par with many classical potentials (see Fig. S1 in the supplemental information). Whereas classical potentials are explicitly *designed* to reproduce experimental elastic constants, in ANI-Al this capability is an *emergent* property. Our active learning sampling discovers the FCC lattice and its properties entirely on its own.

![Transformational energy barriers. We compare ANI-Al and various classical potentials to reference DFT data. (a) Volume-conserving Bain path energies. (b) Generalized Stacking Fault Energy (GSFE) slip path. (c) GSFE twinning path.[]{data-label="fig:bain_sf"}](images/Bain-SF_paper.pdf){width="1.0\columnwidth"}

Predicting energy barriers in crystal
-------------------------------------

The Bain path represents a volume-preserving homogeneous deformation that transforms between FCC and BCC crystals. Starting from the initial FCC cell ($c/a$= 1), we compress along one of the $\langle 100 \rangle$ directions (length $c$) while expanding equally in the two orthogonal directions (lengths $a=b$). The special value of $c/a=$ 1/$\sqrt{2} \approx 0.71$ would correspond to BCC symmetry. Figure \[fig:bain\_sf\]a shows energies along this Bain path, in which $c/a$ varies continuously while conserving volume, $a^2 c$. We compare ANI-Al to DFT reference calculations, as well as seven EAM-based potentials [@Mendelev2008AnalysisCu; @Sjostrom2016MultiphaseTheory; @Liu2004AluminiumEnergy; @Mishin1999InteratomicCalculations; @Zope2003InteratomicSystem; @Winey2009AAluminum; @Pascuet2015AtomicAlloy; @Sheng2011HighlyMetals]. Figure S3 of the SI provides the errors for each method along these paths. The observed maximum at $c/a=$ 1/$\sqrt{2}$ indicates that the BCC structure is unstable to tetragonal deformation.

A stacking fault in FCC represents a planar defect in which the crystal locally is in HCP configuration within the nearest neighbor shell (note that FCC and HCP are competing close packed structures). The generalized stacking fault energy (GSFE) slip path provides an estimate of the resistance for dislocation slip and the energy per unit area required to form a single stacking fault. The GSFE twinning path (also known as the Generalized Planar Fault Energy) is an extension of the slip path and provides an estimate of the energy per unit area required to form $n$-layer faults (twins) by shearing $n$ successive $\{111\}$ layers along $\langle 112 \rangle$. We calculated the GSFE slip path and the twinning paths using standard methods [@Vitek1968IntrinsicCrystals; @Duesbery1998OverviewMetals; @Tadmor2004ACrack; @VanSwygenhoven2004StackingMetals].

Figures \[fig:bain\_sf\]b and \[fig:bain\_sf\]c show energies along the GSFE slip and twin paths, respectively. As before, we compare with seven EAM-based potentials. The ANI-Al potential agrees quite well with the reference DFT data for all measurements in Fig. \[fig:bain\_sf\]. To quantify this agreement, we calculate the root mean squared error (RMSE), formed as an average over the Bain, GFSE slip, and GFSE twinning paths. ANI-Al achieves RMSE values of [${4.5}\,\mathrm{meV/atom}$]{}, [${16.6}\,\mathrm{mJ/m^2}$]{}, and [${11.4}\,\mathrm{mJ/m^2}$]{}, respectively. For predicting these paths, the best classical potential is by Mishin et al.[@Mishin1999InteratomicCalculations], which achieves errors of [${4.3}\,\mathrm{meV/atom}$]{}, [${52.5}\,\mathrm{mJ/m^2}$]{}, and [${15.9}\,\mathrm{mJ/m^2}$]{}. Figure S2 of the SI provides RMSE and mean absolute error (MAE) values for each method compared to the DFT reference computed paths. It is interesting to note that the Winey et al. potential[@Winey2009AAluminum], which does exceptionally well in predicting many FCC properties (see Table S3 of the SI), struggles to accurately predict the Bain and GSFE slip paths. Errors in modeling the BCC and FCC energy barriers can have severe consequences in MD simulations. We will show an example in Sec. \[results:meltcurve\], where the Mendelev et al. potential[@Mendelev2008AnalysisCu] predicts transformation from FCC at BCC at just [${20}\,\mathrm{GPa}$]{}, whereas the physically correct pressure should be hundreds of GPa.

Predicting radial distribution functions
----------------------------------------

![image](images/liq_rdf-rho.pdf){width="\linewidth"}

To validate our ANI-Al model in the liquid phase, we carry out MD simulations to measure radial distribution functions (RDF), densities at various temperatures, and a partial solid/liquid phase diagram. Figure \[fig:rdfdensity\]a compares simulated RDFs with experimental measurements [@Mauro2011HighAl] at 1123K, 1183K, and 1273K. Independent simulations were performed in the isobaric-isothermal (NPT) ensemble to determine equilibrium densities of liquid Al at the relevant (P,T) conditions. The starting configurations for Radial Distribution Function (RDF) calculations were created using these densities and was equilibrated for 50 ps in the NVT ensemble using the Nosé-Hoover-style equations of motion [@Hoover1985CanonicalDistributions] derived by Shinoda et al [@Shinoda2004RapidStress]. Reported RDFs were calculated (bin size of 0.05 Å) by averaging 100 instantaneous RDFs, which were 0.1 ps apart, in the final 10 ps of the NVT equilibration. A timestep of 1 fs was used for these simulations. We find ANI-Al melt simulations to be most accurate near the melting temperature of $T_\textrm{melt} \approx {{933}\,\mathrm{K}}$. At higher temperatures, however, deviations become visibly apparent. This increased error may stem from the fact that the physically correct potential should actually be *temperature-dependent*, whereas our DFT training data was entirely calculated at zero electronic temperature. Figure \[fig:rdfdensity\]b compares ANI-Al predicted densities at various temperatures (still at atmospheric pressure) to multiple experimental values [@Assael2006ReferenceIron; @Smith1999MeasurementTechnique; @YatsenkoSPKononenkoVI1972ExperimentalGa; @CoyW.J.Mateer1955CoyMateer; @LevinE.S.AyushinaG.D.Geld1968LevinDensity; @GebhardtE.BeckerM.andDorner1955GebhardtAl]. All temperatures are above the melting point. The agreement between ANI-Al predictions and experiment is comparable to the variation between different experiments.

Predicting liquid-solid phase boundaries {#results:meltcurve}
----------------------------------------

![(a) Aluminum melt curves calculated from DFT [@Sjostrom2016MultiphaseTheory], ANI-Al, and the Mendelev et al. EAM potential [@Mendelev2008AnalysisCu], compared with experimental data [@Hanstrom2000HighAluminium]. Below [${120}\,\mathrm{GPa}$]{} we show FCC-liquid coexistence. Above [${240}\,\mathrm{GPa}$]{} we show BCC-liquid coexistence. The inset zooms to pressures from 0 to [${20}\,\mathrm{GPa}$]{}. (b) Errors in predicting the melt temperature at atmospheric pressure.[]{data-label="fig:melt_curve"}](images/melt_curve_melt_point.pdf){width="\linewidth"}

Figure \[fig:melt\_curve\] shows the liquid-solid coexistence line in the pressure-temperature plane. At each pressure, we calculated the coexistence temperature by performing simulations with an explicit solid-liquid interface [@Morris1994MeltingPhases; @Morris2002TheSimulations; @Espinosa2013OnModel]. The details of these simulations are provided in Section 1.4.1 of the supplemental information. We studied FCC-liquid coexistence at pressures up to about [${100}\,\mathrm{GPa}$]{} to match the range of experimental data [@Hanstrom2000HighAluminium]. The ANI-Al simulations are compared with prior DFT calculations [@Sjostrom2016MultiphaseTheory] and a classical MD potential. For the latter, we used the Mendelev et al. potential [@Mendelev2008AnalysisCu], which was explicitly parameterized to model the melting point of aluminum, $T_\textrm{melt} \approx {{933}\,\mathrm{K}}$ at atmospheric pressure. At this pressure, both Mendelev and ANI-Al potentials predict an FCC melting point of about [${925}\,\mathrm{K}$]{}, in good agreement with experiment. The ANI-Al model accurately predicts the melt curve up to [${50}\,\mathrm{GPa}$]{} and [${3000}\,\mathrm{K}$]{}, where it begins to underestimate the melting temperature. Figure S2 of the SI shows the distribution of pressures in the active learning training dataset. The Mendelev model begins to underestimate the melting temperature at around [${5}\,\mathrm{GPa}$]{}. Surprisingly, at pressures above about [${20}\,\mathrm{GPa}$]{}, the Mendelev simulations are unstable to nucleation into the BCC crystal. According to prior DFT-based studies [@Tambe2008BulkStudy; @Sjostrom2016MultiphaseTheory], and experiment [@Fiquet2019StructuralGPa], the solid-to-solid transition out of FCC should require hundreds of GPa.

At higher pressures, BCC becomes the preferred crystal phase. Figure \[fig:melt\_curve\] shows the liquid-BCC coexistence curve pressures around [${250}\,\mathrm{GPa}$]{}. The ANI-Al predictions are within 5% of prior DFT calculations [@Sjostrom2016MultiphaseTheory]. Interestingly, in performing these simulations, the ANI-Al potential is extrapolating well beyond its training data. During active learning, we collected DFT training data only over a limited region of atomic densities (yellow region of Fig. \[fig:coldcurveAl\]a) reaching up to about [${60}\,\mathrm{GPa}$]{} (Fig. S1 of the SI).

Phase transition dynamics {#results:phasechange}
-------------------------

![image](images/melt-freeze_MD.pdf){width="100.00000%"}

Next we carry out a nonequilibrium MD simulation to observe both freezing and melting dynamics. Our intent is to validate the ANI-Al predicted energies and forces at snapshots along the dynamical trajectory. Along the trajectory the temperature is slowly increased from [${300}\,\mathrm{K}$]{} to [${1500}\,\mathrm{K}$]{}, then cooled back to [${300}\,\mathrm{K}$]{}. The details of these simulations are provided in Section 1.4.2 of the supplemental information.

Figure \[fig:almeltfreeze\] shows the potential energy, mean force magnitude, and pressure for both ANI-Al and DFT along this trajectory. Melting from FCC to liquid occurs at around [${300}\,\mathrm{ps}$]{} and freezing occurs around [${700}\,\mathrm{ps}$]{}. Pressure was calculated using the method of Ref. . The inset images in the middle panel of Figure \[fig:almeltfreeze\] show the composition of the system before and after melting, and after refreezing. Compositional information was obtained using the Common Neighbor Analysis as implemented in the OVITO visualization software [@Stukowski2010VisualizationTool].

Every [${2.5}\,\mathrm{ps}$]{} along the trajectory we sampled a frame to perform reference DFT calculations. The error between ANI-Al and DFT is generally small. Over the full trajectory, the MAE for energy and mean force magnitude are just [${0.84}\,\mathrm{meV/atom}$]{} and [${4.1}\,\mathrm{meV/{\textup{\AA}}}$]{}, respectively. The MAE for ANI-Al predicted pressure is [${0.36}\,\mathrm{GPa}$]{}. Interestingly, there is a systematic tendency for ANI-Al to overestimate pressure, especially at negative pressures. This could reflect the fact that the majority of ANI-Al training data was sampled mostly at positive pressures (see Fig. S1 in the SI). Another possible source of error is Pulay stress due to incompleteness of our $3 \times 3 \times 3$ $k$-space basis for the DFT calculations.

Figures S4 and S5 of the SI further validate the ANI-Al model over large a range of temperatures and densities.

Importance of active learning
=============================

The success of ANI-Al hinges on the diversity of the active learned dataset. To demonstrate this, we compare ANI-Al against an ML model trained on a much more limited dataset. We will call this baseline dataset “FCC/Melt,” as it consists only of samples from the FCC and liquid phases. Specifically, the FCC/Melt dataset is constructed by taking regular snapshots from near-equilibrium MD trajectories. For each snapshot, we perform a DFT calculation to determine the reference energy and forces.

The first such MD trajectory is shown in Fig. \[fig:almeltfreeze\]. There, 108 atoms were initialized to FCC, heated from [${300}\,\mathrm{K}$]{} to [${1500}\,\mathrm{K}$]{}, and cooled back to [${300}\,\mathrm{K}$]{}. We take 300 snapshots from this trajectory, equally spaced in time, to add to the FCC/Melt dataset. For increased variety, the FCC/Melt dataset contains an additional 250 DFT calculations taken from the liquid phase over a range of temperatures and pressures (Sec. 1.4.3 in the SI contains details). In sum, the FCC/Melt dataset contains 550 DFT calculations for near-equilibrium FCC and liquid configurations.

[l l l]{} Model type & FCC/Melt test & AL test\
&\
FCC/Melt trained & 2.0/4.0 & 40/110

------------------------------------------------------------------------

\
AL trained & 1.4/1.9 & 1.3/1.9\
&\
FCC/Melt trained & 0.04/0.07 & 0.49/1.53\
AL trained & 0.03/0.04 & 0.04/0.06\
\[1ex\]

\[table:disordervsfccmelt\]

Table \[table:disordervsfccmelt\] compares our ANI-Al model, trained on the full active learned (AL) dataset, to an ANI model trained on the much more restricted FCC/Melt dataset. The two model types are compared by testing on held out portions of both datasets. Figures S6 and S7 in the SI show the associated correlation plots.

A conclusion of Table \[table:disordervsfccmelt\] is that both the AL trained and FCC/Melt trained models have comparable errors when predicting on the held out FCC/Melt test data. However, when testing on the held out AL data, the FCC/Melt trained model does quite poorly. This failure casts doubt on the ability of the FCC/Melt trained model to study *any* new dynamical physical process: Will a rare event occur that pushes the FCC/Melt trained model outside its range of validity? To mitigate this danger it is essential to make the training dataset as broad as possible, which is our aim with active learning.

Outlook
=======

We showed that machine learning enables automated construction of general purpose interatomic potentials. Such ML-based potentials may significantly outperform “classical” potentials in their accuracy, often achieving errors of just a couple meV per atom, as benchmarked over a wide variety of ordered and disordered atomic configurations of aluminium.

Neural networks provide a highly flexible and computationally efficient fitting function, suitable for training on extremely large *ab initio* datasets. For example, our optimized C++/CUDA implementation of ANI, running on a modern GPU, can calculate all forces for a box of 1,000 atoms in about [${20}\,\mathrm{ms}$]{}. This is many orders of magnitude faster than applying *ab initio* methods such as DFT to metallic systems.

Because ML models are so flexible, the quality and diversity of the training dataset is crucial. Counterintuitively, providing too much human guidance to the data collection procedure may actually be detrimental to the design of accurate and general models. A human-guided approach tends to bias models towards things we already know, leading to blind spots when pushed to unknown (but physically relevant) regimes of configuration space. For example, in prior work, datasets have often been designed to include specific crystal structures, and specific defect configurations. But in real-world MD simulations of solids, we often care about conditions beyond isolated defects. Defect *dynamics* are essential, and these may involve the complex collective motion of many atoms (consider, e.g., driven interactions between dislocation cores). The details of such processes seem extremely difficult to anticipate, but must somehow be captured in the dataset; if the dataset is lacking, then the ML model will often fail to generalize correctly.

In this work, we have developed an active learning framework that allows building ML based potentials in a fully automated way. The required human input is quite limited (apart from ML hyperparameters, we must select physical parameters such as the temperature and density ranges over which to sample, and the interaction cutoff distance). Using this framework, we automatically collected a dataset of DFT calculations for bulk aluminum, and used it to train the ANI-Al potential. Active learning automatically *discovered* that FCC is the ground state of aluminum, and learned to predict other crystal energies as well (Fig. \[fig:coldcurveAl\]). Furthermore, the ANI-Al model accurately predicts energy barriers for transforming between different crystal structures (Fig. \[fig:bain\_sf\]) and melt curves. (Fig. \[fig:melt\_curve\]).

A challenge for active learning of interatomic potentials is the tremendous demand on computational resources. Our final active learned dataset contains over 6,000 DFT calculations; each calculation was performed on a (typically disordered) supercell containing up to 250 atoms. For future work, it would be interesting to explore whether smaller supercells could be used. We also plan to employ *ab initio* calculations with greater accuracy, and this will make the accessible system sizes much smaller.

The sampling procedure used to find new atomic configurations is another significant expense. During active learning, MD simulations must run until a configuration arises for which the ML uncertainty exceeds a threshold. As the ML models improve, the required MD trajectories become longer (reaching hundreds of picoseconds) and eventually grow to become the dominant computational cost. It will be very useful to develop better sampling methods that quickly identify physically relevant atomic configurations for which the ML model is uncertain or to apply methods that automatically diversify the sampling space of interest [@Karabin2020AnPotentials]. These approaches can be extended to applications of AL guided ML in the general area of materials and chemical processes.

This work was partially supported by the LANL Laboratory Directed Research and Development (LDRD) and the Advanced Simulation and Computing Program (ASC) programs. Within ASC, we acknowledge support from the Physics and Engineering Modeling (ASC-PEM) subprogram and the Advanced Technology Development and Mitigation (ASC-ATDM) subprogram. We acknowledge computer time on the Sierra supercomputing cluster at LLNL, Institutional Computing at LANL, and the CCS-7 Darwin cluster at LANL. JSS was supported by the Nicholas C. Metropolis Postdoctoral Fellowship. NM and JSS were partially supported by the Center for Nonlinear Studies (CNLS). This work was performed, in part, at the Center for Integrated Nanotechnologies, an Office of Science User Facility operated for the U.S. Department of Energy (DOE) Office of Science.

[^1]: J. S. Smith, N. Lubbers, and K. Barros, in preparation
