A KIM-compliant potfit for fitting sloppy interatomic potentials: application to the EDIP model for silicon

Fitted interatomic potentials are widely used in atomistic simulations thanks to their ability to compute the energy and forces on atoms quickly. However, the simulation results crucially depend on the quality of the potential being used. Force matching is a method aimed at constructing reliable and transferable interatomic potentials by matching the forces computed by the potential as closely as possible, with those obtained from first principles calculations. The potfit program is an implementation of the force-matching method that optimizes the potential parameters using a global minimization algorithm followed by a local minimization polish. We extended potfit in two ways. First, we adapted the code to be compliant with the KIM Application Programming Interface (API) standard (part of the Knowledgebase of Interatomic Models project). This makes it possible to use potfit to fit many KIM potential models, not just those prebuilt into the potfit code. Second, we incorporated the geodesic Levenberg–Marquardt (LM) minimization algorithm into potfit as a new local minimization algorithm. The extended potfit was tested by generating a training set using the KIM environment-dependent interatomic potential (EDIP) model for silicon and using potfit to recover the potential parameters from different initial guesses. The results show that EDIP is a ‘sloppy model’ in the sense that its predictions are insensitive to some of its parameters, which makes fitting more difficult. We find that the geodesic LM algorithm is particularly efficient for this case. The extended potfit code is the first step in developing a KIM-based fitting framework for interatomic potentials for bulk and two-dimensional materials. The code is available for download via https://www.potfit.net.


Introduction
An interatomic potential (IP) is a model for approximating the quantum-mechanical interaction of electrons and nuclei in a material through a parameterized functional form that depends only on the positions of the nuclei.IPs such as the Lennard-Jones [1,2] and Morse [3] potentials were initially introduced as part of theoretical studies of material behavior in the early 20th century.Interest in IPs renewed in the 1960s following the development of Monte Carlo and molecular dynamics (MD) simulation methods.This led to the development of IPs for a variety of material systems, a trend that accelerated in the 1980s with the formulation of a large number of increasingly complex IPs.(See [4] for more on the history and scientific development of IPs and MD.)Even complex IPs are far less costly to compute than a first principles approach that involves solving the Schrödinger equation for the quantum system.This enables simulations of large systems (millions, billions and even trillions of atoms depending on the IP) for very short times or moderate-size systems for longer times (tens or hundreds of nanoseconds).Such simulations can tackle problems that are inaccessible to quantum calculations, such as plastic deformation, fracture, atomic diffusion, and phase transformations [5].
Obtaining an accurate and reliable IP constitutes a challenging problem that involves choosing an appropriate functional form and adjusting the parameters to optimally reproduce a relevant training set of experimental and/or quantum data.Traditionally IPs were fitted to reproduce a set of material properties considered important for a given application, such as the cohesive energy, equilibrium lattice constant, and elastic constants of a given crystal phase.However experience has shown that the transferability of such IPs (i.e.their ability to accurately predict behavior that they were not fitted to reproduce) can be limited due to the small number of atomic configurations sampled in the training set (although recent work [6] has shown that this approach can be effective in some cases).Further, as the complexity of IPs increases (both in terms of the functional forms and the number of parameters) it can be difficult to obtain a sufficient number of material properties for the training set.This is particularly true for multispecies systems like intermetallic alloys.
To address these difficulties, Ercolessi and Adams [7] proposed the force-matching method in which a training set containing material properties is augmented with the forces on atoms obtained by first principles calculations for a set of atomic configurations.These can be configurations associated with important structures or simply random snapshots of the crystal as the atoms oscillate at finite temperature.By fitting to this information the transferability of the potential is likely enhanced since it is exposed to a larger cross-section of configuration space.The issue of insufficient training data is also resolved since as many configurations as needed can be easily generated.This makes it possible to increase the number of parameters and in fact in the original Ercolessi-Adams potential for aluminum [7] the functional forms were taken to be cubic splines with the spline knots serving as parameters.This gives maximum freedom to the fitting process (at the expense of a clear physical motivation for the functional form).The force-matching method has been widely used since its introduction and there are a number of open source implementations available.The ForceFit [8] and ForceBalance [9] packages target organic-chemistry applications, while potfit [10][11][12] focuses on solid-state physics.
IP fitting programs (such as those mentioned above) define a "cost function" quantifying the difference between the training set data and IP predictions and use global and/or local minimization algorithms to reduce the cost as much as possible.For example, the potfit program uses simulated annealing for global minimization followed by a local polish using a conjugate gradient method.A difficulty associated with this procedure is that IPs are nonlinear functions that are often "sloppy" in the sense that their predictions are insensitive to certain combinations of their parameters [13].These soft modes in parameter space can cause the minimization algorithms to fail to converge.Recently an understanding of sloppy models based on ideas from differential geometry has led to efficient methods for fitting such models [14] The basic idea is that the parameters of a model (like an IP) define a manifold in the space of data predictions.Fitting the IP then corresponds to finding the point on the manifold closest to the training set.Using these ideas Transtrum et al. [14] augmented the Levenberg-Marquardt algorithm with a geodesic acceleration adjustment to improve convergence.The new geodesic Levenberg-Marquardt (LM) algorithm is likely to be more efficient and more likely to converge for sloppy model systems than conventional approaches.Geodesic LM has been applied to a variety of applications in physics and biology, but until now not to the IP fitting problem.
Once an IP is developed it must be integrated into an MD code (or other atomic-scale simulation package).Traditionally this was done on a code-by-code basis.This creates a barrier for the universal usability of IPs.Users are typically limited by the IPs that happen to be available in the package that they choose.The same applies to IP fitting programs, which support a limited set of functional forms that they can fit.Extending codes to support new IPs is a time consuming and error prone process.It also leads to an inevitable uncontrolled mutation of IPs as modifications to functional forms and/or parameters are introduced (often without adequate documentation).More recently, emphasis have been placed on archival storage of IPs within projects such as the NIST Interatomic Potential Repository (IPR) [15] and the Open Knowledgebase of Interatomic Models (KIM) project [16][17][18].The NIST IPR stores tabulated parameters for a limited number of functional forms (mostly the embeddedatom method (EAM) [19][20][21]).OpenKIM hosted at https://openkim.org is a more general effort in which the IP implementation itself is stored in a format that conforms to an application programming interface (API) standard developed as part of the KIM project [22].The KIM API enables any IP stored in KIM to work seamlessly with any KIM-compliant simulation code including ASAP [23], DL POLY [24,25], GULP [26,27], IMD [28][29][30], LAMMPS [31,32], libAtoms+QUIP [33] and Quasicontinuum [34].In addition IPs stored in OpenKIM are exhaustively tested by computing their prediction for a large number of material properties.These calculations are performed by "KIM Tests" also archived in OpenKIM.Users can extend the OpenKIM system by uploading their own Tests.The predictions of the IPs can be explored using text-based searches and a user-extendible visualization system [35].
In this paper we describe a KIM-compliant implementation of the potfit program designed for fitting sloppy IPs.The work has two main objectives.First, we extend the potfit program to support the KIM API making it possible to fit IPs with arbitrary functional forms that are portable to a large number of KIM-compliant simulation codes.Second, we explore the efficiency of the geodesic LM algorithm for fitting IPs using the forcematching method within the potfit framework.We test the new implementation on a KIMcompliant implementation of the environment-dependent interatomic potential (EDIP) for silicon [36][37][38].We use the EDIP model to generate a training set of force configurations and explore the ability of potfit to converge back to the original EDIP parameters from different initial guesses.We find that EDIP has the properties of a sloppy model and that the geodesic LM algorithm can converge to the correct results in situations where other minimizers in potfit fail.
The paper is organized as follows.In sections 2 and 3, we briefly describe the potfit program and the KIM project, respectively.This is followed by a discussion of the modifications that have been made to potfit in order to support the KIM API and the geodesic LM algorithm.In section 5, we present the application of the extended potfit to the EDIP model.We conclude in section 6 with a summary and a discussion of future directions.

The potfit program for fitting interatomic potentials
The potfit program [10][11][12] (http://www.potfit.net) is an open source implementation of the force-matching method.It optimizes an IP's parameters by minimizing a cost function defined as the weighted sum of squares of deviations between calculated and reference quantities for forces, energies and/or stresses.(The particular form of the cost function used in this paper is described in section 5.) potfit consists of two largely independent components: IP implementations, and least-squares minimization.This program architecture makes it comparatively easy to add either an IP model or a minimization algorithm, which facilitated both aspects of the current work.
potfit supports both analytic and tabulated IPs.An analytic IP has a fixed functional form with certain free parameters.Consider for example the simplest pair potential form where the total energy of n interacting atoms is given by in which r αβ is the distance between atoms α and β , and φ (r αβ ) is the energy in the bond connecting them.An example of an analytic function, is the Lennard-Jones potential for which where ε and σ are the IP parameters.In contrast, for a tabulated IP, the functional form is arbitrarily defined as an interpolation between sampling points that serve as the IP parameters.For example, for a cubic spline, given a set of distances (r 1 , r 2 , • • • , r N ) and a corresponding set of function values (φ 1 , φ 2 , • • • , φ N ), the functional form in the interval [r i , r i+1 ] is defined as where the coefficients (a i , b i , c i , d i ) are obtained from the conditions that φ (r i ) = φ i and φ (r i+1 ) = φ i+1 (with different conditions at the ends that define the type of cubic spline).
In this scenario the tabulated function values (φ 1 , φ 2 , • • • , φ N ) are the IP parameters.The form of interpolation used (cubic, quartic, quintic, or otherwise) can have a strong influence on the predictions of the IP [39].This is why the interpolation is considered part of the IP by OpenKIM (see section 3) and stored along with the parameters.The potfit program supports various functional forms including a variety of pair potentials, embedded atom method (EAM) potentials [19][20][21], angular-dependent potentials (ADP) [40], bond-order potentials (Tersoff potential [41,42]), and induced dipole potentials (Tangney-Scandolo [43]).Once an IP is selected, different local and global minimization algorithms are available within potfit to find the optimal set of parameters that minimize the cost function.The local minimization algorithm is Powell's method [44], which is a variant of the conjugate gradient (CG) method.Standard CG requires the derivative of the cost function with respect to the IP parameters.The advantage of Powell's method is that it constructs conjugate directions without requiring derivatives.Powell's algorithm begins with an initial guess for the IP parameters and converges to a nearby minimum.For highly nonlinear IPs it is possible that different initial guesses will lead to different solutions.To address this, global minimization algorithms attempt to explore the parameter space to find the deepest minimum.potfit includes simulated annealing and differential evolution in this class of methods.See the original potfit publications [10][11][12] for detailed descriptions of the optimization algorithms.

The Open Knowledgebase of Interatomic Models (KIM) project
The KIM project [16][17][18] is an international effort aimed at improving the reliability and accuracy of molecular simulations using IPs.
A schematic of the OpenKIM cyberinfrastructure is displayed in figure 1.First and foremost OpenKIM provides a repository for IPs at https://openkim.org.Within the OpenKIM system IPs are referred to more generically as models.A KIM Model consists of the computer implementation of the IP along with any parameters.A KIM IP can be a stand-alone model, or a model driver that reads in different parameter files to define different models.All content in the OpenKIM system is archived subject to strict provenance control and can be cited using permanent links.For example, later in this article we will be analyzing the EDIP model for silicon [36][37][38] archived in OpenKIM [16,45,46].The citations contain a unique 12-digit KIM identifier and 3-digit version number.This makes it possible to access the actual IP used in this publication at any later date to reproduce the calculations -an ability lacking prior to OpenKIM archiving.All content within the OpenKIM repository is citeable in this manner and is accessible to external users through web queries.
In addition to archiving IPs, OpenKIM is tasked with determining the suitability of an IP for a given application.Since IPs are approximate, their transferability to phenomena they were not fitted to reproduce is limited.To help define a rigorous method for IP selection, OpenKIM includes a user-extendible suite of KIM Tests for characterizing an IP.Each KIM Test computes the prediction of an IP for a well-defined material property (such as the surface free energy for a given crystal face under specified temperature and pressure).Each IP in the system is mated with all compatible KIM Tests using an automated system called the Processing Pipeline [47].Whenever a new Model or Test is uploaded to OpenKIM by members of the community, the pipeline launches to perform the new calculations.The results are stored in the OpenKIM repository and can be explored and compared with first principles and experimental reference data using text-based and graphical visualization tools [35].Machine learning based tools are being developed to use the accumulated information on the accuracy of an IPs predictions to help select for applications of choice [48].
The operation of the pipeline involves the coupling of pairs of computer codes, Model and Test, which may be written in different computer languages and need to exchange information.Specifically, the Test repeatedly calls the Model with new input (typically positions and species of atoms) and requires the generated output (typically energy and forces on atoms).This is made possible because all codes in OpenKIM conform to the KIM API [22].This is a standard developed as part of the KIM project for atomistic simulation codes to exchange information.The KIM API is a lightweight efficient interface where only pointers are exchanged between the interfacing codes, which can be written in any supported language pointers in API OpenKIM Repository   It allows researchers to select the appropriate IP for their application and use it with their code of choice, rather than being limited to the IPs that happen to come with it.

Extensions to potfit: Geodesic LM and KIM-compliance
An extended version of the potfit code has been developed that includes a geodesic LM optimizer and support for the KIM API.The new code is available via https://www.potfit.net.Details of the extensions are given below.

Geodesic Levenberg-Marquardt (LM) algorithm
As explained in the introduction, IPs often have the characteristics of sloppy models.This can cause the standard minimization algorithms in potfit to fail to converge.As an alternative we have incorporated into potfit the newly developed geodesic LM algorithm [14,[49][50][51].
Geodesic LM is an efficient method offering a higher likelihood of convergence for sloppy models.The performance of the geodesic LM method is discussed in section 5.4.The method is described next.
The least-squares optimization problem is to find the parameters θ = {θ 1 , . . ., θ N } that minimize the cost function C : R N → R: where the residual r : R N → R M is an M-dimensional vector function of N parameters.The gradient of C is where J is the Jacobian of r, and the Hessian of C is, where K is the second derivative of r with respect to θ (a third-order matrix).Before introducing geodesic LM, let us consider the more basic Gauss-Newton method.This approach is based on a linear approximation of r [52], Substituting ( 4) into (1) we see that, Taking the first derivative of ( 5) with respect to δ θ and setting it to zero, we have, In ( 6) δ θ is a local decent direction of C, thus we can obtain the local minimum of C by solving (6) iteratively.This is the Gauss-Newton method.
Levenberg [53] suggested a "damped version" of ( 6) where the parameters are iteratively updated according to Here λ is a (non-negative) damping parameter, and D = I is the identity.The Levenberg method is an interpolation between the Gauss-Newton algorithm and the steepest decent algorithm.For small λ , (7) reduces to the Gauss-Newton equation in (6); whereas for large λ , δ θ ≈ −λ −1 J T r, which lies along the gradient (i.e. the steepest descent direction) in (2).
In order to overcome slow convergence in directions of small gradients, Marquardt [54] proposed taking D = diag(J T J ), where diag(A) returns a diagonal matrix with elements equal to the diagonal elements of the square matrix A).The downside of this approach is that it greatly increases the susceptibility for parameter evaporation [50], i.e. the algorithm pushes the parameters to infinite values without finding a good fit around saddle points.
To improve the efficiency of the Levenberg-Marquardt (LM) method, Transtrum and Sethna [14,[49][50][51] proposed the geodesic LM algorithm in which the minimization step size is constrained and the residual is modified to include a second-order (geodesic acceleration) correction, The expansion of the cost function is then where a Lagrange multiplier enforces the constraint on the step size δ θ.Taking the first derivative of (9) with respect to δ θ and setting it to zero, we have where the indices are included explicitly to avoid ambiguity and Einstein's summation convention over repeated indices is adopted.The step size δ θ i is written as a sum of two parts: where δ θ (1) i is given by ( 7) and δ θ (2) i includes the remaining terms.Substituting ( 11) and ( 7) into (10), δ θ (2) i can be solved as, where some small terms are ignored (see [50] for details) and the directional second derivative r p = K pkl δ θ k δ θ (1) l is used.As seen from ( 12), the geodesic acceleration correction only depends on the directional second derivative oriented along the first order correction δ θ (1) .This feature is very important because the cost to compute the directional second derivative is reasonably small and hence will not overly add to the computational burden.
Finally, in order to prevent parameter evaporation in the geodesic LM algorithm and to increase the likelihood of convergence, allowable step sizes must satisfy: 2 δ θ (2)  δ θ (1) ≤ α, where α is some parameter (usually smaller than 1) and whose optimal value depends on the specific problem.
An open source implementation of the geodesic LM algorithm has been made available by Transtrum [55].In addition to the geodesic acceleration correction, the geodesic LM package includes different options for updates to the damping parameter and the damping matrix, different conditions for accepting uphill steps (i.e.steps that increase C), and whether or not to use Broyden's method [56] to lower the Jacobian update frequency.In this work, the damping matrix is set to I, Umrigar and Nightingale's method [50] is used to update the damping parameter and to decide whether to accept uphill steps, and we do not employ Broyden's method to update the Jacobian because the computational cost in this problem is not too high.

KIM-compliance
The original version of potfit is limited to using IPs built into the program.The philosophy underlying the KIM API is different.It is based on the idea of separation between the simulation code (called a "simulator" in KIM parlance) and the IP (called the KIM Model).Whenever the simulator (potfit in this case) requires the energy or forces for an atomic configuration, it calls a KIM Model with the necessary input (e.g.coordinates and species of the atoms) and receives back the required output (e.g. the forces on the atoms).Enforcing this separation makes it possible for the simulator to work with any compatible KIM Model in plug-and-play fashion.§ The KIM API is the standard for the exchange of information between the communicating codes.In practice this works by having all data exchanged between the simulator and IP stored in a data structure called the "KIM API Object."The simulator and IP can access and change this data using KIM API library routines.The only direct data transfer required between the programs is for the simulator to pass to the KIM Model the pointer to the location in memory of the KIM API Object.This makes the KIM API extremely efficient and provides it with cross language compatibility (as described in section 3).This means that potfit, which is written in C, can work efficiently and seamlessly with KIM Models written in C, C++, Fortran 77/90/95/2003, and Python.
The potfit program has an unusual requirement for a simulator in that it needs to modify the parameters of the IP that it is optimizing.This is supported by the KIM API through the ability of KIM Models to "publish their parameters."A KIM Model that does so has in its KIM descriptor file a list of its parameters identified by name with the prefix "PARAM FREE ".The potfit user can select any of these parameters to optimize.This is specified in the input file to potfit along with the initial guesses for the parameters, and the KIM Model identifier (see section 3).
The operation of the KIM-compliant implementation of potfit is described in figure 2. The program reads in the input defining the optimization problem and then generates one or § A KIM Model is compatible with a simulator if it supports the required species and boundary conditions that the simulator requires.This is automatically determined by the KIM API using metadata in "KIM descriptor" files that accompany the simulator and KIM Model.When the optimization is complete, potfit outputs the final set of parameters and writes out a new KIM Model, which can be uploaded to https://openkim.orgif desired.

Application of KIM-potfit to EDIP
In this section we test the performance of KIM-potfit by applying it to the optimization of the EDIP model for silicon [36][37][38] archived in OpenKIM [16,45,46].The training set consists of the energy and forces computed with the original EDIP model for a periodic configuration of 5 × 5 × 5 conventional unit cells for a total of 1000 atoms perturbed off the ideal diamond structure positions.Our objective is to see whether starting from different initial guesses for the parameters, the optimization procedure is able to recover the EDIP parameters.

Cost function
As explained in section 2, IP optimization corresponds to minimizing a cost function quantifying the difference between the IP's predictions and a training set of reference data.
Our training set consists of a single atomic configuration with n = 1000 atoms.The cost function used in our analysis is where f α = f α (θ) and f 0 α are the predicted and reference force on atom α, E = E(θ) and E 0 are the predicted and reference energy, and w 1 and w 2 are weights.In our example we take w 1 = 1 ( Å/eV) 2 , w 2 = 1 (1/eV) 2 to give a unitless cost function.The reference forces and energy are defined as f 0 α = f α (θ 0 ) and E 0 = E(θ 0 ), where θ 0 are the original EDIP parameters.EDIP's functional form and parameters along with the training set used are given in the supplementary information accompanying this article.
The Hessian of the cost function at the original parameters θ 0 follows from (3) as where we have used (f α − f 0 α ) θ 0 = 0 and (E − E 0 ) θ 0 = 0. We will use the Hessian in the next section to evaluate the sensitivity of EDIP predictions to its parameters.

EDIP sensitivity analysis
It is a common feature of IPs that the prediction of a model (including the cost function) are weakly dependent on certain directions in parameter space, i.e. the models are sloppy [57].This is demonstrated schematically in figure 3 showing the contour plot of a cost function with only two parameters.Clearly the function varies more slowly along the Λ 2 eigendirection than along Λ 1 .As a result of this structure of C(θ), convergence along certain directions can be incomplete.For example in the figure, the minimum is located at O but a convergence could terminate at a set of parameters offset to the optimal ones by ∆θ.An indication of whether this is occurring can be obtained by computing the angle between ∆θ and the eigenvectors of the Hessian ∂ 2 C/∂ θ i ∂ θ j evaluated at O, Point A and point B (corresponding to another fit) are equally close to the eigendirections because they form the same angles.
In general, the eigendirections associated with smaller eigenvalues are more "sloppy" and we would expect larger scatter in the cos α i values computed in (16) for an ensemble of fits.Thus qualitatively we expect a behavior similar to that shown in figure 4, where the scatter in cos α i obtained from a large number of fitting attempts with different initial guesses would occupy a region similar to Region I in the figure with larger scatter along directions associated with smaller eigenvalues.This would be the signature of a sloppy model.
We perform this analysis for EDIP.EDIP has 13 parameters, but only 11 of them are used in this study with the two cutoffs excluded.The training set is described in section 5.1.We take an ensemble of 100 initial guesses obtained by adding to the original EDIP parameters random numbers ρ drawn from a Gaussian distribution with average zero and a standard deviation of 0.15, i.e. θ i := θ 0 i × (1 + ρ), i = 1, . . ., 11, and perform the fitting procedure for each one.The relative errors ∆θ i /θ 0 i = (θ * i − θ 0 i )/θ 0 i (for i = 1, . . ., 11) between the parameters θ * obtained in the fitting procedure and the original EDIP parameters θ 0 for all 100 optimizations are shown in figure 5a.No clear trend is apparent.Instead let us plot the cosine defined in (16) as done schematically in figure 4. To improve the scaling of the plot let : Contour plot of a cost function with two parameters.Λ 1 and Λ 2 are the eigendirections of the Hessian of the cost function corresponding to the larger and smaller eigenvalues, respectively.Points A and B represent two fits.
Figure 4: Expected region for normalized projected parameter errors if the fits are aligned along the sloppy direction(s).The horizontal axis corresponds to the eigenvalues of the Hessian sorted in increasing order.us modify the definition to use the logarithm of the parameters, θi = log θ i .(This is important when different parameters have very different magnitudes.)The definition equivalent to ( 16) is then where ∆ θ = θ * − θ0 and Λi is the ith eigendirection of the Hessian of the cost function C in the logarithm parameter space, The results for EDIP are plotted in figure 5b, and are very much in line with the expected behavior in figure 4. In particular, we see that the spread is mostly accounted for by the first four eigenterms each of which corresponds to a particular combination of parameters as quantified by the eigendirections.

Cost along eigendirections
The results in the previous section support the conjecture that the observed scatter in results is due to sloppy directions.However, since the Hessian is a local measure of the structure of C near the optimal solution, it is possible that there are other minima further away that can trap the optimization process.To explore this possibility, we plot the cost function along the different eigendirections, where as before θ 0 corresponds to the original EDIP parameters and s is a scalar factor that can be positive or negative.The results are plotted in figure 6 for all eigendirections from the smallest eigenvalue 1 to the largest eigenvalue 11.The cost C varies more slowly along the eigendirections associated with the five smallest eigenvalues than along other directions, i.e. the cost function is more flat in these five directions.(Note that the scale of the y-axis is not the same in the plots, and the first five plots have a significantly smaller range than the others.)This is in agreement with the previous observation that the spread of the fitted parameters is larger along these directions.We also see that an additional local minimum (aside from the one associated with original parameters at s = 0) is only apparent along the positive 11 direction at a significantly higher cost function.This minimum can trap the optimization algorithm for initial guesses in its vicinity, however over a large range near the original parameters there is only a single minimum.This suggests that a local minimizer is sufficient for this IP (using different initial guesses), which can reduce the fitting time drastically compared with global minimization efforts.

Performance of minimization algorithms
Having established in the previous section that a local minimization is sufficient for EDIP (at least in the vicinity of the original parameters), we explore the efficiency of Powell's method, the geodesic LM algorithm added to potfit, and the standard LM algorithm described in [58] We verified that this is a local minimum by checking the eigenvalues of the Hessian at the minimum and finding that they are all positive.(see section 4.1 for details of LM and geodesic LM).The three minimizers are applied to fit EDIP.The initial guesses for the parameters are obtained in the same way as discussed in section 5.2.Three sets of 100 initial guesses are constructed by perturbing the original EDIP parameters by random numbers drawn from a Gaussian distribution with zero average and standard deviations of 0.1, 0.2 and 0.3.The results are given in table 1 showing for each minimizer and perturbation magnitude the number of attempts (out of 100) that converged below the cost level indicated on the left.For example for Powell's method with a perturbation magnitude of 0.2, the algorithm converged to a cost below 10 −3 for 46 out of the 100 initial guesses.We see that geodesic LM outperforms the other two algorithms for all perturbation  amplitudes.In all cases it is able to converge from more initial guess to below the specified tolerance.The standard LM algorithm is superior to Powell in all cases except for the smallest perturbation.
We are also interested in the computational cost of the different minimization methods.The Powell algorithm requires only evaluations of the cost function, whereas LM and geodesic LM require the gradient of the cost function with respect to the parameters.The gradient is computed using numerical differentiation, which only requires evaluations of the cost function.Therefore the number of cost function evaluations (N cost ) performed during the minimization is a good metric for comparing the performance of all three methods.Figure 7 shows for each perturbation amplitude, the normalized cost function C = C/C init (where C init is the cost at the start of the minimization) versus N cost for the three methods.The initial flat region is associated with the setup stage of the methods (e.g.building an initial estimate for the Hessian in Powell's method).Examining the plots, we see that Powell and LM are comparable in performance, whereas geodesic LM is about an order of magnitude slower.Thus geodesic LM's superior convergence properties come at the expense of a significant increase in the computational burden.For the current problem this tradeoff is acceptable since geodesic LM converges in most cases in less than 1000 cost function evaluations, which is quite reasonable for the current problem.However, when computation of the cost function is expensive (due to a larger training set or more complex IP) or when more steps are needed to converge as would be expected for IPs with more parameters, the computational cost of geodesic LM may be an issue.

Summary and future directions
The accuracy of an IP is critical to atomic-scale simulations.Designing a high-quality transferable IP is a very challenging task due to the complexity of the quantum mechanical behavior that the IP is attempting to mimic.The developer of an IP must select a functional form with adjustable parameters and a suitably-weighted training set of first principles and experimental reference data.The IP is then optimized with respect to its parameters to reproduce the training set as closely as possible in terms of specified cost function.It is then validated against an independent test set of reference data.If the validation results are inadequate, the functional form and/or the training set must be amended and the process repeated until a satisfactory IP is obtained.This process is made more difficult by the fact that many IPs are sloppy models in the sense that their predictions are weakly dependent on certain combinations of their parameters.This can cause optimization algorithms to fail to converge.
potfit is an IP fitting program based on the force-matching method.The idea is to include in the training set the forces on atoms computed from first principles for a set of atomic configurations of interest.This allows the developer to train the IP with a large amount of reference data spanning a larger portion of configuration space.We have extended potfit in two ways: (i) We have incorporated a new optimization method based on the geodesic LM algorithm.
Although slower than some other local minimizers it has improved convergence properties for sloppy models.This is the first application of geodesic LM to IP optimization.(ii) We have adapted potfit to be compatible with the KIM API.This allows potfit to optimize KIM Models stored in the KIM Repository at https://openkim.org.This greatly expands the functional forms that potfit can optimize.The resulting IPs are portable and can be used seamlessly with many major software packages that support the KIM API.
We study the effectiveness of the new code by using it to optimize the EDIP model for silicon.This IP is not available within the original potfit code, but is available as a KIM Model.
We use EDIP to generate a training set of atomic configurations and test whether potfit can reproduce the EDIP parameters from different initial guesses.We make several observations: (i) EDIP displays the characteristics of a sloppy model.For different initial guesses we find a large scatter in the parameter values.An analysis shows that the scatter is maximal along the eigendirections of the Hessian of the cost function associated with the smallest eigenvalues.Thus the IP's behavior depends primarily on a few stiff combinations of parameters.
(ii) At least in the vicinity of the original parameters the cost function is largely convex with only two local minima.This suggests that a local minimization method with different random initial guesses is sufficient.
(iii) The geodesic LM algorithm was on average twice as likely to converge to the correct solution from different initial guesses than standard LM or Powell's method that tended to get trapped along sloppy directions.However, the improved performance comes at the cost of about an order of magnitude increase in computational expense for geodesic LM.
For the current problem the computations were not prohibitive and geodesic LM was the preferred method.
The extended potfit code is the first step in developing an IP fitting framework based on the OpenKIM project.Such a framework would make it possible for the developer to draw upon the large selection of KIM Tests (standardized codes for computing material properties) and reference data in OpenKIM.Models developed using this framework and the training set used in the fit will then be archived at https://openkim.org.The latter is an important advantage as it means that other developers will be able to draw upon previous work and adapt a Model to a new application starting from the training set used in its development.This framework is currently under development for use in fitting IPs for 2D materials [59].

Training set
The training set consists of the energy and forces computed with EDIP for a periodic configuration of 1000 atoms perturbed off the ideal diamond structure positions.The training set file "supplementary traning set.txt"accompanies this article.The header of the file specifies the supercell lattice vectors and the average energy per atom.This is followed by a list giving for each atom the atomic species code, atom coordinates, and force on the atom.See the potfit webpage http://www.potfit.net/wiki/doku.phpfor details on the format of the training set file.

Figure 2 :
Figure 2: Flowchart of the KIM-compliant potfit program

Figure 5 :
Figure 5: Errors in fitting the EDIP parameters.(a) Relative errors.(b) Normalized error projected on the eigendirections of the Hessian of the cost function in the logarithm parameter space.

Figure 6 :
Figure 6: Cost along eigendirections.C 1 is associated with the smallest eigenvalue, and C 11 is associated with the largest.For C 10 and C 11 the insets show a zoomed view of the flat portions of the curves.

Figure 7 :
Figure 7: Normalized cost C versus the number of cost function evaluations N cost .The initial guesses for the parameters are perturbations off the original EDIP parameters drawn from a Gaussian distribution with zero average and standard deviations of (a) 0.1, (b) 0.2, and (c) 0.3.
Fortran 77/90/95/2003, Python).A number of major simulation codes currently support the KIM API (listed in section 1).Anyone of these codes can work seamlessly (without having to be recompiled) with any KIM Model stored on the same computer as a dynamically linked library.(This is akin to a CD being placed in a player from any manufacturer and having its content played.)This greatly increases the portability of IPs.

Table 1 :
Statistics on the optimization of EDIP.The number of minimization attempts (out of a 100) that converge to a cost function below the tolerance specified on the left for the three minimizers studied and the three perturbation amplitudes (0.1, 0.2 and 0.3) used to generate the initial guesses.

Table 1 .
EDIP Parameter values for silicon