\maketitle
Keywords: Kohn-Sham density functional theory, neural networks, reactive force fields, potential energy surfaces, machine learning
Growing interest in nanoscale phenomena has driven the need for accurate computational models to study chemical properties at an atomistic level. Advances in high performance computing now enable powerful simulation tools such as Born-Oppenheimer molecular dynamics (BO-MD) and extended Monte Carlo techniques that employ first principles quantum chemistry (QC) methods (e.g. Kohn-Sham density functional theory (KS-DFT)) on increasingly complex reactions. However, many genuine nanoscale systems of interest are currently too large to be studied with first principles QC and thus require the use of computationally efficient atomistic potentials.
Atomistic potentials approximate the potential energy surfaces (PES) for atomic systems by mapping potential energies and forces as functions of atomic positions. Physical potentials are parameterized to fit analytical expressions for the known physics of pairwise and many body interactions, e.g.: Lennard-Jones potentials, cite:jones-1924-deter-molec-field classical force fields, cite:brooks-1983-charm,rappe-1992-uff,casewit-1992-applic,cornell-1995-secon-gener,cornell-1996-secon-gener the modified and standard embedded atom methods (MEAM and EAM) cite:daw-1983-semiem-quant,baskes-1992-modif, reactive force fields (ReaxFF) cite:duin-2001-reaxf, and charge-optimized many-body (COMB) potentials cite:liang-2012-variab-charg. Another intriguing category are mathematical potentials that “learn” the PES directly through a minimization of residuals with no a priori knowledge of the underlying physics, e.g. Gaussian regression functions cite:rasmussen-2004-gauss-proces and artificial neural networks (NN) cite:haykin-2009-neural-networ. We refer to these as mathematical potentials.
Physical potentials have been used for decades to capture the underlying physics of complex systems. It is substantially more challenging to develop so-called reactive potentials that can accurately model the underlying physics of systems where bonds are broken and formed. One approach is to develop potentials that use chemical bond-order dependent functional terms that are fit exclusively (or primarily) to QC datasets cite:liang-2013-react-poten. A challenge with these approaches is knowing what data should be used for training and how to parameterize potentials that best fit the data.
In contrast, mathematical potentials are becoming more popular in chemical applications, and recent work has shown reasons for excitement.cite:behler-2011-neural,behler-2014-repres. Specifically, recent descriptive models from Behler and Parrinello cite:behler-2007-gener-neural have expanded the applications of neural networks to “high-dimensional” systems that can account for variable numbers of atoms, multiple compositions, and reactions involving thousands of atoms. Such networks have already been implemented on a large range of systems, including: Si bulk structures cite:behler-2008-press, water clusters cite:morawietz-2013-densit-funct, Cu surfaces cite:artrith-2012-high, ZnO cite:artrith-2011-high, and even a quaternary system of Cu/Au/H/O cite:artrith-2015-grand-cu. This opens the door for mathematical potentials to be developed that are accurate and transferable across diverse bulk, surface, and cluster regimes.
In this work, we compare the performance of a widely used physical potential, ReaxFF, with the recently developed Behler-Parrinello neural network (BPNN) potential. We have trained both to subsets of a full dataset comprised of
We have benchmarked both potentials to our QC dataset that contains information from KS-DFT bulk equation of state (EOS) data, vacancy formation energies, surface energies, adatom diffusion profiles, slipping barriers, and cluster binding energies. Parameterization of ReaxFF potentials were automated using the Monte Carlo Force Field optimization (MCFFopt) tool in ADF cite:velde-2001-chemis-adf,iype-2013-param-monte. Our BPNN was parameterized using the Atomistic Machine-learning Potentials (AMP) code from the Peterson group at Brown University cite:peterson-amp. This allows feed-forward neural networks to be developed inside the atomic simulation environment (ASE) cite:atomic-simul-envir. All details of the trained BPNNs are stored in a JSON file which can be found in the supporting information (SI) file.
KS-DFT calculations for training set data were performed using the Vienna ab initio simulation package (VASP) cite:kresse-1993-ab,kresse-1994-ab,kresse-1996-effic,kresse-1996-effic2 with the Perdew-Burke-Ernzerhof generalized gradient approximation (GGA-PBE) cite:perdew-1996-gener-gradien,perdew-1997-gener-gradien exchange-correlation functional. Core electrons were described using the projector augmented wave function (PAW) cite:blochl-1994-projec-augmen,kresse-1999-from-ultras. k-points were represented using Monkhorst-Pack grids cite:monkhorst-1976-special-point with a density of at least 14 × 14 × 14 for a single atom of Au in the primitive ground state configuration. Kohn-Sham orbitals were expanded up to energy cutoffs of at least 300 eV to attain an energy convergence of at least 5 meV/atom. All calculations involving relaxations were completed with relaxation criteria of
The full KS-DFT training set contained 9,972 calculations that included 905 bulk, 1,022 surface, and 8,045 cluster configurations. The majority of these calculations (9,076 calculations) were taken from coordinate relaxation steps performed by VASP. These structures are the incremental steps taken from its initially guessed positions to the ground state configurations predicted by GGA-PBE. Each of the structures in a particular relaxation are very similar from one relaxation step to the next. The remaining 896 calculations are either the local ground state configurations or images from optimized NEB calculations. Our bulk Au data were obtained by plotting EOS data for a variety of bulk structures. Vacancy formation and diffusion calculations were also included in the bulk dataset. Our surface dataset includes calculations on fcc(111) surfaces as well as a variety of fcc(100) surface diffusion pathways that were originally generated in previous work by Pötting et. al. cite:potting-2010-self-diffus. The training set used single-point energies on the latter coordinates (without geometry relaxations) calculated using the methods listed above. Our cluster data include various 3D ordered, planar, and disordered structures that contain up to 126 atoms. The SI file has further details about the data.
Bond order based reactive force fields, such as Tersoff cite:tersoff-1988-new, Brenner cite:brenner-1990-empir, and ReaxFF cite:nielson-2005-devel-reaxf,duin-2001-reaxf potentials, differ from classical force fields, such as UFF cite:casewit-1992-applic,rappe-1992-uff, CHARMM cite:brooks-1983-charm, or AMBER cite:cornell-1995-secon-gener,cornell-1996-secon-gener, which require that defined bonds remain fixed over the course of a simulation. ReaxFF potentials developed for Au and other metals normally employ three separate energy terms as seen in Equation ref:eqn-base-reax. cite:jarvi-2008-devel-reaxf,keith-2010-react,cabrera-trujillo-2015-theor
\begin{eqnarray} Etotal = Ebond + Eover + Evdw \label{eqn-base-reax} \end{eqnarray}
$Ebond$ is for bond energies of atom pairs, $Eover$ is an energy penalty to prevent overcoordination, and $Evdw$ accounts for van der Waals interactions and interatomic repulsions when interatomic distances are too small. ReaxFF potentials can also be parameterized to include 3-body terms which provide energy contributions from valence angles between sets of three Au atoms. Backman et. al. developed a Tersoff potential for Au that involves 3-body terms cite:backman-2012-bond, but these terms are not always added to ReaxFF potentials for metals due to increased computational cost. Our 3-body terms have the same form as valence angle interactions in hydrocarbon ReaxFF potentials cite:nielson-2005-devel-reaxf. We report comparisons in equations of state and timings for ReaxFF using 3-body and 2-body terms in the SI file. Future work will discuss these topics in greater detail.
We parameterized our Au ReaxFF using the MCFFopt tool implemented in ADF cite:velde-2001-chemis-adf,iype-2013-param-monte. MCFFopt seeks to minimize an objective function by randomly changing force field parameters within a predefined range. The Monte Carlo nature of this process allows some parameter changes that increase the objective function. This “annealing” allows the optimizer to sample a larger parameter space and potentially produce multiple distinct parameter sets. This approach can also find parameter sets with less total error than the traditional parabolic search parameter optimization cite:iype-2013-param-monte. Further information on running the MCFFopt procedure and optimized force field parameters are available in the SI file.
Au ReaxFF potentials appear to have an optimal training set size. Fitting to larger training sets does not always improve the quality of the ReaxFF potential, and this overfitting is found to bias predictions toward certain geometry types. As a result, the ReaxFF training set was constructed using the 848 ground state geometries from within the training set. Out of these geometries, the number of calculations classified as bulk, surface, and cluster structures are roughly equal. During ReaxFF parameterizations, each geometry in the training set is also assigned a weight depending on its relative importance in the overall fitting procedure. Our goal was to produce a ReaxFF potential with reasonable accuracy across these three different structure regimes, so most of the geometries were given a weight of one (specific details are given in the SI file). In principle, one could increase weights to parts of the PES so that properties, such as desired lattice constants, bulk moduli, or barrier heights would be accurately reproduced. However, weighting a potential in this way will affect its ability to make accurate predictions in less-weighted regions of the PES.
Figure ref:fig-reax-train shows the error distribution of residual error between the trained ReaxFF and KS-DFT training set data labeled by geometry type. Errors in bulk data greater than 0.2 eV stem from an unphysical convex region in the ReaxFF functional form which causes bulk EOS data to significantly deviate from the KS-DFT data at atomic volumes ranging from 60-200 Å3/atom. Since these atomic volumes fall outside those found in most simulations involving bulk and surface structures of Au, these inaccuracies are not a cause for significant concern. However, large errors in bond energies for pairs of atoms at intermediate distances may be problematic for molecular clusters. Images of the entire EOS for each bulk structure can be found in the SI.
A predefined validation set consisting of 238 calculations (out of the total 9,972 KS-DFT calculations) was set aside to test the transferability of predictions from our ReaxFF and BPNN potentials. This validation set was chosen to represent a variety of different Au structure types which are represented in the results section of this work. By reporting probability distributions for both the training and validation sets, we can determine the degree that our potentials show selection bias. For an ideal fitting procedure, the probability distributions for both the training and validation set would match, and any differences between the two would signify an over- or under-sampling. Figure ref:fig-reax-valid shows the residual error for the validation calculations labeled by geometry type. Significant deviations were found in bulk and cluster calculations from the validation and training set data.
The NN is a machine learning algorithm. Unlike ReaxFF, potentials constructed from the NN have no physical basis, making them highly flexible, but also unsuitable for extrapolation. Nevertheless, NN potentials are growing in popularity due to their abilities to accurately characterize a PES from QC calculations. BPNN potentials were selected for this study because they can be trained to relatively large training sets cite:behler-2007-gener-neural. BPNN potential fitting is facilitated through symmetry functions and utilization of multiple feed-forward NNs, one for each chemical species in the system. These modifications eliminate many of the shortcoming of traditional Cartesian NN potentials that are only applicable to systems with a fixed number of atoms.
Unlike other mathematical potentials, BPNN potentials cite:behler-2011-atom,behler-2007-gener-neural utilize a cutoff radius,
The general structure of a feed-forward NN consists of nodes in an input layer, one or more hidden layers, and an output layer. In our BPNN potential, the nodes of the input layer are the Cartesian-coordinates of each atom in the unit cell. Each hidden layer is a linear combination of the values of the nodes from the previous layer. Each layer is also multiplied by an activation function (which often has a bounded non-linear form) to allow the NN potential to fit to arbitrary functions. The most accurate BPNN potential that we produced, and report in this work, utilizes four hidden layers with 40 nodes per layer and a hyperbolic tangent activation function. These specifications make our BPNN potential large compared to other BPNN potentials and more at risk for overfitting. However, the root-mean square error (RMSE) of the validation set is similar to that of the training set, showing that overfitting has not occurred cite:behler-2015-const. For further details on the theory behind BPNNs, we refer to previous work cite:behler-2007-gener-neural,behler-2011-atom.
We trained BPNN potentials using AMP, a code produced by the Peterson group at Brown University cite:peterson-amp. This software conveniently interfaces with the Atomic Simulation Environment (ASE) software package cite:atomic-simul-envir for ease of reusability and reproducibility. The trained calculator parameters used by AMP are included in the SI file.
Of the 9,972 total calculations, 9,734 were used for training the BPNN potential. Figure ref:fig-neural-train shows the error distribution from the training set. The mean,
Figure ref:fig-neural-valid shows the error distribution for the validation dataset. Overfitting can be identified by a divergence between the RMSE of the training set and validation set data. In this case, the distribution is clearly not normal and arises from some underrepresented data in the training set, notably the fcc(100) terrace and dimer diffusion pathways (discussed below).
We now benchmark the performance of the BPNN and ReaxFF potentials against KS-DFT energies across three different material regimes: bulk, surface, and molecular cluster structures. Both of our generated potentials can provide reasonably accurate descriptions of Au in the different material regimes. In general we find that ReaxFF potentials are more readily overfit, less transferrable to applications involving clusters of 126 atoms or fewer, and overall less accurate than the BPNN. However, ReaxFF potentials demonstrate a notable strength by predicting barrier heights that resemble those found in their training sets. BPNN potentials in general are significantly more accurate than ReaxFF potentials, but they require significantly larger training sets to ensure well-balanced fitting. As explained below, they also currently bring substantially higher computational cost than ReaxFF potentials.
EOS data for face centered cubic, simple cubic, and diamond structures are shown in Figure ref:fig-bulk-eos. All training and validation calculations are fit to a 3rd order inverse polynomial cite:alchagirov-2003-reply-commen. The metrics for each fit are included in Table ref:tbl-eos. Results for the body centered cubic and hexagonal close packed EOS data are similar to the face centered cubic curve. Fits to all curves can be found in the SI file.
Figure ref:fig-bulk-eos shows that the EOSs are very well represented by our BPNN potential. Validation set data are also well behaved, indicating that overfitting has not occurred. Metric data shown in Table ref:tbl-eos shows excellent agreement in the minimum volume, minimum energy, and bulk modulus found using KS-DFT results. Data for the hcp and bcc structures shown in the SI file are reproduced similarly well.
We find that ReaxFF potentials with 3-body terms have substantially better fits compared to force fields which do not include 3-body interactions (see cite:keith-2010-react). However, in both cases ReaxFF exhibits an unphysical convexity of the bond energy curve that creates problems manifested by large residual errors that can reach as high as ± 1 eV/atom at volumes away from the minimum energy volume. Many simulations sample regions in the vicinity surrounding the minimum volume, so these deviations are not shown in Figure ref:fig-bulk-eos. Data from Table ref:tbl-eos shows reasonably good agreement for the equilibrium volume and minimum energy of the three structures. Bulk moduli are underpredicted by
| Structure | Minimum volume (Å3) | Minimum energy (eV) | Bulk Mod. (GPa) |
|---|---|---|---|
| KS-DFT-fcc | 17.97 | -3.23 | 147 |
| BPNN-fcc | 17.99 | -3.23 | 145 |
| ReaxFF-fcc | 17.60 | -3.22 | 122 |
| KS-DFT-sc | 20.73 | -3.02 | 110 |
| BPNN-sc | 20.66 | -3.02 | 110 |
| ReaxFF-sc | 21.29 | -2.96 | 84 |
| KS-DFT-diam | 29.04 | -2.51 | 56 |
| BPNN-diam | 28.98 | -2.51 | 57 |
| ReaxFF-diam | 31.92 | -2.54 | 37 |
Vacancy formation energies (
\begin{eqnarray} E_v = E_f - \frac{n_0 - 1}{n_0} E_i \label{eqn-vac} \end{eqnarray}
Our BPNN vacancy formation predictions are systematically overestimated by
The residual errors for structures with concentrations below 0.04 vacancies/atom are very low (less than 0.006 eV/atom, even for the point in the validation set having
Figure ref:fig-vacancy-diffusion shows the calculated bulk vacancy diffusion barrier using a vacancy concentration of
The training set for the ReaxFF potential in Reference citenum:keith-2010-react contains 166 surface diffusion barrier calculations from GGA-PBE using the SEQQUEST code cite:schultz-2002-seqques. NEB calculations with VASP were not used to recalculate the minimum energy pathways, but we recalculated single point energies on these structures using GGA-PBE in VASP to be consistent with the rest of our training set. Since NEB calculations were not done, there are significantly fewer points sampling the PES for these pathways compared to other pathways (8-10× fewer in most cases). Consequently, our BPNN fits to these pathways are expected to be less accurate compared to other pathways obtained from NEB calculations.
Figure ref:fig-full-diffusion contains recreations from Figure 2 (a \& b) in Ref. citenum:keith-2010-react using the BPNN potential and ReaxFF potential. Note that the terrace and dimer diffusion pathways are not included in the training set for either potential, and they represent predictions by both potentials. For the terrace diffusion pathway, the ReaxFF potential performs quite well and shows that the ReaxFF potential can provide very accurate predictions of barrier heights when the training set contains similar pathways. The BPNN potential, which contains more than 10× the training set data as the ReaxFF potential, can reasonably produce this adatom diffusion barrier but residuals fall between 0.2-0.3 eV. On the other hand, for a different adatom diffusion barrier, the BPNN potential predicts the dimer diffusion pathway quite well while the ReaxFF potential has higher residual errors between 0.1-0.2 eV. Larger training sets can be expected to reduce errors in both potentials, but reparameterization of these potentials with a larger training set will undoubtedly impact the accuracy when predicting other pathways.
To assess the performance of these potentials under a wide range of adatom diffusions, Figure ref:fig-barrier-residuals shows the residuals for all 144 fcc(100) surface diffusion calculations. Solid shapes represent training set data and hollow shapes represent validation set data. Residuals are the same as those shown in Figure ref:fig-full-diffusion. Our ReaxFF potential (which has roughly 1/3 of its training set devoted to surface calculations) has 86.1% of these structures falling within a ± 0.1 eV tolerance of error. For the BPNN potential (with roughly 1/10 of its training set devoted to surface calculations), has 52.1% of these structures fall within a ± 0.1 eV tolerance of error.
Many of the calculations from the BPNN potential are underestimated compared to the reference KS-DFT data, signifying (as stated above) that these structures come from a poorly sampled region of the PES and improvements could be attained with more training. For the ReaxFF potential, errors appear to be less systematic, showing improved accuracy would require more training to specific pathways. In practice, both ReaxFF and BPNN potentials are normally trained with a specific application in mind, and so training sets, particularly those for ReaxFF potentials, can be smaller.
A slipping barrier is the minimum energy pathway required for a certain number of mono-layers of atoms to move from their ground state site to the next most adjacent site of the same kind. Slipping barriers were performed on fcc(100) and fcc(111) surfaces for one and two layers in a five layer slab. Figure ref:fig-111-slipping shows the single-layer slipping barrier for the fcc(111) surface. Both models find almost identical energies as KS-DFT (within 0.05 eV). We can see that the ReaxFF potential finds a metastable intermediate instead of a single barrier as found by KS-DFT and the BPNN potential. This ReaxFF potential also finds metastable intermediates when slipping in a different direction primarily over bridge sites (see the SI file), but residual errors are even lower. The very small difference in energies makes it difficult to assess if these are due to fitting errors or an unphysical component within the ReaxFF potential. Either way, both potentials can reproduce low energy slipping barriers within 0.05 eV with sufficient training.
Calculations on clusters up to 126 atoms make up
Figure ref:fig-6atom-md depicts the path taken by the BPNN BO-MD simulation (red) over the course of 2,000 time steps. Once every 100 steps we validated the energy using KS-DFT. The residuals are less than 0.05 eV/atoms for the BPNN potential, including the structure of the global minimum. We re-ran this simulation several times throughout development of the database. The first attempt at performing the described BO-MD simulation was with a dataset of
The 2,000 structures generated from the BO-MD run with our BPNN potential were then calculated using the ReaxFF potential. In this case, the ReaxFF potential did not identify the same minimum energy configuration of the six atom system. However, the cohesive energies of structures resembling the planar cluster are fairly consistent with KS-DFT data. Although the presented data shows situations where ReaxFF is not accurate, we note that this may signify an area where ReaxFF could be extended with additional functionality. For example, metal-metal bonds in small clusters could be treated with functional forms different than those used for bulk metal-metal bonds. This would likely correct systematic deviations, but such re-parameterizations may also adversely affect other structure types and/or increase computational cost. We note that Narayanan et. al. have reported a hybrid bond order potential that uses a screened Lennard-Jones term for bulk structures in combination with a highly trained Tersoff potential for smaller regimes cite:badri-2015-descr. This is a possible work-around to make other physical potentials accurate across different size regimes.
A similar exploration for multiple local minima was implemented on a 38 atom cluster using minima hopping techniques cite:goedecker-2004-minim. This exploration of minimum energy structures works through a series of fixed temperature NVT BO-MD simulations followed by geometric optimization requiring a significant number of calculations between each iteration. After each iteration, the minimum geometry is stored and perturbed before restarting its search. The resulting minima predicted from 125 such iterations are shown in Figure ref:fig-38atom-minima.
Again, this approach located a lower energy minimum than the starting point geometry. The largest energy difference between minima occurred during the first iteration of the process. After this initial step, the energies do not change as dramatically. This can be interpreted as a shift into a local minimum energy basin (a group of configurations with similar atomic positions and energies) which the BPNN potential proceeds to explore in the next 124 minima. A more complete analysis of the 38 atom Au cluster space would be time consuming and is beyond the scope of this work. Despite demonstrating low residual errors, the BPNN potential does not correctly predict the lowest energy structure determined by KS-DFT in this set of minima. Regardless, it is still capable of distinguishing between configurations in different basins, and thus could be a valuable tool for exploring minimum energy structures in conjunction with KS-DFT calculations.
Residual errors for the ReaxFF potential are consistently lower by -0.11 eV/atom compared to KS-DFT. If energetics are shifted by this amount (as show in the top of Figure ref:fig-38atom-minima) one finds that the trend in relative energies is in reasonable agreement with KS-DFT, although our ReaxFF potential does not correctly predict the lowest energy configuration either. The performance of the ReaxFF potential for clusters could always be improved by adding more cluster data to its training set, but we found that doing so rapidly deteriorates its ability to calculate bulk and surface properties. As a result, we do not recommend using ReaxFF in its standard formalism for applications involving clusters with fewer than 126 atoms.
An important aspect of these modeling approaches is their computational cost. This includes the time needed to produce the necessary QC training sets, train the potentials, and the time needed to run the calculations. Implementation and training of parameters for both the ReaxFF and BPNN potentials can be automated using instructions in our SI file, thus reducing the time needed to learn how to train potentials. The generation of meaningful QC data is also a significant bottleneck in time, particularly for BPNN potentials that require large training sets to be accurate. This is simplified in part by generating NEB data and geometry optimizations which contain many valuable calculations on which physical and mathematical potentials can be trained. One of the best ways to speed the progress of developing accurate and transferrable potentials is to make data and methods freely available and easily accessible.
A fair comparison between calculation times between ReaxFF and BPNN potentials is not currently possible. The BPNN potential we developed used a Python code that is still in early stages of development cite:peterson-amp. In comparison, ReaxFFs and other force field codes have been implemented in the LAMMPS program which is already a high performance computing code. cite:plimpton-1995-fast-paral. Using available open source tools, BO-MD simulations on the 6 atom cluster using the C-compiled ReaxFF code performs
We have trained ReaxFF and BPNN potentials using subsets of
In applications on bulk structures, our BPNN almost exactly reproduces reference QC data of equations of state, while the ReaxFF potential is less accurate, particularly at atomic volumes that extend far beyond the equilibrium structures. When modeling surface structures and adatom diffusions, both the ReaxFF and BPNN potentials perform quite well with sufficient training, but obtaining a BPNN potential having comparable or higher accuracy than ReaxFF for adatom diffusions requires substantially larger training sets. For clusters, the BPNN potential exhibits essentially negligible residual errors compared to the KS-DFT calculations it was trained to, while the ReaxFF potential exhibits sizable systematic errors of 0.11 eV. This highlights the challenge of developing a physical potential that is accurate across bulk, surface, and cluster data. Increasing the size of the training set for the ReaxFF potential to include more cluster data was found to be detrimental to the accuracy of bulk and surface data, thus showing an area needing improvement in terms of ReaxFF functionality.
Although BPNNs can be trained to the desired level of accuracy, the computational cost, both upfront in the form of training set data and during calculation time, are currently substantially higher than ReaxFF potentials. Nevertheless, BPNN potentials are very promising if trained for specific applications (hence requiring smaller training sets) and they will be highly intriguing as computational developments enable faster runtimes. Since accurate BPNN potentials contain substantially larger numbers of parameters than most physical potentials, it is unlikely that BPNN potentials will ever be as fast as ReaxFF potentials, but we have demonstrated that BPNN potentials can be trained to be substantially more accurate.
JRK and JRB gratefully acknowledge support from the National Science Foundation under grant number CBET-1506770. JAK and MCG gratefully acknowledge financial support from the R.K. Mellon Foundation and the University of Pittsburgh’s Department of Chemical & Petroleum Engineering. They also thank the University of Pittsburgh Center for Simulation and Modeling for computational support.
bibliography:./manuscript.bib