Uncertainty quantification for classical effective potentials: an extension to potfit

Effective potentials are an essential ingredient of classical molecular dynamics (MD) simulations. Little is understood of the consequences of representing the complex energy landscape of an atomic configuration by an effective potential or force field containing considerably fewer parameters. The probabilistic potential ensemble method has been implemented in the potfit force matching code. This introduces uncertainty quantification into the interatomic potential generation process. Uncertainties in the effective potential are propagated through MD to obtain uncertainties in quantities of interest, which are a measure of the confidence in the model predictions. We demonstrate the technique using three potentials for nickel: two simple pair potentials, Lennard-Jones and Morse, and a local density dependent embedded atom method (EAM) potential. A potential ensemble fit to density functional theory (DFT) reference data is constructed for each potential to calculate the uncertainties in lattice constants, elastic constants and thermal expansion. We quantitatively illustrate the cases of poor model selection and fit, highlighted by the uncertainties in the quantities calculated. This shows that our method can capture the effects of the error incurred in quantities of interest resulting from the potential generation process without resorting to comparison with experiment or DFT, which is an essential part to assess the predictive power of MD simulations.


Introduction
Our understanding of the physics underlying material properties relies on verification from computational models of materials and molecules. Such materials simulations also allow us to predict new properties and structures which can then be reproduced experimentally. In order to facilitate the modelling, interatomic potentials have long been used to circumvent the limitations in timescale and simulation size of costly first principles calculations by specifying the energy dependant on the atomic positions. This functional representation is a key constituent of molecular dynamics simulations, where the quality of the output relies predominantly on the potential employed. Currently, the systematic error incurred in using an interatomic potential is generally unknown, as is the resulting effect on quantities of interest it is used to predict, therein forming the motivation for this work.
The general idea behind interatomic potentials is that the energy of a collection of atoms can be represented by an explicit functional form or model, dependant on the atomic separation. These analytic potential functions encode the physics into the system and contain a limited number of additional parameters which are adjusted to reproduce desired quantities. The first potentials used intuitive functional forms, fitted to experimental data, however new potentials are frequently fitted to ab-initio data such as atomic forces, energies and stresses.
We use the potfit open source force matching code to fit interatomic potentials to density functional theory data, and subsequently quantify the uncertainty incurred in simulations using the potential, in lieu of a first principles approach.
There has been a significant recent effort in the literature to quantify uncertainty in this multiscale modelling process, with Bayesian frameworks proposed for a variety of interatomic models and force fields [1][2][3]. There has also been work toward the quantification of uncertainty due to the potential fitting reference set [4], as well as a proposed framework to efficiently propagate parameter uncertainties to MD outputs [5]. More specifically, quantification of parameter uncertainty for single potentials has been undertaken in a handful of cases [6][7][8][9]. We add to the growing body of uncertainty quantification (UQ) work in potential development and application by providing an open source implementation of the framework in [1] for use in future potential development projects.
We have implemented a new module in potfit which adds a potential ensemble functionality to the potential fitting work flow. The ensemble of potentials generated by the potential ensemble inform of the uncertainty in the parameters of a potential fitted with potfit. Section 2.1 introduces the potfit code, followed by an outline of the potential ensemble method in section 2.2, with implementation specific details in section 3. We then show in section 4, how an ensemble can be propagated through molecular dynamics simulations to illustrate the incurred uncertainties in both equilibrium and non-equilibrium quantities of interest (QoI): the equilibrium lattice constant, elastic constants C 11 , C 12 and C 44 , and thermal expansion coefficient at 300 K.

The potfit force matching code
The potfit code [10][11][12][13][14] is an open source package implementing the force matching method [15], where the parameters of an interatomic potential are adjusted to optimally reproduce forces, energies and stresses typically obtained from DFT calculations. The potential parameters θ = {θ 1 , ..., θ N } either belong to an analytic potential, or are the values of the potential function at sampling points for tabulated potentials.
In the force matching method, the deviation from the reference data is quantified by the N-dimensional cost function or kernel, C(θ), with (1) The first sum runs over all M force components of the reference configurations in the training set, where F 0 k is the set of individual atomic forces, with weighting a k . Here F k (θ) represents the corresponding set of forces from the potential for each atom in the configuration. The fit can be enhanced with additional information about the target system energies and stresses (and optionally other quantities), represented by A 0 r . The quantities can be obtained through first principles calculations and can be given weights b r , depending on the importance of accuracy in their descriptions given by the potential, A r (θ). The best fit potential parameters θ * produce the minimum cost value, C(θ * ).
The fitting procedure, by default, uses a combination of simulated annealing followed by a gradient descent method to minimise the cost function, although there are a variety of space searching algorithms implemented. Recently potfit was also updated to work within the OpenKIM framework [16][17][18], providing users with easy access to a repository of fitted and tested potential models to utilise and expand.

The potential ensemble method
In the cost landscape defined by the potential fitting procedure, there typically is significant covariance between potential parameters, hence the eigen-directions are used to define the basin curvature. The relative degrees of curvature are given by the eigenvalues of the Hessian at the cost minimum where θ = {θ i } N i=1 represents a set of interatomic potential parameters. Using the information about the curvature of the minimum, we can generate an ensemble of candidate potentials of varying suitability. This ensemble inherently describes the robustness of each parameter fit to the reference data, and hence their uncertainty.
The eigenvalues of the minimal cost space Hessian often have a large spread in their magnitudes, representing a best-fit basin with vastly differing steepness along eigendirections. The hallmark of a sloppy model is that the the basin encapsulating the minimum has significantly differing degrees of curvature along the pricipal axes, therefore the majority of interatomic potentials fitted in potfit fall into the category of sloppy models. Investigation into the sampling of such sloppy models has been extensively undertaken by Brown and Sethna [19], with a focus on their occurrence in systems biology. The approach has since been outlined for interatomic potentials [1] and forms the basis for the implementation of the uncertainty quantification in potfit. The approach relies on generating a potential ensemble to represent the uncertainty, by drawing samples scaled in parameter directions using curvature information from the Hessian in (2).
Markov chain Monte Carlo (MCMC) is used to draw samples from the minimum basin. Candidate steps are generated using random displacements, taking into account information about the curvature from the eigenvalues of the minimum Hessian. In this way, larger steps are taken in sloppy directions (i. e. those associated with small eigenvalues), with smaller steps proposed in stiffer directions.
Steps are taken in cost space, starting from the best fit parameter set, θ * , by proposing a simultaneous perturbation to each parameter of the form where ∆θ i is the proposed change to each potential parameter i, R is a system dependent scaling factor and V i j the matrix of normalized eigenvectors of the Hessian. The parameter λ j is j-th Hessian eigenvalue and r j a normally distributed random number. The acceptance criteria for a MCMC step is set by a temperature, T , where the cost minimum is analogous to sampling at a temperature of T = 0. The sampling temperature is by default chosen to be the natural temperature T 0 = 2C(θ * )/N, which is dependant on the dimensionality of the cost landscape. This follows the assumption that the increase in cost space dimensionality due to each additional parameter requires a reduction of sampling height of the basin to ensure sample errors are captured within the errors of C(θ * ).
The acceptance probability of each the Monte Carlo move is calculated as This ensures downhill moves are always accepted, and that MC moves to higher cost potentials are accepted with a probability decreasing exponentially with the increase in cost between potential parameter sets.

Extensions to potfit: uncertainty quantification
The new functionality in potfit allows for the quantification of uncertainty in any potential fitted using the program, as well as any quantities predicted by it. The uncertainty quantification extension produces a set of accepted MCMC potentials with their attributed cost and weighting. A decorrelated subset forms an ensemble of potentials which can be used to quantify the uncertainty in the potential parameters, and the resulting uncertainties incurred in using the potential to predict quantities of interest. A demonstration using the ensembles for three analytic potentials fitted for nickel to quantify uncertainties in elastic constants and thermal expansion coefficient is detailed in section 4. The uncertainty quantification component requires few external parameters to run. The number of potentials is specified, and the R value in (3) is tuned. It is recommended that an initial investigation trialling a variety of R values is undertaken, with the objective being to accept approximately 23 % of moves for optimal sampling of high dimensional cost landscapes [20]. There is also the option for the user to alter the sampling temperature; T = αT 0 where α is a scale factor. This alteration can be used to improve the convergence of the ensemble to the underlying distribution [1] and is demonstrated in section 4 for the EAM potential. The calculation of the Hessian curvature, which relies upon a finite difference calculation of the cost space minima, also has a tunable perturbation parameter. The percentage perturbation to each parameter in the finite difference calculation can be tuned to ensure the curvature on the scale of the minimum basin is probed. In the case that a new minimum is found at any point in the process, the implementation outputs this new optimal parameter set and restarts the procedure.
The potential ensemble implementation can be used as part of an end-to-end potential fitting workflow, or can be instigated from a previously fitted potfit potential as a stand-alone analysis. Full documentation for the implementation, as well as the code is available on the potfit website ‡.

Demonstration of uncertainty quantification for three analytic models
We evaluate three analytic potential models of increasing complexity and compare their performance and uncertainty against DFT and experimental values. The potentials investigated were the Lennard-Jones (LJ), Morse and embedded atom method (EAM) potentials. A hard cut-off of 10 Å was used for each, with the tails smoothed to converge to zero over 0.75 Å. Two pair potentials (LJ and Morse) were chosen as they contain two and three parameters respectively, thus the ensemble parameter output could be easily visualized in initial investigations. The 10-parameter EAM potential is chosen, with the large parameter space (and encoded physics) expected to best reproduce the desired QoI with the least error. Furthermore, there are a variety of high quality EAM potentials for nickel in the literature [11,[21][22][23]. For the pair potential part of the EAM we have used a Morse function, for the embedding function the universal form proposed by Banerjea and Smith [24] is used, and finally the transfer function is of the oscillatory form necessary for cubic metals as detailed in [25].
The potentials are fitted using potfit to a nickel reference set comprising of 23 atomic configurations of 108-atom fcc (3 × 3 × 3 unit cell) DFT and MD snapshots at a variety of temperatures and stresses. CASTEP [26] was used for the DFT simulations, and MD was performed using LAMMPS [27]. For each MD snapshot in the reference set, a DFT

Lennard-Jones
Morse EAM calculation was performed on the configuration to obtain the forces, energy and stresses used to fit the potentials. Figure 1 shows the spread in the eigenvalues obtained for the three fitted analytic forms. With a difference of up to six orders of magnitude in the degree of curvature as denoted by the magnitudes of the eigenvalues, this illustrates the necessity for a sampling procedure which accounts for this variation in order to efficiently sample the underlying distribution. It is of note that the reference dataset was not particularly tailored to predict the elastic constants or thermal expansion coefficient. When fitting a potential for production use, the process typically involves multiple iterations of the reference dataset and potential fitting to tailor the regions of cost space explored for its intended use. The potentials generated are solely an illustration of the newly implemented potential ensemble method and are by no means suggested as new production-grade potentials to be used outside of this work in the modelling of nickel.
For each of the three fitted analytic potential forms, a 500-member ensemble of potentials were obtained from MCMC samples output from the potfit ensemble implementation. To ensure uncorrelated ensemble members, starting from the best fit potential, each sample was drawn after 50 000 accepted steps with roughly 23 % [20] of steps accepted for each analytic form through tuning of the value of R in (3).
Sampling from the EAM potential landscape was performed at a reduced cost temperature of 0.05 T 0 in accordance with [1], due to the Markov chain leaving the minimum basin for higher temperatures. Outside, the Hessian calculated at the minimum is no longer valid, resulting in inefficient sampling. Reducing the temperature avoids these issues. We believe that this behaviour is due to the quality of the reference data the potential is fitted to. In the fitting of a potential for production use, the reference data is typically weighted and complemented with configurations generated using iterative improvements of the fitted potential, which might alleviate the issue. In case that does not provide a remedy, a different Table 1.
Comparison of results for QoI with their associated uncertainties (IQR). The equilibrium lattice constant is denoted a. C 11 , C 12 , C 44 are the elastic constants in Voigt notation, and α is the linear thermal expansion coefficient at 300 K. The DFT value for a is from a geometry optimized cell included in the reference data.

QoI
DFT sampling strategy would need to be employed.
We have chosen to investigate the performance of the potentials in reproducing the equilibrium lattice constant, the elastic constants C 11 , C 12 and C 44 , and the thermal expansion coefficient at 300 K. Each analytic potential from each of the three ensembles was propagated through MD simulations to obtain the uncertainties for each potential model.
The resulting uncertainties displayed in figures 2 and 3 are compared using box-whisker diagrams, illustrating the uncertainty in each quantity by the inter-quartile range (IQR). The box denotes the IQR, with whiskers extending 1.5×IQR beyond each quartile. The values obtained from the best fit (minimum cost) potential are shown in dashed purple. The notches indicate the confidence interval in the ensemble median (green) and the ensemble means are indicated as dotted red. Due to the presence of skewness and outliers, the uncertainties are reported using the resistant measures of sample median and IQR. Reporting uncertainties via the sample means and standard deviations will bias the reported quantities towards misleading ensemble outliers which tend to result from unlikely higher cost potential ensemble members.

Equilibrium lattice constant
The lattice constants were found by minimising an fcc nickel lattice with a starting lattice parameter guess of 3.5 Å for each potential ensemble member. Table 1 reports the uncertainties for each fitted analytic potential.
The Lennard-Jones ensemble is unable to capture the correct lattice constant within the IQR. This is largely due to higher temperature configurations in the reference data to which the potential was initially fit. The Morse ensemble does manage to capture the correct value, but the larger spread in uncertainty for a simple 0 K quantity again implies the potential is limited by the higher temperature configurations in the reference data. It is important to note that poor potential performance due to reference data selection is not always the culprit, an insufficient choice of model may mean the capturing the correct physics is not possible. The performance of both the Lennard-Jones and Morse clearly illustrate the case of insufficient model choice.
The EAM potential accounts for a non-linear dependency on the local environment through the embedding term. It appears that this is essential for a realistic prediction of the lattice constant given the set of reference data used: The EAM ensemble, albeit sampled with a reduced temperature, not only captures the DFT value within its error but also has a considerably constricted spread of the ensemble predictions compared to the pair potentials.

Elastic constants
The elastic constants are investigated to compare the restoring forces of the potentials to small perturbations of atomic positions. Calculation of the elastic constants for pair potentials is known to be unable to resolve the differences in C 12 and C 44 due to insufficient parameters to describe the off-diagonal tensor components. This is illustrated in figure 2, with only the EAM potential having distinct values for these elastic tensor components. It is noticeable that the best fit potential does not necessarily lie near the centre of the prediction interval. This highlights that there exist competing potentials of comparable suitability which may shift predictions in a particular direction away from the initially obtained best fit value. The performance of the Lennard-Jones potential in the prediction of C 11 and C 12 = C 44 is poor as expected. A very large spread in the predicted values for both, all of which fail to reproduce the experimental values is expected of an ill-fit two parameter potential. This demonstrates the known limitations of such a simple potential and illustrates results in line with an insufficient choice of model.
The Morse potential is able to capture the expected experimental values for both C 11 and C 44 , although still unable to resolve the C 12 , C 44 difference due to the pair potential nature. This is another example of poor model selection despite promising initial predictions of the diagonal elasticity elements.
The EAM potential is able to capture the expected trend for the three constants but does not capture the expected value for C 44 in the uncertainty. Again the initial failure of the fitted EAM potential to achieve the expected off-diagonal C 44 component could be rectified through improvements to the reference data. Figure 2 focuses on demonstrating the output of the potfit uncertainty quantification, for a comparison of production grade EAM potentials in predicting the elastic constants the reader is referred to [30].

Thermal expansion coefficient
Finally the thermal expansion coefficient for nickel at 300 K was calculated. The linear expansion of the solid is investigated to compare the energetic contributions to the system provided by the potentials. The thermal expansion coefficient was calculated by evaluating the length change of crystalline nickel for five temperatures, symmetrically distributed around 300 K at 20 K intervals. A curve is then fitted to the results using regression to find the thermal expansion coefficient at 300 K. On first inspection the results in figure 3 may misleadingly imply that the pair potentials outperform the EAM, due to the small spread in the ensemble members. However upon closer inspection, the predicted pair values in-fact fail to capture the correct value within the IQR and even within the tails of the distribution. The EAM potential does manage to bound the correct value within the uncertainty, although a large spread in the uncertainty is again likely a result of the choice of reference data to which it has been fit.
The uncertainties for the thermal expansion coefficient demonstrate the importance of looking at the predictions of a selection of relevant QoI when evaluating the suitability of a potential model. In table 1 the predicted lattice constant and uncertainty bounds for the Lennard-Jones potential fail to capture the simple equilibrium quantity despite a large uncertainty. This unsuccessful prediction, combined with the high relative uncertainty in the thermal expansion coefficient, illustrate an example of poor model selection.

Conclusion
The potfit potential fitting workflow has been enhanced to generate an ensemble of potentials encapsulating the uncertainty of the correct parameters. This allows a propagation of this uncertainty to quantities of interest of molecular dynamics simulations. The ensemble is generated by sampling the cost landscape using sampling techniques developed for sloppy models. In this work we demonstrated a preliminary uncertainty quantification of quantities of interest for three distinct effective potential models. The implementation enables users to quantify the uncertainty of simulation values incurred by the choice of potential parameters. This is a further puzzle piece towards reproducible and transparent MD simulations -an effective potential should not exist on its own, but rather together with its implementation (as e.g. provided by the OpenKIM framework [14]), its reference data and its uncertainties (potfit + UQ). This integration is also a step toward predictive simulations, i.e. with error bounds determined a priori.
A potential future line of enquiry opened by this work is to investigate how the ensemble information may be used efficiently in production simulations. While ensemble simulations are trivially parallel and scale perfectly with the number of ensemble members, reducing the number of simulations performed may still be desirable, e.g. by intelligently selecting and weighting potentials in the ensemble. Similarly the determined uncertainties could inform the choice of reference data. Closing this feedback loop may lead to further improved classical effective potentials and trustworthy simulation results with quantified uncertainties.