Database optimization for empirical interatomic potential models

# Database optimization for empirical interatomic potential models

Pinchao Zhang    Dallas R. Trinkle Department of Materials Science and Engineering, University of Illinois, Urbana-Champaign
July 9, 2019
###### Abstract

Weighted least squares fitting to a database of quantum mechanical calculations can determine the optimal parameters of empirical potential models. While algorithms exist to provide optimal potential parameters for a given fitting database of structures and their structure property functions, and to estimate prediction errors using Bayesian sampling, defining an optimal fitting database based on potential predictions remains elusive. A testing set of structures and their structure property functions provides an empirical measure of potential transferability. Here, we propose an objective function for fitting databases based on testing set errors. The objective function allows the optimization of the weights in a fitting database, the assessment of the inclusion or removal of structures in the fitting database, or the comparison of two different fitting databases. To showcase this technique, we consider an example Lennard-Jones potential for Ti, where modeling multiple complicated crystal structures is difficult for a radial pair potential. The algorithm finds different optimal fitting databases, depending on the objective function of potential prediction error for a testing set.

34.20.Cf

## I Introduction

Atomic-scale simulations have the capability to predict the properties of defect structures that are often inaccessible by experimental techniques.Daw and Baskes (1984); Baskes (1992); Xu and Moriarty (1996); Daw et al. (1993); Wirth et al. (2000); Lilleodden et al. (2003); Mishin et al. (2001); Woodward et al. (2008) These predictions require accurate and efficient calculations of energies and forces on atoms in arrangements that sample a variety of atomic environments, and may represent even different binding configurations. Accurate quantum mechanical methods are difficult to scale to large systems and long simulation times, while empirical interatomic potentials offer increased computational efficiency at a lower level of accuracy. Maximizing the efficiency of computational material science studies requires the development of potentials that are transferrable, i.e., capable of predicting properties outside their fitting range, and accurate for static and dynamic calculations.

However, without direct transferable derivations of interatomic potentials from quantum mechanical methods, empirical interatomic potentials require high-dimensional non-linear fitting. Many functional forms for empirical potentials have been proposed, including embedded-atom method (EAM)Daw et al. (1993); Li et al. (2003), modified embedded-atom method (MEAM)Baskes et al. (1989); Baskes and Johnson (1994); Baskes (1992) and charged-optimized many-body potential (COMB).Yu et al. (2007); Shan et al. (2010) There have been multiple implementations of different potential functional forms for various materials.Mishin et al. (2001); Lenosky et al. (2000); Hennig et al. (2008); Park et al. (2012); Liu et al. (1996); Li et al. (2003); Fellinger et al. (2010); Sheng et al. (2011) Even for the same type of materials, such as CuFoiles et al. (1986); Zhou et al. (2001); Mishin et al. (2001) and SiBaskes (1987); Baskes et al. (1989); Tersoff (1986); Balamane et al. (1992); Lenosky et al. (2000), different empirical interatomic potential models are proposed for different applications with different transferabilities. There are advanced techniques to optimize the potential parameters based on a weighted least-squares regression to a fitting database of experimental or quantum mechanical calculation data,Foiles et al. (1986); Daw et al. (1993) including the force-matching methodErcolessi and Adams (1994) for empirical interatomic potential parameter optimization. In force-matching, a fitting database includes quantum mechanical force calculations for diverse atomic environments to obtain realistic empirical potential models. To study the transferability of the empirical potential model, Frederiksen et al. applied Bayesian statistics to empirical interatomic potential models: instead of using the best fit, an ensemble of neighboring parameter sets reveal the flexibility of the model.Frederiksen et al. (2004) They showed that the standard deviation of the potential prediction of structure property function is a good estimate of the true error. However, even with these advances, the determination of empirical interatomic potentials relies on the selection and weighting of a fitting database without a clear, quantitative guide for the impact on predictions.

To address the issue of fitting database selection, we present an automated, quantitative fitting-database optimization algorithm based on prediction errors for a testing set using Bayesian statistics. We construct an objective function of the prediction errors in the testing set to optimize the relative weights of a fitting database. This includes the addition or removal of structures to a fitting database when weights change sign. We demonstrate the viability of the optimization algorithm with a simple interatomic potential model: Lennard-Jones potential fitting of Ti crystal structures. We choose this example as a radial potential has difficulty describing the stability of different crystal structures of a transition metal. The new algorithm also helps to understand the transferability of the empirical potential model for the structures in the testing set.

We start with a brief review of empirical potential models and parameter optimization using a fitting database. Next, we discuss Bayesian error estimation as it applies to our problem. Then we define an objective function with a testing set, and use this quantitative measure to devise an algorithm to optimize a fitting database. Lastly, we demonstrate this new approach on an example system with clear limitations: Lennard-Jones potential for titanium.

## Ii Interatomic potential models

The total energy and forces for the structure of interest are the most basic quantities to calculate since they determine the structural properties. In particular, we are interested in predictions that are derived from energies of atomic arrangements. In atomic-scale simulations, a structure is a set of atomic positions with chemical identities : . The total energy of a structure is with forces . Density-functional theory (DFT)Hohenberg and Kohn (1964); Kohn and Sham (1965) calculations can provide accurate structural energies and forces, but their computational cost limits them to simulation system with at most a few thousand atoms. Other structural properties are derived from energies and forces, and so without loss of generality, we develop our approach based on energies and forces.

Parameterized empirical interatomic potentials offer a computationally efficient alternative to DFT. Potentials provide approximate energies and forces for atomic configurations that are inaccessibly large for DFT calculations. Generally, an empirical interatomic potential functional can be written as

 Eα(θ)≡E({(→Ri,χi)};θ)=12!∑mnVχmχn2(→Rm−→Rn;θ)+13!∑mnlVχmχnχl3(→Rm−→Rn,→Rn−→Rl;θ)+⋯, (1)

where are parameters, and is an interatomic potential function between atoms of chemical identity . Symmetries of the potential functional form, such as permutation symmetry, translational symmetry, rotational symmetry, etc., can simplify the functional form. A general empirical interatomic potential that reproduces all DFT energy calculations accurately is computational intractable, since it would require a large number of many-body terms. Rather, we are interested in simpler potentials that provide accurate results for a smaller range of atomic configurations including perfect crystals and defect structures under various thermodynamic conditions; this includes potentials that may not be easily written in the form of Eqn. (1), such as EAM and MEAM potentials.

The optimal potential parameters derive from comparison to predictions of a database of DFT calculations, and the performance of the potential is evaluated by a testing set of structure property predictions. A fitting database is a set of structure property functions with an associated structure and positive (relative) weight . While a single structure will often have multiple structure property functions with unique weights for fitting, we simplify our notation by indexing on the structure; in what follows, sums over structures indicate sums over all members of the database . The structure property function may be a scalar such as the energy (relative to a reference structure), vectors such as forces on the atom of the structures, stress tensors and more complicated structure property functions such as lattice constant, bulk modulus or vacancy formation energy or anything that can be defined from the energy . In the fitting database, we will compare the structure property functions evaluated using an empirical potential, with the values from DFT, though other choices are possible, such as experimental data. In the weighted least-squares (described later), we impose the trivial constraint , as only relative values of are important. A testing set is a set of structures with structure property functions . In a testing set, we will compare structure property functions evaluated using an empirical potential with either values from DFT, or using Bayesian sampling of the empirical potential, following Frederiksen et al.,Frederiksen et al. (2004) which we will discuss in Section III. There are no relative weights for structures in a testing set; rather, these represent a set of predictions whose errors we will evaluate.

In order to assess the prediction errors of the structure property functions, we define the prediction error function as

 ϵ2α(θ)=∥Aα(θ)−Aα∥2, (2)

where is the structure property function from the empirical atomic potential with parameters , is the structure property function from DFT, and denotes the 2-norm of a -dimensional vector

 ∥x∥2=d∑a=1|xm|2. (3)

We will take the error evaluation of the energy differences between two structures and forces as examples. For energy calculations, the structure property function is

 Aα=Eα−E0, (4)

where is the energy of a reference structure, . The potential energy prediction error is

 ϵα(θ)=|(Eα(θ)−E0(θ))−(Eα−E0)|. (5)

The force predictions errors are

 ϵ2α(θ)=∥fα(θ)−fα∥2. (6)

Then, the weighted summed squared error function for a fitting database is

 S(θ,F)=∑α∈Fwαϵ2α(θ), (7)

## Iii Bayesian Error Estimation

We introduce Bayesian sampling to estimate the errors of structure property function predictions and quantitatively analyse the relative weight values in the fitting database. Given a fitting database, we calculate the prediction of structure property function and the error of the structure property function. We then derive the analytical expression of the gradient of the Bayesian errors with respect to the weights, . These gradients provide quantitative information on how structure property functions in the fitting database influence the Bayesian predictions of structure property functions in the testing set though weight change.

Bayesian statistics treats model parameters as random variables with a probability distribution given by a posterior distribution.Bolstad (2007) According to the Bayes’ theorem, the posterior distribution of the parameters is a product of the prior distribution and the likelihood function ,

 P(θ;F)∝π(θ)×L(θ;F), (8)

where the prior distribution includes the information about the potential model before the we take the fitting data into account. Here we use the maximally unbiased prior distribution of a uniform distribution over a measurable set of allowed parameters sets,

 π(θ)=[∫Hdθ]−1, (9)

though other choices are possible. Assuming the errors are independent and identically normally distributed, the likelihood function isHogg et al. (2004); Frederiksen et al. (2004)

 L(θ;F)∝exp(−1W∑α∈Fwαϵ2α(θ)), (10)

where

 W=infθS(θ,F). (11)

The log-likelihood function is proportional to the squared error function, .

 logL(θ;F)=−S(θ,F)W=−1W∑α∈Fwαϵ2α(θ), (12)

Since the logarithm is a monotonically increasing function, minimizing is equivalent to maximizing the log-likelihood function. The maximum likelihood estimate (MLE) of the parameters is a function of the fitting database , and .

The Bayesian prediction of a function is the mean

 ⟨A(θ)⟩F=∫P(θ;F)A(θ)dθ∫P(θ;F)dθ=∫HL(θ;F)A(θ)dθ∫HL(θ;F)dθ. (13)

All the averages are implicit functions of the relative weights in the fitting database. The Bayesian error is the mean squared error of the Bayesian prediction:

 ⟨ϵ2A(θ)⟩F=|⟨A(θ)⟩F−A|2+varF(A(θ)), (14)

where is the variance of the Bayesian prediction. The covariance of two functions and represents the correlation between two functions:

 covF(Aα(θ),Aβ(θ))=⟨Aα(θ)Aβ(θ)⟩F−⟨Aα(θ)⟩F⟨Aβ(θ)⟩F. (15)

The derivative of a Bayesian prediction with respect to weight is

 ∂⟨A(θ)⟩F∂wα=∂∂wα∫HL(θ;F)A(θ)dθ∫HL(θ;F)dθ=covF(A(θ),∂logL(θ;F)∂wα). (16)

Note that

 ∂logL(θ;F)∂wα=−ϵ2α(θ)+∂W∂wαlogL(θ;F)W. (17)

The derivative of with respect to weight is found using the chain rule,

 ∂W∂wα=∂S(θ;F)∂wα∣∣∣θ% MLE+∑n∂S(θ;F)∂θn∣∣∣θMLE∂θMLEn∂wα=ϵ2α(θMLE), (18)

as is an extremum of . Applying Eqn. (16)–Eqn. (18) to the Bayesian error yields

 ∂⟨ϵ2β(θ)⟩F∂wα=covF(ϵ2β(θ),−ϵ2α(θ)+ϵ2α(θMLE)logL(θ)W)=−1W[CFαβ−ϵ2α(θMLE)WcovF(ϵ2β(θ),S(θ;F)]=−1W[CFαβ−ϵ2α(θMLE)W∑γwγCFβγ], (19)

where

 CFαβ=CFβα=covF(ϵ2α(θ),ϵ2β(θ)). (20)

We define the error of a structure property function in a testing set without DFT calculations by approximating the unknown DFT calculations of the structure property function with its Bayesian prediction. Based on Eqn. (14), is approximated by , and so a testing set can include structures in the absense of DFT calculations.

We need to evaluate the integral in Eqn. (13) to calculate the Bayesian predictions and Bayesian error estimation. For complicated high-dimensional, non-linear models such as empirical potentials, the integral cannot be evaluated in closed form, and the high-dimensionality makes direct numerical quadrature converge slowly. We instead use Markov Chain Monte Carlo (MCMC) to numerically integrate. The chain of will contain a set of independent samples (where is the autocorrelation length), and the numerical estimate of the mean is

 ⟨Aα(θ)⟩F=∫P(θ;F)Aα(θ)dθ∫P(θ;F)dθ≈1NN∑n=1Aα(θn), (21)

with a sampling error of . Hence, once fitting is complete, the “best” set of parameters defines the empirical potential for predictions, while the ensemble of parameters allows the estimation of errors on those predictions.

## Iv Database Optimization

We define an optimal fitting database based on Bayesian errors in the testing set. An empirical potential model should reproduce DFT calculations for a set of atomic environments described by structures in a testing set, and so the Bayesian errors of structure property functions in the testing set are the quantities of interest. Because different types of structure property functions have different units, different error magnitudes, and different degrees of freedom, we need an unbiased choice of objective function to evaluate different fitting database performances based on the Bayesian errors for the same testing set. Here, we consider the difference of the logarithm of the Bayesian errors for one structure property function for two different fitting databases, and

 ln⟨ϵ2β⟩F1−ln⟨ϵ2β⟩F2=ln⟨ϵ2β⟩F1⟨ϵ2β⟩F2=ln(1+⟨ϵ2β⟩F1−⟨ϵ2β⟩F2⟨ϵ2β⟩F2). (22)

If is small, then , and

 ln⟨ϵ2β⟩F1−ln⟨ϵ2β⟩F2≈⟨ϵ2β⟩F1−⟨ϵ2β⟩F2⟨ϵ2β⟩F2≈2⟨ϵ2β⟩F1−⟨ϵ2β⟩F2⟨ϵ2β⟩F1+⟨ϵ2β⟩F2. (23)

Then the right side of the equation is a relative difference in errors. We propose the objective function of a fitting database with testing set ,

 O(F;T)=∑β∈Tln⟨ϵ2β(θ)⟩F, (24)

so that is approximately the sum of relative differences in error. Then, the optimal fitting database minimizes the sum log Bayesian errors for a testing set . The objective function is implicitly dependent on the relative weights in the fitting database through the Bayesian error. The gradient of the objective function with respect to weight is analytically calculable (c.f., Eqn. (19)). We obtain the optimal weights in the fitting database by minimizing the objective function. Hence we will be able to compare potentials fitted with different fitting databases with respect to the same testing set.

However, the minimum of the objective function can be trivial for pathological fitting databases and testing set combinations. A pathological fitting database and testing set combination is an underdetermined fitting database, where the MLE predictions can match the true values of a DFT calculation in both the fitting database and testing set. Thus, if for any structure , then approaches logarithmically. In order to eliminate the trivial minimum of pathological databases, we introduce a threshold function ,

 t(x)={x:x≥2ϵ20x2/4ϵ20+ϵ20:x<2ϵ20, (25)

that creates a finite minimum of at . We can choose different error tolerances for each testing set structure property function. The objective function is then

 O(F;T)=∑β∈Tlnt(⟨ϵ2β(θ)⟩F), (26)

and the derivative of the objective function is

 ∂O(F;T)∂wα=∑β∈Tt′(⟨ϵ2β(θ)⟩F)⟨ϵ2β(θ)⟩F∂⟨ϵ2β(θ)⟩F∂wα, (27)

where the derivative calculations are from Eqn. (19) and Eqn. (28). The derivative of the threshold function is

 t′(x)={1:x≥2ϵ20x/2ϵ20:x<2ϵ20. (28)

Finally, note that as our likelihood function is independent of , so

 ∑αwα∂O(F;T)∂wα=0. (29)

The optimal weights are found by minimizing , and this includes the addition and removal of structures from the fitting database. According to the definition of the likelihood function, Eqn. (12), the fitting database could include any structures with DFT calculations with a non-negative relative weight value. Structures with positive weight values are structures to fit, and all the other structures that do not contribute to the fitting will have a weight of zero. The optimal weight value can be determined even for structures not presently in the fitting database. A structure is added to the fitting database if its optimal weight value is positive, as inclusion of that structure decreases the relative error in the testing set. A structure is removed from the fitting database if its optimal weight value is zero or negative, since removing the structure decreases the relative error in the testing set.

Fig. 1 outlines the new interatomic potential development algorithm based on Bayesian statistics. It starts with a conventional interatomic potential fitting procedure by selecting a potential functional form and DFT structural energies and forces forming a DFT data set. We build the fitting database with a set of structures from the DFT data set and assign each structure with a relative weight. The testing set also contains a set of structures from the DFT data set that test the ability of the potential to model atomic environment outside of the fitting data. We apply non-linear weighted least squares regression to obtain the MLE of parameters of the empirical potential model, and use Markov Chain Monte Carlo (MCMC)Givens and Hoeting (2005) sampling of the posterior distribution to generate the ensemble of parameters around the MLE. We calculate the mean-squared errors of the testing set structures using the parameter ensemble, and construct the objective function and its gradients. Next, we apply a conjugate-gradient method to optimize the objective function and obtain the optimal weights of the fitting database; we can also determine if structures should be added or removed from the fitting database. This step can take advantage of structural searching methods,Feng et al. (2008) for example, to identify candidate structures, though we do not do so here. We repeat the circuit with the modified relative weight set of the fitting database until the optimal weights converge.

The testing set is the key component of this approach not only because the objective function consists of the mean squared errors of the testing set structures, but also the empirical potential predictions for structures in the testing set should have small errors—whether that is known from comparison with DFT calculations or estimated from Bayesian sampling without DFT. With the relative errors in the testing set minimized, any weight deviation from the optimal will result in an increase in relative errors. This means that while we could choose weights to reduce the error of one or several testing set structure property function predictions, it will worsen the predictions of other structures and the trade-off is not worthwhile. Although we are able to optimize the fitting database of the empirical potential models, an optimal fitting database does not guarantee a reliable empirical potential model. The optimization algorithm provides the best possible empirical interatomic potential for a given a fitting database and a given testing set, but it has no judgment on whether the optimal Bayesian errors are acceptable; they can, in fact, be quite large. This can occur if the empirical potential model does not contain the relevant physics to describe the atomic environments in the testing set, which produces reduced transferability. Then, we must—for predictive empirical potential methods—decide to improve the potential model itself to increase transferability or remove structures from the testing set to optimize for reduced transferability.

## V Implementation on Lennard-Jones Potential fitting of Ti

### v.1 Potential form and calculation details

We apply the database optimization algorithm to a simple empirical interatomic potential model, the Lennard-Jones potential. The Lennard-Jones potential is a two-parameter pair potential:

 V2(r;r0,Eb)=⎧⎪⎨⎪⎩4Eb[(r0r)12−(r0r)6]:r≤rcutoff0:r>rcutoff−V2(rcutoff;r0,Eb) (30)

where is the binding energy of a dimer with a separation of . We choose the cutoff radius , and the allowable parameters are .

The DFT data set contains six different crystal structures of Ti and the energy versus volume data for all six structures. The DFT calculations are performed with vaspKresse and Hafner (1993); Kresse and Furthmüller (1996), a plane-wave density functional code. We apply a Ti ultrasoft Vanderbilt type pseudopotential,Vanderbilt (1990), with a plane-wave cutoff energy of 400eV for energy convergence of 0.3meV/atom.Hennig et al. (2008) The k-point meshes for different structures are, for hcp, for bcc, for fcc, for hexagonal, for A15 and for , with Methfessel-Paxton smearing parameter of 0.2eV to obtain an energy accuracy of 1meV/atom.Hennig et al. (2005, 2008) The energy versus volume data includes four different structures with volume of the unit cell as and , where is the unit cell volume of the equilibrium structure. The fitting databases are built from various energy differences and energy versus volume data combinations among the six crystal structures. We generate the Markov chain of the potential parameters using the Metropolis-Hasting algorithm.Givens and Hoeting (2005) The ensemble of the potential parameters contains independent parameter sets from the MCMC simulation with attempted steps, with an auto-correlation length of approximately 100. We use a reweighting scheme discussed in Appendix A to approximate the objective function values for all possible sets of weights with only one sampling run. Since a radial potential model does not describe the physics of metallic bonding, we expect that the Lennard-Jones potential will not be transferable for testing sets containing many different structures. Our goal is for the algorithm to identify this lack of transferability in the optimization. We systematically consider different types of fitting databases and testing sets with this in mind.

### v.2 Two structured fitting database

We start with a simple fitting database that contains two structural energy differences, and a testing set with the same structures. Fig. 2 shows two typical objective function behaviors considering all possible weight combinations of two different two-structured fitting database and testing set combinations. Fig. 2 shows the behavior of the objective function of a fitting database with and . The objective function has a unique minimum with a specific relative weight ratio of the two structures. Moreover, if we calculate the derivative of the objective function with respect to weight at endpoints (where one weight is zero), we can see that each derivative of the objective function with respect to weights indicates that the other structure should be added to the fitting database. Therefore the optimal weight value for both fitting database structures are positive, and we refer this as a “mixed” fitting database. On the other hand, Fig. 2 shows a fitting database with and . The objective function has minima at the endpoints, which means that a fitting database containing both and has higher relative errors for the testing set than a “pathological” fitting databases with only one structure. This is due to the non-transferability between hcp and fcc structures. We refer to these pathological cases as “unmixed” fitting databases.

Fig. 3 shows the result of optimizing all possible combinations of two-structured fitting databases with two energy differences sharing a common reference structure. Databases with physical MLEs with positive and are either mixed or unmixed two-structured fitting databases. Most mixed fitting databases include fcc, bcc, hex and A15 structures and most unmixed fitting databases includes hcp or energy differences. By exploring a wide phase space (six crystal structures of Ti) of Lennard-Jones potential fitting, we have shown that the database optimization algorithm offers an automated, systematic and quantitative way of analyzing empirical potential model fitting with different fitting and testing sets.

### v.3 Three-structured fitting database

We next apply the database optimization algorithm to the three-structured fitting databases where the testing set contains the same structures as the fitting database. Fig. 4 is the Gibbs triangle (so that ) contour plot for the objective function for a fitting database constructed from , and . All three of the simpler two-structured fitting databases contained in the three-structured fitting database are mixed fitting databases, and there exists a global minimum in the interior of the Gibbs triangle. The optimal weights for all three fitting database structures are positive, with an optimal weight set , , . If we start with any two-structured fitting database and consider the optimal weight value for the third structure, the gradients of the objective function for the third structure are negative. This indicates that the inclusion of the third structure will decrease the relative errors. Fig. 5 shows the comparison of the prediction distributions with equal weights and optimal weights, where the optimal weight set provides reasonable predictions for all three testing set structures. Note, however, that we are only testing energies at one volume for each structure; we will consider the addition of volume changes in Section V.4.

Fig. 6 is a Gibbs triangle for a three-structured fitting database constructed from , , and . While all three of the two-structured fitting databases are mixed fitting databases, the minimum occurs between and . The optimal weight values are and . The gradient of from the two-structured fitting database and is positive meaning fitting to will increase the relative errors in the testing set. Fig. 7 shows comparison of the prediction distributions evaluated at equal weights and optimal weights. The prediction distribution shows that the optimal testing errors for all three structures are all about 5meV. Therefore, it suggests that for Lennard-Jones potential, one can fit bcc and hex structures to predict A15 structure well.

Fig. 8 is a Gibbs triangle for a three-structured fitting database constructed from , and . Now, all three of the two-structured fitting databases are unmixed fitting databases. Based on the Gibbs triangle contour plot of the objective function, the optimal database includes only one structure . Adding any of our candidate structures to this will increase the relative errors and worsen the predictions. Fig. 9 shows comparison of the prediction distributions evaluated at initial weight and optimal weight. The potential yields a very good prediction for , but poor estimates for the other two structures. If we use this optimal Lennard-Jones potential to predict and , the optimal distributions show that the probability of getting the true values are very low. For , the true value is 58meV with Bayesian errors of 93meV, and for , the true value is 192meV with Bayesian error of 30.8meV. It reveals that the optimal Lennard-Jones potential is not transferrable for this testing set, as we expect.

### v.4 Structural energy differences and volume changes

Next we apply the algorithm to larger fitting database and testing set combinations. The testing set includes the energy differences between hcp and the other five structures and hcp energy versus volume data. We use four hcp structures with unit-cell volumes of and , for the hcp equilibrium volumes . The fitting database starts with the same set of structures and the hcp energy versus volume data. For the four structures representing the hcp energy versus volume data, we constrain their weights to be equal during weight optimization. The and hcp energy versus volume data have optimal weight values: and , and all other weights are zero. Fig. 10 shows that after weight optimization, the inclusion of hcp energy versus volume data improves the prediction of the shape of the energy versus volume curve. It also shows the automatic removal of structures from the database by the optimization algorithm.

We next expand our testing set to include energy versus volume data for all six structures along with the structural energy differences. Our fitting database starts with all structural energy differences, and the hcp energy versus volume data—but not the other energy versus volume data. Now, the energy difference and hcp energy versus volume data have optimal weight values and . In Fig. 11, after weight optimization, the predictions for hcp energy versus volume data improve significantly compared to the initial equal weight guess. Fig. 12 shows that the optimal fitting database offers a close prediction of the shape of fcc energy versus volume curve, which is expected since the fcc and hcp structures have the same first nearest neighbor atoms. Bayesian error of the four bcc energy differences is too large to have a good predictions for the bcc energy versus volume curve. Similarly, in Fig. 13, sloppy predictions for hex and energy versus volume data are obtained from the optimized fitting database. Hence, the optimal Lennard-Jones potential with the given fitting database and testing set combination does not have enough flexibility to predict energy versus volume curve for bcc, hex and .

### v.5 Defect structures without DFT data

Finally, we demonstrate the inclusion of a structure without DFT calculations in the testing set. As we explained in Section III, we assign the structure property prediction as the mean value of the ensemble of the prediction. The “unknown” structure we add in the testing set is an unrelaxed hcp supercell containing one vacancy, and the structure property function is the vacancy formation energy. As we are not including DFT data, the comparison is to the mean value of the vacancy formation energy from our empirical potential. The fitting database includes , and and hcp energy versus volume data. The testing set consists of all hcp energy difference, all six energy versus volume data and the single hcp vacancy configuration. Fig. 14 shows the prediction distributions of the vacancy formation energy before optimization and after optimization. The energy difference and hcp energy versus volume data have optimal weight values: and . The result shows a significant variance reduction for the vacancy formation energy prediction. It suggests that the prediction will be accurate if the DFT calculation locates within the high likelihood parameter neighborhood of the empirical potential prediction.

## Vi Summary

We combine conventional potential parameter optimization methods and the Bayesian sampling technique to propose a new definition of optimal fitting database of empirical interatomic potential models. We choose an objective function as a function of prediction errors for the testing set, and show that minimizing the objective function is equivalent to minimizing the sum of relative errors in the testing set. We optimize the relative weights in the fitting database to minimize the objective function and quantitatively determine the inclusion and removal of structures in the fitting database. Moreover, the performances of two different fitting databases applied on the same testing set can be compared. The algorithm is demonstrated by a simple classical potential model, Lennard-Jones potential fitting of Ti. We go through all possible combinations of two-structured and three-structured fitting databases and analyze the behavior of the objective function with respect to weight change. The new algorithm leads to the best possible empirical interatomic potential model based on the current fitting database and testing set combination.

## Vii Acknowledgements

The authors thank Yuguo Chen, Henry Wu, and Michael Fellinger for useful comments. This work was supported by the Office of Naval Research through ONR Award No. N000141210752.

## Appendix A Reweighting of the sampling chain

When weights are changed in the database, we need to reevaluate the mean in Eqn. (13) using Eqn. (21). A change to the weights in the fitting database changes the likelihood function to . We rewrite as

 ⟨A(θ)⟩F∗=∫P(θ;F∗)P(θ;F)A(θ)P(θ;F)dθ∫P(θ;F∗)P(θ;F)P(θ;F)dθ=∫L(θ;F∗)L(θ;F)A(θ)P(θ;F)dθ∫L(θ;F∗)L(θ;F)P(θ;F)dθ≈∑Ni=1A(θi)L(θi;F∗)L(θi;F)∑Ni=1L(θi;F∗)L(θi;F). (31)

Thus a reweighting term is assigned to the original data, and provides new predictions without requiring a new sampling chain.

## References

You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters