Active Learning of Uniformly Accurate Inter-atomic Potentials for Materials Simulation

Linfeng Zhang, De-Ye Lin, Han Wang, Roberto Car, Weinan E

I Introduction

The inter-atomic potential energy surface (PES) plays a central role in the molecular modeling of materials. Obtaining an accurate and efficient representation of the PES is a central issue in molecular simulation. In this context, one faces the dilemma that ab initio methods are accurate but highly inefficient, while empirical force fields (FFs) are efficient, but there is a limited guarantee for their accuracy. Thus, there is a great demand for an efficient and uniformly accurate PES model that can be used to compute a broad range of atomistic properties for most material compounds of practical interest.

Developing empirical FFs has been challenging due to the high dimensionality and many-body character of the PES. Usually, empirical FFs parameterize the PES by assuming an analytical functional form in terms of relatively simple functions based on physical/chemical intuition, and by fitting the model parameters against a bundle of experimental properties and/or microscopic quantities from ab initio calculations. Some popular examples are the Lennard-Jones potential jones1924LJ, the Stillinger-Weber potential stillinger1985SW, the embedded-atom method (EAM) potential daw1984EAM, the CHARMM mackerell1998charmm/AMBER wang2000amber FFs, the reactive FFs van2001rff, etc. Representability and transferability are two main issues faced by empirical FFs. By representability, we mean the ability of the assumed functional form to reproduce accurately the target properties. By transferability, we mean the ability of a PES model to describe properties that do not belong to the set of fitting targets. Due to the physical/chemical knowledge encoded in the functional form, we expect the empirical FFs to be qualitatively transferable to a moderate range of thermodynamic conditions beyond those adopted for the fitting. However, as a consequence of assuming relatively simple functional forms, empirical FFs usually face a severe representability problem. Moreover, a substantial human effort in tuning the model parameters is often required to achieve the best balance in fitting the target properties.

Recent progress with machine learning (ML) methods is changing the outlook behler2007generalized; bartok2010gaussian; rupp2012fast; montavon2013machine; botu2016machine; chmiela2017machine; schutt2017schnet; bartok2017machine; smith2017ani; han2017deep; chen2018atomic; zhang2018deep; zhang2018end. ML models, being capable of learning complex and highly nonlinear functional dependence, are excellent in their representability. It is now possible, using modern ML approaches, to parametrize the PES using data from ab initio calculations to obtain models that have ab initio accuracy and are, at the same time, competitive regarding efficiency against empirical FFs. In spite of the remarkable success of these ML methods, there is no guarantee for the quality of ML models when they are used to predict the properties of a configuration that is far from the training dataset lecun2015deep. In addition, since the training data is usually generated with expensive first-principle calculations, one would like to obtain good ML models without having to rely on very large ab initio datasets. These questions arise not only for PES modeling, but in many other contexts when ML methods are applied to problems involving physical models.

To address this issue, we get inspiration from active learning settles2012active; rubens2015active, an area of supervised learning whose aim is to learn general purpose models with a minimal number of training data. A training data point involves an input and an output. For example, in an image recognition task whose goal is to judge whether a cat is in an image or not, the input is an array of digits that represents the image, and the output is a boolean proposition. Usually the output is called a label and the term labeling is used to denote the creation of a label. In the context of active learning, one typically faces a situation in which unlabeled data are abundant, but labeling is expensive. Therefore, an interactive algorithm is required to efficiently explore unlabeled data, collect feedbacks on-the-fly, and actively query the teacher for labels on data points with negative feedbacks. Along this line of thinking, at an abstract level, one can formulate an active learning procedure for PES modeling that involves three steps: exploration, labeling, and training.

Exploration requires an efficient sampler and an informative indicator. The sampler uses the current PES model to quickly explore the configuration space. The indicator monitors on-the-fly the configurations explored by the sampler, selects those with low prediction accuracy, and sends them to the labeling step.

Labeling means generating reference ab initio energies and forces for the selected configurations. Labeling can be done by a code that implements high-level quantum chemistry, quantum Monte Carlo, or density functional theory (DFT) methods. The labeled configurations are then added to the existing dataset and used in the new iteration for training.

Training requires a good model, or PES representation, which can fit the ever-increasing dataset with satisfactory accuracy. Such a representation should be efficient and should satisfy certain physical constraints like the extensive and symmetry-preserving properties of the PES.

The whole scheme falls into a closed loop: One starts with a relatively poor approximation of the PES and uses it to explore different configurations. Then a selected set of new configurations is labeled, and a new approximation of the PES is obtained by training. These three steps are repeated until convergence is achieved, i.e., the configuration space has been explored sufficiently, and a minimal set of data points have been accurately labeled. At the end of this procedure, a uniformly accurate PES model is generated.

In this work, our first goal is to translate the general proposal described above into a practical scheme for modeling the PES. In this scheme, for the PES representation, we use an advanced version of the Deep Potential (DP) model zhang2018end, which has shown great promise in learning the PES of a broad range of systems, such as insulators, molecular crystals, and a 5-component high entropy alloy, etc. See, e.g., Fig. 1 of Ref. zhang2018end. For the sampler, we use molecular dynamics (MD) based on the DP model. Thereafter DP based MD will be referred to as DPMD. At the same time, we introduce an indicator that we call the model deviation. This is done as follows. We train an ensemble of DP models using the same dataset but different initialization of the DP parameters. For each new configuration that is explored by DPMD, these models generate an ensemble of predictions. For each configuration, the model deviation is defined as the maximum standard deviation of the predicted atomic forces. A high model deviation indicates low quality in the model prediction and is proposed for labeling. In this work, we use in the labeling stage DFT within the generalized gradient approximation perdew1996generalized; kresse1996efficient; kresse1996efficiency; monkhorst1976special, which works well in the chosen testing examples. We will see that sampling is much cheaper than labeling, and only a very small fraction of the explored configurations is selected for labeling. We call the methodology introduced here the Deep Potential Generator, abbreviated DP-GEN.

Our second goal is to demonstrate the uniform accuracy of a PES model obtained in this way. To this end, we consider the example of Al, Mg, as well as Al-Mg alloys. Using DP-GEN, we construct a model that can accurately describe these systems at different compositions and thermodynamic conditions. The resulting PES model is evaluated from the point of view of a material scientist. We calculate several statical, dynamical, and mechanical properties, such as radial distribution functions (RDF), phonon spectra, elastic constants, etc. Some of these properties are compared with DFT results. We also compare DP calculated properties directly with experimental results when these are available. To further test the quality of the PES model, we introduce an automatic procedure based on the Materials Project (MP) database jain2013commentary. In this procedure, one searches the database by entering a material composition, such as Al-Mg in the present case. The database will then return a large number of locally stable structures, including many structures of potential practical interest. Based on these structures, we evaluate several equilibrium properties and compare the DFT predictions with those of the PES model. In addition, for each one of these structures, we automatically generate unrelaxed vacancy and interstitial defects as well as the set of surfaces corresponding to a range of Miller indices ong2013python. We then compare the relaxed formation energies of the defects and the unrelaxed formation energies of the surfaces predicted by DFT and by the PES model. We stress that these structures, i.e., crystals, defects, and surfaces, were not explicitly included in the training data. We find that our PES model can achieve uniform accuracy in the prediction of all of these structural properties.

We notice that there is a difference between active learning in conventional ML problems and the active learning we pursue here. This difference lies in exploration or sampling. Conventional active learning problems in ML typically deal with an existing unlabeled dataset. Here our dataset is generated on the fly via sampling. This means that we need to have an efficient sampling method.

We should mention that related work can be found in the literature podryabinkin2017active; smith2018less; herr2018metadynamics; bonati2018silicon; musil2019fast. In particular, Smith et al smith2018less utilized an active learning scheme to model the PES of organic molecules based on an existing large database smith2017ani-data. Moreover, Bartok et al bartok2018machine constructed a kernel based general purpose PES model for pure silicon, wherein they exhaustively enumerated possible structures for labeling. Finally, the principle of active learning was also used in the reinforced dynamics scheme zhang2018reinforced for enhanced sampling and free energy calculation.

II Methodology

In this section, we introduce the three essential components of the DP-GEN scheme: the model, the sampler, and the indicator. Fig. 1 shows a schematics of DP-GEN. To initialize the procedure, we label a small set of initial structures introduced in Fig. 1(a) and train an ensemble of preliminary DP models. More details on the simulation protocol and the iterative process are reported in the supplementary materials (SM).

Model. The DP scheme assumes that the potential energy $E$ can be written as a sum of atomic energies, i.e., $E=\sum_{i}E^{i}$ . Each atomic energy $E^{i}$ is a function of $\mathcal{R}^{i}$ , the local environment of atom $i$ in terms of the relative coordinates of its neighbors within a cut-off radius $r_{c}$ . The dependence of $E^{i}$ on $\mathcal{R}^{i}$ embodies the nonlinear and many-body character of the inter-atomic interactions. Therefore, we use a deep neural network function (DNN) to parameterize it, i.e., $E^{i}=E^{\bm{w}_{\alpha_{i}}}(\mathcal{R}^{i})$ . Here $\alpha_{i}$ indicates the chemical species of the $i$ -th atom; $\bm{w}_{\alpha_{i}}$ denotes the parameters of the DNN we call network parameters, that are determined by the training procedure. A vital component of the DP model is a general procedure that encodes $\mathcal{R}^{i}$ into the so-called feature matrix $\mathcal{D}^{i}$ . This procedure guarantees the conservation of the translational, rotational, and permutational symmetries of the system, without losing coordinate information in the local environment. Derivatives of the energy with respect to the atomic positions give the forces. During the training process, the network parameters evolve in order to minimize the loss function, a measure of the error in the energies and the forces predicted by DP relative to the labels, i.e., the corresponding DFT predictions kingma2015adam. Upon convergence, the model can match the labels within a small error tolerance. The details of the architecture of the DP model and the training process are given in Ref. zhang2018end.

Sampler. The goal of the sampler is to explore the configuration space in a range of thermodynamic variables, say temperature and pressure. Ideally one should develop an automatic/adaptive procedure for this purpose. However, since exploration is relatively cheap compared to labeling, we adopt a more heuristic approach in which the exploration is done through: (1) carefully selecting the initial configurations, and (2) exploring the volume-temperature space. We use a variety of crystal structures as our initial configuration, as in the procedure illustrated in Fig. 1(a). To explore the volume-temperature phase space, we adopt a temperature increasing scheme, in which the temperature of the DPMD simulations is increased systematically with the iteration index in the range 50-2000 K. We notice that many structures constructed in this way are far from equilibrium structures so that the subsequent DPMD simulations in the 50-2000K temperature range produce a large sample of configurations that may differ substantially from the initial structure. More details on the initial structures and the thermodynamic conditions in each iteration are summarized in Tables S1–S4.

Indicator. It is well-known that neural network models are highly nonlinear functions of the network parameters $\bm{w}_{\alpha_{i}}$ . The loss function, as a function of $\bm{w}_{\alpha_{i}}$ , is highly non-convex, i.e., several local minima exist in the landscape of the loss function. In the current work, we initialize the $\bm{w}_{\alpha_{i}}$ randomly according to the standard normal distribution. As a result, different initializations often lead to different minimizers of the loss function. These minimizers fit well the training data, so in the configurational region belonging to the neighborhood of the training data, they generate equally accurate energies and forces and show small deviations in their predictions. However, for snapshots “far” from the training data, these minimizers usually predict inaccurate values that show significantly larger deviation. This property of neural network models motivates us to define the indicator as the deviation of the predictions generated by an ensemble of DP models trained with the same dataset but with different parameter initializations. In practice, we define the model deviation, denoted as $\mathcal{E}$ , as the maximum standard deviation of the predictions for the atomic forces, i.e.:

where $i$ runs through the atomic indices in a configuration, and the ensemble average $\langle\cdots\rangle$ is taken over the ensemble of models. We find that using the predicted forces to evaluate the model deviation is generally better than using the predicted energies. The force is an atomic property and is sensitive to a failure in local predictions, while the energy is a global quantity and does not seem to provide sufficient resolution in this regard. Moreover, we find that a failure in local predictions can be better signaled by using the maximum over $i$ in Eqn. 1, instead of the average over $i$ ( $\frac{1}{N}\sum_{i}$ ).

III Results

As examples, we report the results of the DP-GEN scheme for Al, Mg and their alloys. At the end of the DP-GEN scheme, we collect a set of labeled data and obtain a DP model for the Al-Mg system. As shown in Table S4, about 650 million configurations were explored by DPMD, but only 0.0044% of them were selected for labeling. To get an idea of the usefulness of the resulting DP model for materials science applications, we compare the accuracy of the DP model in predicting important material properties with a state-of-the-art empirical FF like the modified embedded atom method (MEAM) baskes1992modified. MEAM adopts a more general definition of embedding than EAM in order to improve the description of directional bonding and of alloy systems. In this work, we compare our method with a very recent version of the Al-Mg MEAM potential that is available in the literature jelinek2012modified. We used DeePMD-kit wang2018kit in the training step, LAMMPS plimpton1995lammps in the exploration step, and VASP kresse1996efficiency; kresse1996efficient in the labeling step.

The equilibrium properties of pure Al are presented in Table 1, including the atomization energy and equilibrium volume per atom, defect formation energies, elastic constants and moduli, stacking fault energies, melting point, enthalpy of fusion, and diffusion coefficient. The defect formation energy is defined as $E_{\textrm{df}}=E_{\textrm{d}}(N_{\textrm{d}})-{N_{\textrm{d}}}E_{0},$ $d=v(i)$ indicating vacancy (interstitial) defects. $E_{\textrm{d}}$ denotes the relaxed energy of a defective structure with $N_{\textrm{d}}$ atoms and $E_{0}$ denotes the energy per atom of the corresponding ideal crystal at $T=0$ K. To compute the defect formation energies, we use a supercell in which we replicate $7\times 7\times 7$ times the primitive FCC cell. We estimate the melting temperature ( $T_{m}$ ) by simulating with DPMD coexisting crystal and liquid phases in a supercell containing 8000 atoms within the isothermal-isobaric ensemble at standard pressure. To estimate the liquid diffusion coefficient ( $D$ ), we perform DPMD simulations on large supercells (6912 atoms) for which finite size effects are negligible. For all the properties in Table 1, the DP predictions are in satisfactory agreement with DFT and/or experiment. Notice that MEAM reproduces quite accurately the solid state properties in Table 1, particularly when compared to experiment, which is not surprising since the basic experimental solid state properties have been used to tune the parameters of this FF. However, the vibrational properties at short wavelength, particularly the zone boundary phonons, are not reproduced well by MEAM in contrast to DP, as shown in Fig. 2. MEAM fails even more dramatically in predicting the properties of the liquid: the MEAM liquid is largely overstructured (see Fig. 3). Its diffusion coefficient is one order of magnitude smaller than in experiment or DP, and its enthalpy of fusion is also significantly smaller than in experiment or DP (see Table 1).

DFT, DP, and MEAM predictions for the equation of state (EOS) of Al are reported in Fig. 4. DP reproduces well the DFT results for all the crystalline structures considered here, i.e., FCC, HCP, double-hexagonal-closed-packed (DHCP), body-centered-cubic (BCC), SC and diamond. Interestingly, the range of DP accuracy extends well beyond the volume interval that was included in the training data, which is indicated by the yellow shaded area in the figure. As shown in the inset of Fig. 4, the energy difference between FCC and DHCP, and the one between DHCP and HCP is small, only 12 meV/atom and 19 meV/atom, respectively, yet DP reproduces accurately the relative stabilities. The MEAM potential performs well for FCC, HCP, DHCP, and SC, but shows significant deviations from DFT for diamond and BCC. DP and MEAM predictions for the phonon dispersion relations are compared with experimental results in Fig. 2. DP results agree very well with experiment.

The promise of ML potential models is to retain the accuracy of ab initio molecular dynamics (AIMD) at the cost of FF simulations. Therefore, ML potential models can be used to simulate much larger systems for much longer times than possible with AIMD. This is illustrated by our calculations for the diffusion coefficient and the radial distribution function (RDF) of the liquid, which were performed on large cells with 4000 atoms with very modest computational resources when using DP. Thus, the DP model opens opportunities for extending the power of ab initio methods.

The DP method gives similarly good results for the corresponding properties of pure Mg, which are reported in the SM.

Finally, we examine the surface formation energy $E_{\textrm{sf}}((lmn))$ , which describes the energy needed to create a surface with Miller indices $(lmn)$ for a given crystal, and is defined by $E_{\textrm{sf}}((lmn))=\frac{1}{2A}(E_{\textrm{s}}((lmn))-N_{\textrm{s}}E_{0}).$ Here $E_{\textrm{s}}((lmn))$ and $N_{\textrm{s}}$ denote the energy and number of atoms of the relaxed surface structure with Miller indices $(lmn)$ . $A$ denotes the surface area. We enumerate all the non-equivalent surfaces corresponding to Miller index values smaller than 4 for Al, and smaller than 3 for Mg. As shown in Fig. 5, the surfaces formation energies predicted by DP are close to DFT tran2016surface, and those predicted by MEAM are worse in all cases. We report in detail the values of surface formation energies for Al and Mg in Tables S6 and S7, respectively.

III.2 Mg-Al Alloys

The formation energy of an Mg-Al alloy system is defined as

In almost all tested cases, we observe an overall satisfactory agreement between DP predictions and DFT reference results. The accuracy of DP is significantly better than that of MEAM. We stress that the DP-GEN procedure is blind to the alloy structures used to compute the properties reported in Fig. 6, because these structures were not explicitly included in the training data. The number of atoms in the unit cell of 6 MP structures is larger than 32, which was the maximum number of atoms in the unit cell of the structures belonging to the training dataset. This suggests that in the case of Mg-Al alloys the DP model trained with relatively small periodic structures can, to some extent, be used to predict the properties of larger structures. Some structures tested have little in common with the initial training data. Yet the DP model produced satisfactory results, suggesting that it could work for a broader range of materials.

IV Summary and Outlook

The DP-GEN scheme is general, practical, and fairly automatic. To generate the DP model for the Al-Mg system, we did not use any existing DFT database (the MP database was only used for testing), nor did we use an exhaustive list of possible structures based on physical and chemical considerations. Instead, we explored the space of configurations using computationally efficient DPMD simulations. DFT calculations were only performed on a small subset of the configurations that showed large model deviation. This made possible to progressively improve the DP model.

The DP-GEN scheme is quite flexible. The three components, training, exploration, and labeling, are highly modularized and can be implemented separately and then recombined. This makes it easy to incorporate additional functionalities. For example, enhanced sampling techniques bonati2018silicon or genetic algorithms hajinazar2017stratified can be incorporated with minimal effort in the exploration module. We expect that the modular structure of DP-GEN should make possible to use this method to generate models for a variety of important problems, such as finding transition pathways for structural transformations and chemical reactions. The outcome of DP-GEN include the model and the accumulated data, which could be used for further applications. For example, if a rare-earth species is added to the Al-Mg system, one does not need to start the DP-GEN scheme from scratch. Instead, one could restart the DP-GEN scheme with the current model and data, and continue with the exploration of the configuration space involving the new species.

Besides alloys dominated by metallic bonding, it would be interesting to use the DP-GEN scheme to study other materials, such as ceramics, polymers, etc., which include different types of bond interactions. This should be possible because the applicability of DP-GEN relies on three main points: the representability of the model, the validity of the indicator, and the capability of the sampler. Several investigations suggest that the first two issues should be relatively independent of the details of the microscopic interactions. Indeed, our earlier studies zhang2018deep; zhang2018end indicate that the DP model can represent equally well the PES of systems that differ significantly in their bonding character, such as organic molecules, molecular crystals, hydrogen bonded systems, semiconductors and semimetals. In addition, extensive observations by our group show that the DP-GEN indicator, which derives from the variance of the predictions within an ensemble of DNN models, works equally well for different applications zhang2018reinforced; zhang2018deepcg. These observations are further supported by recent work by other groups who used closely related indicators in applications to a variety of different systems smith2018less; musil2019fast. We are left with the sampler, which may require case specific strategies. We are currently investigating this issue in a range of materials, finding that in all cases the search for optimal sampling strategies is facilitated by the modular structure of DP-GEN. We will present specific examples in future work.

Last but not least, one should be aware that DP-GEN scheme may fail in some circumstances. We think that this should occur most likely when the sampler and/or the indicator fail. For example, the sampler could fail when the configuration space has high dimensionality and large free energy barriers prevent exploring important configurations. In these situations, specifically designed good reaction coordinates might be necessary. Additional difficulties may be due to the indicator. To the best of the authors’ knowledge, a rigorous mathematical theory of the indicator is missing. A large value of the proposed indicator is only a sufficient, not a necessary, condition for poor performance of a DP model. There may be situations in which the physics is poorly described by a model, yet the corresponding ensemble of predictions has small variance. We did not face these difficulties in the present investigation but the reader should be aware that systematic validation tests should always be performed before using a DP model to explore new physics.

V Supplementary Materials

The smooth edition of the deep potential zhang2018end model is adopted in this work. The cut-off radius is set to 9 Å. The $1/r$ terms in the network construction is smoothly switched-off by a cosine shape function zhang2018end from 2 Å to 9 Åso that the discontinuity due to the cut-off is removed. The filter (embedding) net is of size $\{25,50,100\}$ , and the fitting net is of size $\{240,240,240\}$ . A skip connection is built between two neighboring layers, so the architecture of the network is ResNet-like he2016deep. The Adam stochastic gradient descent method Kingma2015adam is adopted to train the models, with a learning rate starting at $5.0\times 10^{-4}$ and exponentially decaying to $1.8\times 10^{-8}$ in 400,000 training steps. Four models with the same data and training setting, but different parameter initializations, are trained to estimate the model deviation in the force prediction. After all the data are collected, the final model is trained with 1,2800,000 training steps.

Exploration

Table S1, S2, and Table S3 report the exploration strategy in each iteration for pure Al, pure Mg, and Al-Mg alloy systems, respectively. During the exploration, if the model deviation of a configuration falls in the range [0.05, 0.15] eV/Å in the case of pure Al and Al-Mg alloy, or in range [0.03,0.13] eV/Å in the case of pure Mg, then the corresponding configuration is selected for labeling. The number of atoms in each crystalline structure, the total number of explored and labeled configurations of each crystal structure are reported by Tab. S4.

Labeling

The DFT simulation is carried out by the Vienna ab initio simulation package (VASP) version 5.4.4 kresse1996efficiency; kresse1996efficient, within the Perdew-Burke-Ernzerhof generalized gradient approximation. The kinetic energy cutoff for the plane wave expansion is set to 600 eV, and the K-points is set with the Monkhorst-Pack mesh monkhorst1976special at the spacing $h_{k}=0.08~\AA^{-1}$ . The order 1 Methfessel-Paxton smearing method with $\sigma=0.25$ eV is adopted. The self-consistent field (SCF) iteration will stop when the total energy and band structure energy differences between two consecutive steps are smaller than $10^{-6}$ eV.