{'fraction_non_alphanumeric': 0.04361841215454668, 'fraction_numerical': 0.015508768766061044, 'mean_word_length': 4.186228664952069, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 18, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A framework for a generalisation analysis of machine-learned interatomic potentials

Machine-learned interatomic potentials (MLIPs) and force fields (i.e. interaction laws for atoms and molecules) are typically trained on limited data-sets that cover only a very small section of the full space of possible input structures. MLIPs are nevertheless capable of making accurate predictions of forces and energies in simulations involving (seemingly) much more complex structures. In this article we propose a framework within which this kind of generalisation can be rigorously understood. As a prototypical example, we apply the framework to the case of simulating point defects in a crystalline solid. Here, we demonstrate how the accuracy of the simulation depends explicitly on the size of the training structures, on the kind of observations (e.g., energies, forces, force constants, virials) to which the model has been fitted, and on the fit accuracy. The new theoretical insights we gain partially justify current best practices in the MLIP literature and in addition suggest a new approach to the collection of training data and the design of loss functions.


Introduction
A key question in machine learning tasks is to understand how well a trained model generalizes to inputs outside of the training data. This is particularly challenging in scientific machine learning where one oftentimes requires generalisation to inputs very far from training data. The present work is concerned with sketching out a multiscale numerical analysis framework suitable to study this scenario.
We will focus in particular on atomistic mechanics simulations using machine learned interatomic potentials (MLIPs) [1,3,4,5,16,36]. The success of molecular simulation relies on the accuracy and efficiency of the interatomic force models. The two main approaches to computing interatomic forces are ab initio electronic structure models [21,23,24,35] and purely mechanistic models [15,20,25,38]. The former are computationally prohibitive while the latter often provide insufficient accuracy. The achievement of MLIPs is to provide a classes of models with tunable accuracy/efficiency ratio which promises to bridge the significant gap in accuracy and capability between ab initio electronic structure models and classical mechanistic models (empirical potentials).
In this work we will study by analytical (as opposed to statistical) methods how the choice of training data and the accuracy of the fit to that training data affect the accuracy of predictions.
The key challenge we hope to better understand is the following: training data for MLIPs is obtained from ab initio electronic structure simulations. Due to the high computational cost of these models, only small computational domains (structures) containing at most hundreds of atoms are used. However, predictions during simulations are performed on much larger domains often containing hundreds of thousand or even millions of atoms.
As a prototype application we consider classes of structures containing crystalline defects, and to keep the notational and technical burden to a minimum we further restrict the present work to point defects only. In this setting, training domains would typically be small cells containing a single defect, while simulations would be performed on much larger domains containing potentially many copies of the defects trained on. In this situation, the basic intution is clear: the trained MLIP has already seen local snapshots of the structure on which it is predicting energies or interatomic forces, and due to its functional form that is indeed restricted to only local interactions it is therefore able to to make accurate predictions.
Yet, the details are subtle and warrant a deeper look: How does the prediction error depend on the size of the training domain? How should the various observations we make (energies, forces, virials) be weighted? Which of these observations provide the dominant contribution to the prediction error? Thus, we see that even this highly simplified setting leads to interesting questions that can significantly inform the design of parameter estimation schemes.
1.1. Outline. We focus on multiple point defects embedded in a periodic homogeneous host crystal, where a rigorous numerical analysis approach is feasible. The atomistic equilibration problem for a single crystalline defect in this context is a well-defined variational problem [9,17]. We review the framework and adapt it to the case of multiple point defects considered in this work along the line of [22] in Section 2. This requires in particular a new existence and stability result (Theorem 2.1) for general configurations of multiple point defects. Our generalisation analysis heavily relies on this result since it characterises the structure of the equilibrium structures.
To propose a framework within which the generalisation can be rigorously understood, we investigate the error propagation from fitting MLIPs on a small training domain to predicting the material properties (e.g., defect geometry and formation energy) on a large simulation domain in Section 3. We demonstrate how the accuracy of the material properties in the simulation depends explicitly on the size of the training structures, on the kind of observations (e.g., energies, forces, force constants) to which the model has been fitted, and on the fit accuracy. Explicit theoretical convergence rates are summarized in Theorem 3.1 and Table 1.
We then propose a concrete implementation of MLIPs inspired by our generalisation analysis to confirm the analytical error estimates on several model problems in Section 4.
Finally, we will discuss further consequences and limitations our work in Section 5. For example, a generalisation to other "simple defects" such as straight dislocation lines appears straightforward. In more complex scenarios, such as curved dislocations or extrapolating on grain boundary structures with distinct coordination environment there are additional challenges that our analysis does not cover even heuristically and requires significant additional ideas.
1.2. Notation. We use the symbol ·, · to denote an abstract duality pairing between a Banach space and its dual space. The symbol | · | normally denotes the Euclidean or Frobenius norm, while · denotes an operator norm. For a finite set A, we will use #A to denote the cardinality of A. For the sake of brevity of notation, we will denote A\{a} by A\a, and {b − a | b ∈ A} by A − a. For E ∈ C 2 (X), the first and second variations are denoted by δE(u), v and The symbol C denotes generic positive constant that may change from one line of an estimate to the next. When estimating rates of decay or convergence, C will always remain independent of the system size, the configuration of the lattice and the the test functions. The dependence of C will be normally clear from the context or stated explicitly. The closed ball with radius r > 0 and center x is denoted by B r (x), or B r if the center is the origin.

Background: Equilibration of crystalline defects
A rigorous framework for modelling the geometric equilibrium of crystalline defects has been developed in [8,9,17,31]. These works formulate the equilibration of a single crystalline defect as a variational problem in a discrete energy space and establish qualitatively sharp far-field decay estimates for the corresponding equilibrium configuration. We will review the framework and adapt it to the case of multiple point defects considered in this work along the lines of [22]. This will provide the analytic foundation of our generalisation analysis. For the sake of simplicity of presentation, we will skip over some technical details but fill these gaps in Section 6.1.
Let d ∈ {2, 3} be the (effective) dimension of the system. A homogeneous crystal reference configuration is given by the Bravais lattice Λ h = AZ d , for some non-singular matrix A ∈ R d×d . We admit only single-species Bravais lattices. There are no conceptual obstacles to generalising our work to multi-lattices, however, the technical details become more involved. The reference configuration with defects is a set Λ ⊂ R d . The mismatch between Λ and Λ h represents possible defected configurations. In this paper, we consider multiple point defects in a finite domain with periodic boundary conditions. To that end, let We denote the continuous cell by Ω N : where Λ N is the periodic computational domain and Λ per N is the periodically repeated domain. We consider n D point defects in Λ N , e.g. vacancies or interstitials, located at i ∈ Ω N for i = 1, . . . , n D . Let D : be a set of the positions of these defect cores in Λ N . We assume that the defect cores are localized, that is, there exists We define the minimum separation distance of D by We assume L D N throughout this paper and we refer to Figure 1 for a two dimensional example with A defined by (4.3) specifying a triangular lattice.

GENERALISATION ANALYSIS OF MLIPS 3
as a variational problem in a discrete energy space and establish qualitatively sharp far-field decay estimates for the corresponding equilibrium configuration. We will review the framework, and adapt it to the case of multiple point defects considered in this work. This will provide the analytic foundation of our generalisation analysis. For the sake of simplicity of presentation we will skip over some technical details but fill these gaps in Section 6.1. Let d 2 {2, 3} be the (e↵ective) dimension of the system. A homogeneous crystal reference configuration is given by the Bravais lattice ⇤ h = AZ d , for some non-singular matrix A 2 R d⇥d . We admit only single-species Bravais lattices. There are no conceptual obstacles to generalising our work to multi-lattices, however, the technical details become more involved. The reference configuration for a system with defects is a set ⇤ ⇢ R d . The mismatch between ⇤ and ⇤ h represents possible localized defected configurations.
For a wide range of atomistic materials simulation tasks, periodic boundary conditions are by far the most popular choice. In this paper, we consider multiple point defects in a finite domain with periodic boundary condition. To that end, let where ⇤ N is the periodic computational domain and ⇤ per N is the periodically repeated domain. We consider n D point defects in ⇤ N , e.g. vacancies or interstitials, located at`i 2 R d . We define the minimum separation distance of a configuration D by L D := inf |`i `j| |`i,`j 2 D, i 6 = j . We refer to Figure 1 for a two dimensional example with A defined by (4.5).  A displacement field of the infinite lattice ⇤ is a map u : ⇤ ! R d . For`2 ⇤ and ⇢ 2 ⇤ `, we define the finite di↵erence D ⇢ u(`) := u(`+ ⇢) u(`). For a subset R ⇢ ⇤ `, we define D R u(`) := (D ⇢ u(`)) ⇢2R , and we consider Du(`) := D ⇤ `u (`) to be a finite-di↵erence stencil with infinite range. For a stencil Du(`), we define the stencil norms where N (`) defined by (6.3) is the set containing nearest neighbours of site`. The site potential is a collection of mappings V : (R d ) Λ− → R, which represent the energy distributed to each atomic site. To simplify the notation we assume that V (0) = 0 for all , which is equivalent to considering a potential energy-difference. We state the assumptions on the regularity and locality of the site potentials in Section 6.1 and refer to [9, §2.3 and §4] for a detailed discussion of those assumptions. If Λ = Λ h , we denote the site potential by We define the space of periodic displacements to be and N sufficiently large, the periodic energy functional reads An equilibrium defect geometry is obtained by solvinḡ where A := {u ∈ U per N : x + u ∈ Adm 0 (Λ per N )} is the admissible set with x : Λ per N → R, x( ) := and Adm 0 (Λ per N ) defined in (6.2) is represents a constraint preventing the collision of atoms.
For the purpose of error analysis we will need a strong stability condition [9,17] as well as qualitative information about the equilibrium. To that end, we assume the following: (S) There exists a single strongly stable coreū core defined by (6.8) in the infinite lattice Λ. For the sake of simplicity of presentation, we postpone the rigorous formulation of (S)to Section 6.1, Equation (6.10).
The following result shows that, given n D point defects in Λ N , there exists a strongly stable equilibrium of (2.3) as long as (S) holds and the core positions satisfy a minimum separation criterion from each other. The detailed proof is given in Section 6.2.
Suppose that (S) holds. Then, for each n D ∈ N, there exists a constant L 0 > 0 such that for any core configuration D satisfying L D ≥ L 0 , there exists a strongly stable equilibrium of (2.3) and it can be written as

4)
where the defect core truncation operator Π R is defined by (A.1) with radius R = L D /3, and where ω ∈ A satisfies with a constant C independent of N, n D , L D . The condition L D ≥ L 0 entails that defect cores do not overlap too strongly. Theorem 2.1 not only gives the existence of the equilibrium of the multiple point defects in a periodic domain, but also establishes its structure: The equilibrium can be decomposed into two parts, a truncated defect core centered at each point defect and a remainder term. Our generalisation analysis in the next section heavily relies on this result.

Main results
3.1. Error estimates. The computational cost of solving (2.3) on a large simulation domain is prohibitive when an electronic structure model is taken into consideration. Surrogate models are therefore introduced, fitted to an electronic structure model, to give reasonable approximations.
If we denote the surrogate site potential by V (resp. V h ), then the corresponding energydifference functionals are given by The resulting variational problem for the equilibration reads 2) The surrogate model enables large-scale simulation but obtaining ab initio training data on a large simulation domain remains intractable. Instead, one normally fits the parameters in the surrogate model to ab initio simulations on very small training domains containing at most a few hundred atoms. These should include all possible local snapshots one expects to encounter in the simulation; in our case, regions of homogeneous crystal or single defects. In the following we will give a rigorous framework within which this intuition is made precise and all the resulting errors are quantified.
To that end, given . Let E L (u) and E L (u) be the energy functionals defined on Λ L , and the equilibrium of the corresponding variational problem with E L (u) is denoted asū L .
Next, we introduce the matching conditions between the reference and the approximated models, in terms of the observations (energies, forces and force constants) of the possible configurations nearū L . To be more precise, let δ > 0 such that B δ (ū L ) ⊂ A L representing atomic displacements nearū L , then the matching condition for the energy is defined by To measure the force error, we introduce where · (U per L ) * is the dual norm of U per L . Furthermore, for the stability analysis, we require the force constant error (3.5) In particular, the force constant error on the homogeneous lattice is given by The following result provides a rigorous a priori error estimate of the geometry error Dū − Dũ 2 N (Λ N ) and of the formation energy error |E(ū) − E(ũ)| in terms of the size of the training domains and in terms of the fit accuracy of the surrogate model. The proof is given in Section 6.3.
Theorem 3.1. Suppose that (S) is satisfied and thatū is a strongly stable equilibrium of (2.3). Then, for L, L D sufficiently large, satisfying L ≤ L D , and ε FC and ε FC hom sufficiently small, there exists an equilibriumũ of the surrogate model (3.2) such that

8)
where both constants C G and C E are independent of N, n D , L. Table 1. Limiting cases of error decay with respect to L and the matching conditions ε E , ε F and ε FC hom .
The error estimates in the foregoing theorem identify how the geometry error and the error in formation energy depend on data-oriented approximation parameters: model accuracy on the training domain and its size, L. If we construct the approximated site potential V such that the matching conditions (from (3.3) to (3.6)) are exactly zero, we obtain rates of convergence in terms of L. Conversely, if L is sufficiently large, the errors then depend only on the matching conditions ε F , ε E . These limiting cases are summarized in Table 1. We will see in Section 4.2 that these rates are indeed sharp.
There are several insights we can gain from our results. First, they clearly tease out an issue that is -to the best of our knowledge -never discussed in the MLIPs literature: the size of training domains significantly affects the quality of the fitted model. Intuitively this happens because one cannot in practise obtain "perfect" snapshots. Secondly, we see the importance of fitting force constants in order to significantly reduce the effect that size of the simulation domain has. Finally, our estimates provide a clear guidance on how energy, force and force-constant observations should be weighted in the least squares loss function, in particular suggesting the optimal balance ε E ≈ (ε F ) 2 suggesting to put signifantly higher emphasis on the energy fit.

Numerical results
We propose a concrete implementation of MLIPs inspired by our generalisation analysis of Theorem 3.1 to confirm the analytical error estimates on model problems.

Constructions of MLIPs.
4.1.1. Parameterisation. First, we need to choose a parameterisation of the surrogate potential V . Although a wide variety of choices is available nowadays, we have opted for the linear atomic cluster expansion (ACE) [1,27,39] which has the advantage of achieving close to state of the art accuracy despite being a linear model [27]. However, there is no reason to believe that the specific choice of MLIP is essential in our tests. Briefly, in the ACE model, the surrogate potentialṼ is written as where B are the ACE basis functions and {c B } B∈B B B are the parameters that we will estimate by minimizing a least squares loss. The basis functions B are invariant under rotations, reflections and permutations of an atomic environment. A more detailed review of the ACE model and in particular its approximation parameters is provided in the Appendix C.

4.1.2.
Training sets and loss. Following our generalisation analysis (Theorem 3.1) we require that the ACE model matches the reference model in the sense of making ε E , ε F , ε FC and ε FC hom small. These measures of fit accuracy are specified in terms of max-norms over an infinite set of displacements, which is clearly computationally not tractable. At this point we make several departures from our rigorous analysis.
We first introduce the training set, R: The complete neighbourhood B δ (ū L ) used in the analysis is replaced with a finite number of random samples taken from B δ (ū L ). We fix the perturbed parameter to be 0.01 and denote the number of the configurations in R as N train := #R. Analogously we also produce a test set. The number of configurations in training and test sets will be specified for each individual example.
Next, we consider the construction of a loss function inspired by our theory. We cannot optimize directly but we propose three ad hoc approximations of the matching conditions: • We replace the max-norm with an 2 -norm to obtain a linear least squares problem.
• The force error ε F is defined in a dual norm in (3.4), but for sufficiently small training domains it is almost indistinguishable from a standard 2 -norm, namely − ∇E • Finally, we drop the the force constant errors ε FC and ε FC hom entirely from the loss function. In fact, we have found that constructing the training set R as described above and only fitting forces and energies already results in a sufficiently good accuracy of ε FC . We will give a numerical verification of this statement in Section 4.2, Table 2.
Given the foregoing approximations, the training set R constructed above and the parameterisation defined by (4.1), we determine the parameters {c B } by minimising the following loss function where W E and W F are additional weights that might depend on the configurations and observations. According to (3.8) (or Table 1), we choose W E W F in practice such that the balance ε E ≈ (ε F ) 2 can be achieved. The details will be provided for different model problems in the next section.
The loss function (4.2) is quadratic in the parameters {c B } and can therefore be minimised by using a QR factorisation. In our implementation we use a rank-revealing QR (rr-QR) factorisation [7] which provides a mechanism analogous to Tychonov-regularisation [39]. The regularisation parameter for rr-QR factorisation is set to be 10 −6 (we tested various values and it gave the best performance for all numerical experiments) throughout this work.

Numerical results. We present numerical tests for two examples:
(1) Toy model: We consider various configurations of point defects (vacancies and interstitials) in a two dimensional triangular lattice. As the reference model we will use an embedded atom model (EAM) [15] instead of an electronic structure model. This highly simplified scenario, and the fact that elastic fields decay more slowly in two dimensions, allows us to more easily perform large-scale simulation in which we can most clearly observe the expected convergence results. (2) Tight-Binding model: We will also perform tests on configurations of multiple vacancies in three dimensional Silicon (Si). For these tests we employ the NRL tight binding model [14,28,34] as the reference model; see Appendix D for a brief review.
All numerical tests are implemented in open-source Julia packages ACE1.jl [18] (for the ACE model) and SKTB.jl [19] (for the NRL tight binding model).

4.2.1.
Two-dimensional toy model. To demonstrate the main theory (Theorem 3.1) most clearly, we first explore two-dimensional multiple point defects systems where the reference model is given by an EAM potential [15]. We consider a two dimensional triangular lattice with r 0 chosen such that the triangular lattice becomes the ground state of the reference EAM potential. Using the triangular lattice and an EAM potential as the reference model instead of an actual ab initio model means that we can more easily perform large-scale tests in a wider parameter regime in order to narrow down the best choices. In this example, the size of the simulation domain Λ N is chosen to be N = 60r 0 . The separation distance L D is identical to the size of training domain L for the sake of simplicity.
Multi-vacancies: We first consider the multi-vacancies case. The corresponding simulation domain and training domain for two separated vacancies (n D = 2) are shown in Figure 2. The illustration of other two cases (n D = 3, 4) considered in this example is provided in Figure 10 in the Appendix B.  .2) as W E = 100 and W F = 1 in order to balance the matching conditions ε E ≈ (ε F ) 2 (cf. Table 1). The fitting parameters, fitting accuracy in terms of the root mean square error (RMSE) on the testing sets and fitting time for various size of training domain L are given in Table 2 in the Appendix B. The accuracy of the force constant on the homogeneous lattice is also verified in Table 2. As we discussed in Section 4.1.2, the force constants are already fitted well (relative RMSE are about 5%) despite only fitting forces and energies.
We first test the convergence of the geometry error Dū − Dũ 2 N and the error in energy |E(ū) − E(ũ)| with respect to the RMSE by studying the case of two separated vacancies ( Figure 2). Figure 3 shows that, for different size of training domain L, the error curves of geometry error and error in energy decrease near linearly and quadratically respectively as RMSE decreases, which perfectly matches our theoretical predictions from Theorem 3.1. Figure 4 plots the convergence of the geometry error Dū − Dũ 2 N and the error in energy |E(ū) − E(ũ)| against the size of training domain L. We observe that the convergence rates of geometry error perfectly match our theoretical predictions from Theorem 3.1 for all multivacancies cases. The only exception is the convergence of error in energy, a faster convergence   rate is observed numerically. We speculate that this is due to the additional symmetry of the defective lattice for the cases we consider here, and that this leads to additional cancellation that our general analysis does not capture. Interstitials-vacancies: Next, we test the interstitials-vacancies case. Figure 5 plots the corresponding simulation domain and training domain for one interstitial and one vacancy (n D = 2). The case of one interstitial and three vacancies (n D = 4) illustrated in Figure 11 in the Appendix B will also be considered in this example. We take N vac train = N int train = 100, and N vac test = N int test = 25 for both cases. To balance the errors (ε E ≈ (ε F ) 2 ), the additional weights in (4.2) are chosen to be: W int E = 100 and W int F = 10 for interstitial while W vac E = 10 and W vac F = 1 for vacancy. Here we put more weights on the interstitial to express the fact that the it generates a larger distortion of the surrounding lattice thus making it more challenging to fit than the vacancy cores.
As presented in the previous case, we take n D = 2 (cf. Figure 5) as an example to study the convergence of geometry error and error in energy with respect to the RMSE on testing sets. Similar to the multi-vacancies case, Figure 6 shows that the errors are reduced as the accuracy of corresponding ACE models is improved.
The decay of geometry error Dū − Dũ 2 N and corresponding error in energy |E(ū) − E(ũ)| against the size of training domain L are shown in Figure 7. We observe that the convergence rates roughly match our theoretical predictions from Theorem 3.1 which are summarized in Table 1.

4.2.2.
Convergence results for the NRL-TB model. We now move to testing our schemes when the reference model is an electronic structure model. We choose NRL-TB [14] as the reference model, which is a successful tight-binding model for Si; see Appendix D for a short review. Our choice of Si as the material is due to the fact that it is a semi-conducting material for which we have also strong theoretical and numerical evidence for the localisation of its interatomic forces [8], which is an essential ingredient in our analysis. The simulation domain is constructed by 10 3 unit cells, which contains 8000 Si atoms with periodic boundary conditions in all three directions. We select two vacancy sites which are separated from each other in Silicon bulk crystal. The MLIPs are fitted by following the construction in Section 4.1, where the parameters in building the basis functions B for Si are taken from [1, Section 7.5]. The total number of the configurations in the training and testing sets and the additional weights in (4.2) are chosen to be the same as that in multi-vacancies presented in the previous section.
The convergence of geometry error and energy error against RMSE for NRL-TB Si model is shown in Figure 8, where the predicted convergence is again observed for this electronic structure model. Figure 9 plots the decay of geometry error Dū − Dũ 2 N and corresponding error in energy |E(ū) − E(ũ)| against the size of training domain L. We observe that the convergence rates again perfectly match our theoretical predictions from Theorem 3.1 and Table 1 for NRL-TB Si reference model.

Conclusion
We proposed a framework for a generalisation analysis in a multi-scale setting, and used it to investigate the error propagation from fitting MLIPs on a small training domain to making predictions on a large simulation domain. As a prototypical example, we apply the framework to the case of simulating multiple (weakly interacting) point defects in a crystalline solid. Our analysis identifies what observations one should acquire from the reference model to obtain accurate predictions in this case. Our theoretical results partially justify existing best practices in the MLIP literature, but also provide more fine-grained qualitative information about how prediction accuracy depends on the choice of training data. This approach also suggests a new perspective on how to approach the collection of training data and the design of loss functions.
Our presentation here is restricted to simple-species Bravais lattices and point defects. Generalisations do require additional technical difficulties to be overcome, but there appears to be no fundamental limitation to extend the method and the results to multi-lattices and a range of other defects in some form. To conclude, we briefly discuss some of these possibilities as well as limitations which deserve further mathematical analysis and algorithmic developments.
• More complex crystalline structures: As mentioned above, the extension to multi-lattices is conceptually straightforward and the necessary technical details should be addressed in depth. • Straight dislocations: The extension to straight dislocations appears straightforward applying the techniques of [13,17,22]. We expect that the error estimates for straight dislocations depend not only on the force error but also the matching condition on the linear elasticity, even the nonlinear elasticity (virials) due to the long-range elastic field. • Grain boundaries or curved dislocation lines: These more complex crystalline defects require a much more significant degree of extrapolation than point defects or straight dislocations. Both simulations and rigorous analysis appear to be both conceptually and technically much more challenging. While our overarching strategies should still apply, it is far less clear whether our methodologies in this paper can be applied directly. • Uncertainty estimation: It is common knowledge that MLIPs have a fundamental limitation in that they lack a physical model for the phenomenon being predicted and thus have unknown accuracy when extrapolating beyond their training set. The uncertainty quantification (UQ) capabilities would be included to address this problem. Analyzing the propagation of uncertainty in the training procedure to predicted properties could be understood from a Bayesian statistics perspective, where some recent works [2,29] should provide appropriate references.
6. Proofs 6.1. Preliminaries. In this section, we introduce the concepts needed in the proofs of the main results. We review the framework of [9,17] to formulate the equilibration of a single point defect (n D = 1, vacancy or interstitial) as a variational problem in a discrete energy space and then give the strong stability assumption (S), which assumes the existence of a single stable core in the infinite lattice Λ. The deformed configuration of the infinite lattice Λ is a map y : Λ → R d . We can decompose the configuration y into where the parameter m > 0 qualifies the accumulation of atoms. For site ∈ Λ, we define the nearest neighbours set N ( ) as which is applied in the definition of the energy norm · 2 N (cf. (2.1)). We introduce the discrete energy space for infinite lattice with the associated semi-norm Du 2 N . We also define the following subspace of compact displacements The associated class of admissible displacements is given by In this paper, we make the following assumptions on the regularity and locality of the site potentials, which has been justified for some basic quantum mechanic models [8,10,12,31]. We refer to [9, §2.3 and §4] for discussions of more general site potentials.
(RL) Regularity and locality: For all ∈ Λ, V Du( ) possesses partial derivatives up to n-th order with n ≥ 3. For j = 1, . . . , n, there exist constants C j and η j such that for all ∈ Λ and ρ ∈ (Λ − ) j .
Although we defined the site potentials on infinite stencils (R d ) Λ− , the setting also applies to finite systems or to finite range interactions. It is only necessary to assume in this case that the potential V (g g g) does not depend on the reference sites g g g ρ outside the interaction range. In particular, we will denote by V Ω the site potential of a finite system with the reference configuration lying in Λ ∩ Ω.
Following the results in [9, Theorem 2.1], the energy-difference functional for infinite lattice reads The corresponding variational problem for the equilibrium state is where "arg min" is understood as the set of local minimizers. An auxiliary energy-functional needed in the following analysis is the energy of the homogeneous (defect-free) lattice We are ready to give the assumption on the existence of a single strongly stable core in infinite lattice. Letū core be the corresponding local minimizer of (6.8). Then we give the rigorous formulation of (S) as follows (6.10) To map the displacements defined on Λ to Λ per N , we introduce the operator T per N : U 1,2 (Λ) → U per N (Λ N ). One possible construction is given in [17,Section 7.3]. Then we briefly establish two auxiliary stability results based on (S) that will be needed in our subsequent analysis. We include them here for the sake of completeness and their proofs can be found in [17, Section B.2 and Theorem 7.7]. Proposition 6.1 (Phonon stability). Suppose that (S) holds, then there exists a constant c hom satisfyingc hom ≥c Λ > 0 such that Proposition 6.2. Suppose that (S) holds and the energy-difference functional E is defined by (2.2). Then, for N sufficiently large, there exists a constantc > 0 such that To conclude this section, we give the decay estimates of the equilibrium state for single point defect [9, Theorem 3.2]: Ifū core ∈ A (Λ) is a strongly stable solution to (6.8) satisfying (S), then there exists C > 0 such that where the truncation operator Π R is defined by (A.1) with radius R = L D /3.
The following Lemma provides an estimate on the residual of such an approximated solution in terms of L D .

Lemma 6.3.
Suppose z is the approximated solution to the variational problem (2.3) as defined in (6.12) with truncation radius R = L D /3. Then, there exists a constant L 0 > 0 such that, for L D > L 0 , where the constant C is independent of N, n D , L D .
Proof. Let r := R + 1 = L D /3 + 1. For any v ∈ U per N , we define v i := Π i r v for i = 1, . . . , n D , and v 0 : where the defect core truncation operator Π i r is defined by (A.1). For each i = 1, . . . , n D , Π i r is extended periodically with respect to Λ N since v ∈ U per N . Lemma A.2 implies that We then decompose the residual into three parts where the operator T per N maps the displacements from U 1,2 (Λ) to U per N (Λ N ). For the term T 1 , we obtain that where r = L D /3+1 and the last inequality follows from the locality of site potentials (RL) with a constant κ > 0.
To estimate T 2 , by using Lemma A.2, it is straightforward to obtain that where M 1 is the uniform Lipschitz constant of δE since E is (n − 1)-times continuously differentiable with respect to the D · 2 N norm [9,13]. The term T 3 can be similarly estimated by using [17, Lemma 7.6] Combing (6.15), (6.16), (6.17) and (6.18), by exploiting the assumption that L D N , for L D sufficiently large, we have which yields the stated result.
We then proceed to prove that δ 2 E(z) is positive, where z is given by (6.12). This result employs the ideas similar to those used in the proofs of [17,Theorem 7.7] and [22,Lemma 5.2], modified here to an periodic setting and extended to cover the case of multiple point defects.

Lemma 6.4.
Suppose z is the approximated solution to the variational problem (2.3) as defined in (6.12). Then, there exists a constant L 0 > 0 such that, for L D > L 0 , there exists c L D ≥c/2 so that Proof. We argue by contradiction. Suppose that there exists no L 0 satisfying (6.20), it follows that there exists a sequence of multiple point defects configurations D k := { k i } n D i=1 such that: (1) L k := L D k → ∞ as k → ∞; (2) for all k, let z k be denoted as the approximated solution defined by (6.12), we havec Hence, let v k ∈ U per N be a sequence of test functions such that Dv k 2 N = 1, we can obtain We now employ the result in [17,Lemma 7.9]. This states that there exists a sequence of radii, r k → ∞, for which we may also assume r k ≤ L k /3, so that for each i = 1, . . . , n D , where the construction follows from (6.14). It follows that where we have ensured that supp{w i k } for i = 1, . . . , n D only overlaps with supp{w 0 k } by choosing r k ≤ L k /3, and hence all other cross-terms vanish.
For term S 1 , we have where the last inequality follows from the Proposition 6.1. The term S 2 can be estimated similarly from Lemma A.3 For the cross-terms S 3 , by assuming r k ≤ L k /3, we deduce that Using the techniques in the proof of [22,Lemma 5.2], we can obtain that, for each i, Combining the estimates from (6.22) to (6.25), we have Furthermore, Lemma A.3 implies that for sufficiently large R we havec R ≥c/2 > 0. Hence, together with (6.21), we can obtain so for k sufficiently large, it contradicts the assumption thatc k <c/2 for all k.
We are ready to prove the Theorem 2.1.
Proof of Theorem 2.1. Applying Lemma A.1 with the results in Lemma 6.3 and Lemma 6.4, we can state that there exist L 0 > 0, where D satisfies L D ≥ L 0 , and z is an approximated solution defined by (6.12) corresponding to D. It follows that for any v ∈ U per N , there exists Writingū := z + ω yields the stated result.
6.3. Proof of the generalisation analysis. In this section, we give the detailed proof of the generalisation analysis (Theorem 3.1), which is the main result in this paper.
Proof. Applying the framework of the a priori error estimates in [11,30,32,33], we mainly divide the proof into several steps in order to apply the inverse function theorem (cf. Lemma A.1). 1. Stability: For any v ∈ U per N , we consider the stability of (6.26) From the results in Theorem 2.1, we can obtain thatū is strongly stable, namely there exists c > 0 such that Similarly as shown in (6.14), we split the test function v. For L ≤ L D , let r := L/3 + 1, we define v i := Π i r v for i = 1, . . . , n D , and v 0 : Lemma A.2 implies that Dv i 2 N (Λ N ) ≤ C Dv 2 N (Λ N ) for i = 1, . . . , n D . The term S 2 can be further split into three parts where we have ensured that supp{v i } for i = 1, . . . , n D only overlaps with supp{v 0 } by the choice of r and therefore all other cross-terms vanish.
For the term S 21 , for L sufficiently large and t ∈ [0, 1], we can Taylor expand the δ 2 E at the reference configuration (6.29) We first estimate S To estimate S 22 , recall the definition of the predictor (6.12) and the construction of v i , for each i = 1, . . . , n D and L sufficiently large, we have , v i , the term S 23 can be estimated similarly by Hence, combining the estimates from (6.26) to (6.33), for L sufficiently large and the matching conditions ε FC , ε FC hom sufficiently small, we have

Consistency:
We estimate the consistency error, for any v ∈ U per N , by where v i , i = 0, . . . , n D , are constructed by (6.27).
To estimate T 1 , for each i = 1, . . . , n D , we denote To estimate T 2 , we note thatū is smooth in supp{v 0 } for L sufficiently large according to Theorem 2.1. Hence, we can Taylor expand δ E and δE at the reference configuration where t ∈ [0, 1]. The term T 21 can be bounded by For T 22 , it can be estimated similarly Hence, combining from (6.35) to (6.40), we can obtain 3. Application of inverse function theorem: With the stability (6.34) and consistency (6.41), we can apply the inverse function theorem (Lemma A.1) to obtain, for L sufficiently large and the matching conditions defined in (3.3) and (3.6) sufficiently small, the existence of a solutioñ u to (3.2), and the estimate This completes the proof of (3.7).
4: Error in energy: Next, we estimate the error in the energy. Recall the definition of the predictor z by (6.12), we first spilt the error in energy into two parts The term E 1 can be bounded by where M 1 is the uniform Lipschitz constant of δE.
To estimate E 2 , by applying the technique used in (6.43), similarly we can obtain where M 1 is the uniform Lipschitz constant of δ E. Combining (6.42), (6.43) and (6.44), we obtain Lemma A.1. Let X, Y be Hilbert spaces, w ∈ X, F ∈ C 2 (B X R (w); Y ) with Lipschitz continuous Hessian, suppose that there exist constants µ, r > 0 such that then there exists a locally uniquew ∈ B X R (w) such that δF (w) = 0, w −w X ≤ 2r µ and Following [17] we define a family of truncation operators {Π i R } i ∈D , which we will apply to the single strongly stable coreū core . Let η ∈ C 1 (R d ; [0, 1]) be a cut-off function satisfying η(x) = 1 for |x| ≤ 4/6 and η(x) = 0 for |x| ≥ 5/6.
We denote T Λ as the canonical triangulation of R d whose nodes are the reference sites Λ (cf. [17, Section 2.1]). Let Iu be the piecewise affine interpolant of u with respect to T Λ , and A R := B 5R/6 \ B 4R/6 be an annulus, then we can define the truncation operator by In addition, we define Π R := Π 0 R . In this paper, we consider multiple point defects in a finite domain Λ N with periodic boundary condition. Hence, for R N , we extend Π R periodically with respect to Λ N . We state the following result concerning the approximation property of the truncation operator Π R , which follows from results in [17].
, there exists C > 0 such that, for R sufficiently large, In particular, ifū core is the core corrector satisfying (S), we have where C is independent of R and D.
Next, we show that the assumption (S) implies that δ 2 E(Π Rū core ) is also positive for sufficiently large R, which follows directly from Proposition 6.2.
Lemma A.3. Supposeū core satisfies (S), then for sufficiently large R, there exist a constant c R > 0 such that andc R →c given in Proposition 6.2, as R → ∞.
Proof. For any v ∈ U per N , where the mapping T per N : U 1,2 (Λ) → U per N (Λ N ) is introduced in Section 6.1 and M 2 is the uniform Lipschitz constant of δ 2 E since E is (n − 1)-times continuously differentiable. As R → ∞, we can obtain the stated result.
Fitting accuracy on measurements and fitting time T  Table 2. 2D toy model: The fitting accuracy (RMSE), force constant error and fitting time for the case of two separated vacancies, where r 0 is the rescaled lattice constant.

Appendix B. Numerical supplements
We give the illustration of the simulation domain for n D = 3, 4 considered in this paper in Figure 10.   Figure 11 shows the illustration of the simulation domain for interstitial-vacancies case with n D = 4 considered in the numerical experiments. 1 N ! j 1 =··· =j N V N (g g g j 1 , · · · , g g g j N ), where the N -body potential V N : R dN → R can be approxi- φ k j j m j (g g g j ) with φ k m (r r r) := P k (r)Y m (r), r r r ∈ R d , r = |r r r|,r = r r r/r, where P k , k = 0, 1, 2, · · · are radial basis functions (for example, Jacobic polynomials), and Y m , = 0, 1, 2, · · · , m = − , · · · , are the complex spherical harmonics. The basis functions are further symmetrised to a permutation invariant form, where the coefficients U k mi are given in [1, Lemma 2 and Eq. (3.12)]. It was shown in [1] that the basis defined above is explicit but computational inefficient. The so-called "density trick" technique used in [3,16,36] can transform this basis into one that is computational efficient. The alternative basis is with the coefficients c B . The corresponding force of this potential is denoted by F ACE . The family of potentials are systematically improvable (see [1, §6.2]): by increasing the body-order, cutoff radius and polynomial degree they are in principle capable of representing an arbitrary many-body potential energy surface to within arbitrary accuracy.
The fitting parameters in the foregoing expressions are determined by fitting to some highsymmetry first-principle calculations: In the NRL method, a database of eigenvalues (band structures) and total energies were constructed for several crystal structures at several volumes. Then the parameters are chosen such that the eigenvalues and energies in the database are reproduced. For practical simulations, the parameters for different elements can be found in [34].