Interatomic-Potential-Free, Data-Driven Molecular Dynamics

We present a Data-Driven (DD) paradigm that enables molecular dynamics calculations to be performed directly from sampled force-field data such as obtained, e.g., from ab initio calculations, thereby eschewing the conventional step of modeling the data by empirical interatomic potentials entirely. The data required by the DD solvers consists of local atomic configurations and corresponding atomic forces and is, therefore, fundamental, i.e., it is not beholden to any particular model. The resulting DD solvers, including a fully explicit DD-Verlet algorithm, are provably convergent and exhibit robust convergence with respect to the data in selected test cases. We present an example of application to C60 buckminsterfullerenes that showcases the feasibility, range and scope of the DD molecular dynamics paradigm.


Introduction
Present-day computational chemistry methodology, including all-electron calculations, Density Functional Theory (DFT), and other ab initio paradigms, enables the generation of vast sets of parameterfree force-field data for complex molecular systems (cf., e. g., [1] and references therein). However, for large systems of atoms molecular dynamics based on empirical interatomic potentials remains the only method of choice. In that setting, the empirical potentials required in the calculations are, inevitably, a main source of empiricism and error. Indeed, at present there does not exist a rigorous theoretical means of generating sequences V h of approximate interatomic potentials that are guaranteed to converge to the underlying-and unknown-exact potential V ab initio in a manner that additionally ensures convergence of trajectories. Instead, empirical potentials based on ad hoc assumptions and parametrizations are contrived and fitted to data without global control of errors or convergence guarantees (cf., e. g., [2] for the state of the art of empirical interatomic potentials).
The availability of vast quantities of high-fidelity force-field data (cf., e. g., [3,4,5]) suggests a game-changing paradigm shift whereby the data themselves are the sole basis of molecular dynamics calculations and the conventional interatomic-potential modeling step is eschewed altogether. Evidently, such a strict Data-Driven (DD) paradigm, if feasible, would have the immediate beneficial consequences of eliminating the biases, loss of information, empiricism and error that inevitably afflict empirical interatomic potentials.
In this work we present a provably convergent strict DD paradigm for molecular dynamics based solely on force-field data and demonstrate the feasibility of the paradigm with the aid of selected examples. The resulting DD solvers return approximate trajectories that converge to those of the underlying-and unknown-exact potential when computed from force-field data sets of increasing fidelity. Such data sets can be generated ab initio, stored, reused, merged and adapted to specific application domains along the tenets of active learning [6].
The specific DD paradigm considered in this work builds on similar approaches for static problems [7,8] and dynamic problems [9]. Specifically, the trajectories of the system are thought to take place in a phase space of atomic configurations and force fields, cf. Section 2.2. The objective is then to determine trajectories that satisfy Newton's laws of motion exactly, viewed as a constraint on forces and accelerations, while remaining as close as possible to a given force-field data set in the sense of some suitable distance, cf. Section 3. We show that the resulting problem for the trajectories can be given the structure of an optimal control problem, cf. Section 3.1, or a game-theoretical problem, cf. Section 3.2. In this latter formulation, we also show that the DD problem can be expressed in terms of an effective, or learned, force field. However, the effective force field is only implied and need not be computed explicitly in calculations.
We note, cf. Section 3.4, that a natural distance between local clusters of atoms that is invariant under relabeling of the atoms is provided by the Wasserstein distance of optimal transport [10]. In addition, the force-field data is invariant under the action of the Euclidean group of translations and orthogonal transformation, which requires the evaluation of distances between entire orbits of 1 the Euclidean group. The calculation of distances between local atomic configurations is one of the main computational bottlenecks of the DD solver and, for large atomic clusters, requires the use of relaxation techniquesà la Kantorovich [10] and possibly interior-point regularizations such as max-ent [11]. In addition, the DD solver entails frequent searches in large force-field data sets. Here again, efficient search algorithms originally developed for Big Data applications are in existence (cf., e. g., [12] and references therein) and can be deployed as part of the DD solver.
The convergence of the DD solvers with respect to the data, including a fully-explicit DD-Verlet algorithm derived by time discretization, cf. Sections 3.5 and 3.6, can be verified mathematically under simple data sampling scenarios, cf. Section 3.3, or numerically based on selected test cases, cf. Section 4.1. As a proof of concept, we also present a simple application of the DD-Verlet solver to C 60 buckminsterfullerenes based on synthetic data sampled from the Stillinger-Weber potential [13], cf. Section 4.2. In all these cases, the ability of the DD solvers to compute qualitatively correct trajectories from relatively small data sets is remarkable.

Classical molecular dynamics
In this section, we define the class of problems under consideration and set forth notational conventions. We begin by considering classical dynamics without any assumptions regarding the structure and properties of the force field such as locality, Euclidean invariance or other symmetries. The consequences of such additional structures are elucidated in subsequent sections. We also concern ourselves with the reformulation of classical dynamics as a variational problem for trajectories in phase space, which sets forth a natural framework for the Data-Driven paradigm developed subsequently.
2.1. General interatomic potentials. We consider a system of N atoms adopting configurations described by coordinates r i = (r α i ) 3 α=1 , i = 1, . . . , N , in three-dimensional Euclidean space, which we identify with R 3 equipped with the standard Euclidean metric. For shorthand, we denote by the collection of all the position vectors of the atoms in the system and refer to R 3N as the configuration space.
The motion of the atoms obeys Newton's second law where t denotes time, m i is the mass of atom i, f i (r) is the force on atom i due to atomic interactions, f ext i (t) are applied forces, r 0 ∈ R 3N is the initial configuration, v 0 ∈ R 3N is the initial velocity field and a superimposed dot denotes time differentiation. The force field of the system is the function . We note that the configuration space, identified with R 3N , is to be regarded as an affine space of points, whereas the force-field space, while also identified with R 3N , is to be regarded as a vector space. This distinction is consequential when metrizing configuration space, where a distance between sets of points must be defined, cf. Section 3.4 1 .
The interaction forces are conservative if there exists an interatomic potential V (r) such that Examples of classical interatomic potentials commonly used in practice are presented, e. g., in [15]. Recent work aimed at formulating a framework for developing physics-based and machine learning interatomic potentials is exemplified by [16].

2.2.
Phase-space reformulation. An alternative set-oriented representation of the force field f (r) is to view it as a graph D in phase space Z = R 3N × R 3N , namely, the space of all pairs (r, f ) of system configurations and forces 2 . In this representation, the force field is regarded as a material-specific 3Ndimensional manifold, or graph, in 6N -dimensional phase space Z characterizing the entire range of possible atomic interactions. 1 More general formulations allow for the configuration space to be a smooth manifold with force fields taking values in the corresponding cotangent spaces [14], but that degree of generality is not required here. 2 It should be noted that the term phase space in classical dynamics is often applied to the space of positions and momenta, whereas in this work we use the term to signify the space of positions and forces.
We shall additionally denote by Z the linear space of all trial trajectories z(·) ≡ (r(·), f (·)) of the system over a given time interval [0, T ]. Thus, the elements of Z are curves z(·) in Z parameterized by time, not necessarily satisfying the equations of motion or compatible with the force field of the system.
For a specific material, the phase-space trajectories of the system must take values in the force-field graph D at all times, i. e., they must be contained in the set Thus, D is the set of phase-space trajectories compatible with the force field of the system and may thus be regarded as a material-specific trajectory set, or material set for short. Euclidean invariance requires that (r, f ) ∈ D if and only if (Q (r − c), Qf ) ∈ D for all translations c ∈ R 3 and orthogonal transformations Q ∈ O(3). Thus, D must contain entire E(3)-orbits and be invariant under the action of the Euclidean group E(3). In addition, we recall that the configuration space consists of point sets and should, therefore, be invariant under relabeling of the atoms, i. e., ( . . , N }. The physically admissible trajectories of the system are additionally subject to the constraint set forth by the equations of motion (2.1). The collection of all such admissible trajectories defines the admissible-trajectory set, or admissible set for short, We note that the admissible set E depends parametrically on the applied force field f ext (t) and the initial conditions (r 0 , v 0 ). However, the trial trajectories (r(·), f (·)) in E need not be solutions of the initial-value problem (2.1), since we do not require r(t) and f (t) to be related by the force field f (r) of the material. Evidently, the actual phase-space trajectories (r(·), f (·)) of the system, if they exist, lie in the intersection D ∩ E, i. e., are the admissible trajectories that are consistent with the force field of the material, or, equivalently, the material trajectories that are consistent with the equations of motion and initial conditions.

Data-Driven reformulation
Suppose that, as is often the case in practice, the graph D of the force field is not known in its entirety, but only through an approximating sequence of data sets D h , h = 0, 1, . . . . For instance, the sequence (D h ) may consist of increasing collections of points (r, f ) in phase space Z obtained by means of ancillary ab initio calculations or by some other means. As noted before, D h must contain entire E(3)-orbits in order to be invariant under the action of the Euclidean group E(3). In addition, the configurations sampled in D h are clusters of points and, therefore, their representation must be invariant under relabeling of the atoms. Example 3.1 (Molecular oxygen O 2 ). A simple notional example concerns molecular oxygen O 2 , also known as dioxygen or diatomic oxygen, vibrating in a bound configuration without rotation. The ground state of O 2 has a molecular weight of 31.9988, bond length of 1.21Å and a cohesive energy of −2.58 eV/atom [17]. The oxygen molecule is held together by a strong O=O double covalent bond, each oxygen atom sharing two of its outer-shell electrons with the other atom. A sequence of three force-field data sets, containing #D 1 = 10, #D 2 = 34 and #D 3 = 100 points, computed using the Vienna ab initio Simulation Package VASP is shown in Fig. 1. The calculations use PAW-PBE as the pseudo-potential and exchange-correlation functional [18,19]. An energy cutoff of 400 eV is employed for the plane wave basis, and reciprocal space is sampled using a Γ-centered Monkhorst Pack grid [20].
The data sets D h define approximate material trajectory sets of the form  Thus, if D h is a sample of points in phase space, the corresponding material trajectory set D h consists of piecewise trajectories that 'visit' the orbits in D h in turn. It is clear that, in general, the intersection between the admissible trajectory set E and the approximating material trajectory sets D h may be empty, in which case no classical approximating solution exists. One way to circumvent these difficulty is to relax the notion of 'approximating solutions' generated from the sequence of data sets D h . In particular, we would like such Data-Driven solutions to converge to the exact trajectory as the data sets D h sample the exact force-field graph D with increasing fidelity.

3.1.
Optimal-control Data-Driven reformulation. In order to derive approximations directly from a sequence of data sets D h , we relax the notion of solution in the spirit of [7,8,9,21]to mean a pair (y h (·), z h (·)) of trajectories in phase space such that y h (·) ∈ D h is closest to E, and, conversely, z h (·) ∈ E is closest to D h , in some metric to be defined. Thus, y h (·) is the material trajectory that is closest to being admissible and z h (·) is the admissible trajectory that is closest to being material. Formally, (3.2) y h (·), z h (·) ∈ argmin dist y(·), z(·) : y(·) ∈ D h , z(·) ∈ E , where dist is some suitable distance between phase-space trajectories, to be defined. We remark that problem (3.2) does not identify, or 'learn', an effective force-field from the data. Instead, it endeavors to minimize discrepancy between admissible solutions, satisfying the equations of motion (2.1), and the data, as measured by a suitably defined distance.
Evidently, if the material data set D h coincides with the entire force-field graph D the solutions of the Data-Driven problem coincide with the exact classical solutions of the initial-value problem (2.1), as required. Most importantly, the Data-Driven problem continues to make sense-and returns approximating solutions-in the case of partial point data, thus setting forth an approximation scheme for trajectories. We expect the approximate Data-Driven trajectories to converge to exact classical trajectories as the data samples the force field with increasing fidelity.
We note that problem (3.2) has the structure of an optimal control problem with differential constraints [22]. To exhibit this structure, suppose, for definiteness, that the distance between trajectories in Z is of the integral form for some local distance dist(y, z) over Z, not renamed. Begin by minimizing out the material trajectories y(·) in (3.2) to define a cost function is the distance from state z ∈ Z to the data set D h . Evidently, the cost function F h (z(·)) measures the discrepancy between a trial trajectory z(·) and the force-field data. With these definitions, the abstract Data-Driven problem (3.2) reduces to As announced, this problem has the structure of an optimal control problem constrained by ordinarydifferential equations if, with z(·) = (r(·), f (·)), we identify f (·) as the control and r(·) as the state variable (cf. [22], Section 3.5). The problem is then to minimize the cost F (r(·), f (·)) with respect to the control f (·), with r(·) taken as a solution of the initial-value problem (2.1) for any trial force-field trajectory f (·). General conditions for existence and convergence of solutions for such problems are given in [22], Section 3.5.
In practice, the constraint (3.6b) can be enforced by means of Lagrange multipliers w(t), subject to initial conditions (3.7) w(0) = 0,ẇ(0) = 0, resulting in the Lagrangian L y(·), z(·), w(·) = 1 2 dist 2 y(·), z(·) to be rendered stationary. Assuming, for definiteness, an Euclidean distance of the form with y = (s, g), z = (r, f ), the Euler-Lagrange equations corresponding to variations in z(·) and w(·) are where we write y h (·) = (s h (·), g h (·)) and z h (·) = (r h (·), f h (·)) for the resulting Data-Driven solutions. In (3.10), (s h (t), g h (t)) is to be chosen at all times t ∈ [0, T ] as the point in D h closest to (r h (t), f h (t)). We note that, by this choice, problem (3.10) is reduced to two coupled second-order ordinary differential equations in time in the unknowns r h (·) and w h (·). This duplicate structure is reminiscent of static Data-Driven problems in which the equilibrium constraint is enforced by means of Lagrange multipliers [7].
Suppose that D h = D, i. e., the data supply a full representation of the force field of the system. Then, in (3.10a) and (3.10b) we have r h (t) = s h (t) and f h (t) = g h (t), which, together with (3.7) give w h (t) = 0. In addition, (3.10c) reduces to (2.1a), as required. In general, for underlying force-field graphs of sufficient regularity we expect the Lagrange multipliers w h (·) to tend to zero and the Data-Driven trajectories to converge to exact classical trajectories when the density of sampling increases and the data sets D h approximate D with increasing fidelity.

3.2.
Game-theoretical Data-Driven reformulation. An alternative Data-Driven paradigm that does learn an effective force-field from the data consists of recasting the Data-Driven problem (3.2) as a game-theoretical problem. To this end, we regard the functional −F h (r(·), f (·)), eq. 3.4 with z(·) = (r(·), f (·)), as the payoff for the force player. Evidently, −F h (r(·), f (·)) which measures the agreement between the trajectory z(·) = (r(·), f (·)) and the data set. In this reinterpretation, the objective of the force player is to determine a strategy f (·) that maximizes its payoff, or, equivalently, minimizes F h (r(·), ·), for given r(·). The force player competes against a second position player, whose objective is to determine a strategy r(·) that satisfies the initial-value problem (2.1) for given f (·). The objective of the position player can be expressed variationally by means of the action functional to be minimized with respect to r(·) at fixed f (·), with initial conditions (2.1b) replaced by the boundary conditions with r 0 and r T given. The minimizing property of the action functional, which normally only attains stationarity in general Lagrangian mechanics, is remarkable and owes to the independence of the force field. The Data-Driven problem then becomes with defines a non-cooperative game between the force and position players [23].
The difference between the game (3.13) and the optimal control problem (3.2) is subtle but significant. Thus, in the optimal control problem the trial position histories r(·) are tied to the trial control histories f (·) through the initial-value problem (2.1). By virtue of this constraint, the cost functional F h (r(·), f (·)) has a double dependence on the trial control histories f (·), once through its direct dependence and twice implicitly through r(·). By contrast, in the game problem (3.13) the same functional F h (r(·), f (·)) is minimized with respect to f (·) at fixed r(·). Therefore, the Data-Driven approximations (r h (·), f h (·)) generated by the two procedures are different in general.
Specifically, suppose that we fix a trial trajectory r(·) in configuration space. From (3.4) and (3.13a) if follows that the optimal force strategy f (·) of the force player is given at every time by the effective force field which assigns to every system position r(t) the force f (t) such that (s(t), f (t)) ∈ D h and s(t) is nearest to r(t). The optimal position strategy r h (·) then follows by solving the initial-value problem formally identical to (2.1) with the exact, but unknown, force-field f (r) replaced by the 'learned' Data-Driven force field f (r; D h ).
It bears emphasis that the 'learned' Data-Driven force field f (r; D h ) is not the product of modeling, but rather a result of analysis. Piecewise constant approximations of the force field have been analyzed by Gonzalez et al. [24], who termed the approximation scheme force stepping and analyzed its convergence properties. As in the optimal control formulation, we expect the approximate trajectories (r h (·), f h (·)) obtained from (3.15) to converge to the exact solution as the data sets D h approximate the graph D with increasing fidelity. General conditions for existence and convergence of solutions of general game-theoretical problems are given in [23].

3.3.
Convergence with respect to the data. The convergence properties of the Data-Driven solutions (r h (·), f h (·)) to the exact trajectory (r(·), f (·)) as the material data set D h samples the force-field graph D with increasing fidelity can easily be verified for simple scenarios. The analysis is greatly simplified by assuming long-term existence of trajectories (r(·), f (·)) and then focusing on the convergence of approximations (r h (·), f h (·)) thereof, in keeping with the main focus of this work.
Suppose that the data sets D h are generated by sampling the exact force field f (r) at points s ∈ R 3N and that there is h ↓ 0 such that for every r ∈ R 3N there is (s, g) ∈ D h with |r−s| ≤ h . Then, the optimalcontrol Data-Driven solutions r h (·) converge to the exact solution r(·) strongly in W 1,∞ ((0, T ), R 3N ).
The assumption of Lipschitz continuity places restrictions on the force field. For instance, for a dimer such as molecular oxygen O 2 , Example 3.1, Lipschitz continuity and inversion symmetry require f (r) to be continuous through the origin r = 0 with f (0) = 0, i. e., the interaction force must vanish when the two particles merge. We also recall that, if f (r) is differentiable, then the Lipschitz constant L is given by the maximum value of |Df (r)|, where Df (r) is the matrix of partial derivatives of f (r), or Hessian. Convergence in W 1,∞ means, in particular, that trajectories r(·) in configuration space are Lipschitz continuous and the velocitiesṙ(·) exist almost everywhere in [0, T ] and are essentially bounded.
The convergence of the game-theoretical Data-Driven solutions can be verified likewise. Distances between point sets are naturally measured by the discrete Wasserstein distance [10], i. e., the cost of transportation of one point set to another. Thus, let r = {r i } N i=1 and s = {s i } N j=1 be two clusters of points in configuration space R 3N (e. g., red and yellow in Fig. 2) and let f = {f i } N i=1 and g = {g i } N j=1 be force fields attached to them, respectively. We write y = (s, g) and z = (r, f ) for the corresponding points in phase space. For simplicity, we consider clusters with the same number of points, though the optimal transport framework can be extended to clusters of different size [10]. The discrete 2-Wasserstein distance between the E(3)-orbits of the point sets ∪ where the minimum is sought over the group S N of all permutations σ of the index set {1, . . . , N }, translations c ∈ R 3 and orthogonal transformations Q ∈ O(3), Fig. 2. Fixing σ ∈ S N and taking variations of (3.18) with respect to (c, Q) gives the Euler-Lagrange equations which, inserted into (3.19b), gives upon a redefinition of U . Evidently, the orthogonal transformation Q results from the polar decomposition of a 3 × 3 matrix that is explicit in the data. Once, Q is determined, the translation c follows from (3.20), and the corresponding trial distance as The sought distance (3.18) then follows by minimizing over the trial distances from all permutations of the indices. Evidently, the exact evaluation of the distance (3.18) is of combinatorial complexity due to the need to examine all permutations of the atoms and it is suitable only for small clusters. An alternative form of the distance is the Kantorovich reformulation [10]  Explicit approximations for the weights can further be obtained by recourse to interior point methods [11]. For instance, a maximum-entropy (max-ent) regularization [26] consists of replacing (3.23) which can be readily solved for the approximate weights µ β ij . A further minimization with respect to (c, Q), exploiting the optimality of µ β ij , results in Euler-Lagrange equations identical to (3.19). The 2-Wasserstein distance dist(y, z) then follows from the regularized, or thermalized, distance dist β (y, z) in the limit of β → +∞.
We note that the game-theoretical Data-Driven problem requires distances in configuration space only. Such distances follow from the preceding expressions as a special case by formally setting κ = 0.
3.5. Time discretization. Suppose now that we wish to approximate trajectories at discrete points t 0 = 0, t 1 , . . . , t n , t n+1 = t n + τ, . . . , T , where, for simplicity, we consider a constant time step τ with integer T /τ . The discrete trajectories are then sequences (r · , f · ) ≡ (r n , f n ) T /τ n=0 of points in phase space. We denote by Z τ the space of such discrete phase-space trajectories.
NB: As in the time-continuous case, Section 2.2, in order to carefully differentiate between points and discrete trajectories, henceforth we shall denote by r · , f · and z · = (r · , f · ) trajectories in Z τ . In particular, r n , f n and z n = (r n , f n ) denote the values of trajectories r · , f · and z · = (r · , f · ) at time t n = nτ ∈ [0, T ], respectively.
In the time-discrete setting, we may approximate the equations of motion by means of a general difference formula. For definiteness, we consider three-point formuale of the form Likewise, we identify the set of discrete material trajectories as (3.28) D τ = y · = (r · , f · ) ∈ Z τ : (r n , f n ) ∈ D, n = 1, . . . , T /τ , and the set of discrete material-data trajectories as The time-discrete versions of the optimal-control and game-theoretical Data-Driven formulations defined in the foregoing follow now mutatis mutandi through the introduction of a sui metric in the space Z τ of discrete phase-space trajectories. For instance, we may replace (3.3) by (3.30) dist τ y · , z · = T /τ n=0 dist 2 (y n , z n ) τ 1/2 . Then, the time-discrete version of the optimal-control Data-Driven problem (3.2) is (3.31) y h,τ · , z h,τ · ∈ argmin dist τ y · , z · : y · ∈ D h,τ , z · ∈ E τ .
As in the time-continuous case, the constraint set forth by the discrete equations of motion (3.25) can be enforced by means of a discrete Lagrange multiplier w · satisfying homogeneous initial conditions. Stationarity then results in two coupled second-order recurrence relations for r · and w · .
Likewise, a time-discrete version of the game-theoretical Data-Driven problem (3.13) can be set forth by introducing the discrete cost function, cf. eq. (3.4), (3.32) F h,τ (z(·)) = inf The objective of the force player is to determine a strategy f · that maximizes its payoff, or, equivalently, minimizes F h,τ (r · , ·), for given r · . The force player competes against a second position player, whose objective is to determine a strategy r · that satisfies the discrete equations of motion (3.25) for fixed f · . The objective of the position player can be expressed variationally by means of the discrete action functional to be minimized with respect to r · at fixed f · with given r 0 and r T /τ . The treatment of the time-discrete problems is otherwise identical to that of the time-continuous problems and details are therefore omitted in the interest of brevity.
Assuming that the time discretization scheme is strongly convergent, i. e., time-discrete trajectories converge to the exact time-continuous trajectories strongly in W 1,∞ ((0, T ), R 3N ) as τ ↓ 0, the joint convergence of the time-discrete Data-Driven trajectories with respect to data and time discretization follows directly from Props. 3.2 and 3.3 by passing to diagonal sequences.

Algorithm 1 Data-driven Verlet algorithm
Require: Number of atoms N ; r n−1 , r n and f n ; force-field data set D h ; applied loads f ext n . Then: i) Compute r n+1 from (3.34). iii) Find (s n+1 , g n+1 ) in D h such that s n+1 is closest to r n+1 . iii) Set f n+1 = g n+1 . Return r n , r n+1 and f n+1 .
3.6. Data-driven Verlet algorithm. An explicit time-discretization scheme commonly used in practice is Verlet's algorithm which is a special case of (3.25) with α −1 = α 1 = 0. It is readily verified that, in this case, the optimal control and game-theoretical Data-Driven problems coincide. The resulting time-stepping scheme is summarized in Algorithm 1. It bears emphasis that the Data-Driven Verlet algorithm retains the explicit character of the classical Verlet algorithm. The determination of (r n+1 , f n+1 ) requires a search over the force-field data set D h , which is the main computational bottleneck of the algorithm. Fast algorithms for searching large data sets may be found in [12] and references therein.

3.7.
Short-ranged force fields. As already remarked, the computation of distances between E(3)orbits of point sets is computationally intensive for large systems. Conveniently, many force fields of practical interest are short-ranged, which allows distance calculations to be restricted to small clusters defined according to the range of interaction.
Suppose, for simplicity, that the system under consideration is conservative and the attendant global force field f (r) derives from a global potential V (r). Suppose, in addition, that V (r) is short-ranged. In order to exploit this property, we define the local interatomic potentials in terms of indefinite integrals defined modulo additive constants. Thus, V i (r) is the primitive function of f i (r) with respect to r i with the coordinates of all other atoms held constant. We note that, by the definition of the local interatomic potentials, we have i. e., the force at atom i follows as the derivative of the local interatomic potential V i with respect to r i , with the coordinates of all other atoms held constant. It bears emphasis that definition (3.35) applies to arbitrary interatomic potentials without loss of generality. In addition, from (3.36) it follows that the collection of local interatomic potentials jointly supply the exact global force field of the system and, therefore, it affords an equivalent representation of said global force field.
For short-ranged interatomic potentials, we may expect the number of atoms involved in the definition of each local interatomic potentials to be much smaller than N . Two atoms i and j are non-interacting if (3.37) ∂f i ∂r j (r) = ∂ 2 V ∂r i ∂r j (r) = 0, and interacting otherwise. We denote by N i ⊂ {1, . . . , N } the subset of atoms that interact with r i , not including r i itself, and by #N i the number of atoms in N i . It then follows from (3.36) that f i (r) depends only on r i and its cohort N i of interacting atoms and that V i (r) can be chosen to have the same dependence up to inconsequential additive terms.
Example 3.4 (Monatomic chain). Consider a monatomic chain consisting of atoms interacting through the harmonic potential The local potentials satisfy the relations Integrating, we obtain up to inconsequential additive constants. We note that, whereas V (r) depends on the coordinates of all atoms, V i (r) depends only on the coordinate r i of the central atom and the coordinates of its nearest-neighbors N i = {i − 1, i + 1}.
We proceed to verify that, if the interatomic potential V (r) is invariant under the action of the Euclidean group E(3), then so are the local interatomic potentials V i (r). To this end, we note from definition (3.35) that, modulo inconsequential additive constants, for any curve γ(s) joining a reference configuration r 0 to r, i. e., such that r 0 = γ(0) and r = γ(1).
Since, for any c ∈ R 3 and Q ∈ O(3), the curve Q γ(s) + c joints the reference configuration Q r 0 + c to Q r + c, it follows that, modulo inconsequential additive constants, Suppose now that V (r) is invariant under rigid-body mappings. Then, modulo inconsequential additive constants, which establishes the invariance of V i (r). Finally, we note that the global initial-value problem (2.1) may be regarded as a game involving N players, one per atom in the system, with strategies r i (t) whose payoff is to satisfy their corresponding equations of motion (2.1a) and initial conditions. In variational form, each player seeks to render stationary its action functional with appropriate boundary conditions. The objective of the ith player in the game is then to determine a strategy r i (·) that renders A i (r i (·)) stationary. We note that, by (3.36), the Euler-Lagrange equations of the game so defined are identical to (2.1a), which shows the equivalence between the two paradigms. Evidently, the Data-Driven reformulation of the global molecular-dynamics problem set forth in the foregoing can be applied, mutatis mutandi, to each of the local problems. In order to formulate the local Data-Driven problems, it suffices to supply local force-field data sets D h,i collecting corresponding values of (r i , r |N i , f i ), i. e., local position r i of atom i, positions r |N i of the atoms in the local neighborhood N i of interaction and corresponding interaction force f i on atom i. The dimensionality of the local material data sets D h,i is thus much reduced with respect to that of the global data set D h and, in particular, is independent of the global size N of the system.

Numerical examples
Finally, we present examples of application aimed at demonstrating the feasibility, range and scope of the Data-Driven molecular dynamics paradigm.  Sequence of data sets sampled from a Morse potential [27] with constants fitted to experimental molecular spectra data [28,29], cf. Table 1 [30]. Because of the complexity of the quantum many-body problem, the accurate calculation of high-resolution vibrational spectra of even small molecules remains a challenge (cf., e. g., [31]). The calculations presented here are elementary and   Table 1. Molecular oxygen O 2 . Morse potential constants fitted to experimental molecular spectra data [29].
with constants fitted to experimental molecular spectra data [28,29], cf. Table 1. The molecule is assumed to undergo anharmonic vibration resulting, e. g., from a resonance induced by an external magnetic or optical field. In calculations, the molecule is initially at the equilibrium bond length r 0 = r e and the atoms are imparted an initial relative velocity v 0 = 50Å/ps. The time step is set to τ = 1 fs throughout. Synthetic force-field data sets are generated by sampling the Morse force field derived from (4.1) over the range of bond lengths from 0.5r e to 1.5r e . The resulting force-field data sets are shown in Fig. 3. In view of the simplicity of the Morse potential, the agreement with ab initio data, Fig. 1, is quite remarkable.
The corresponding DD-Verlet solutions are shown in Fig. 4. The 100-data point solution shows ostensible deviations from the Morse solution. Such deviations notwithstanding, it is remarkable that a qualitatively correct solution, showing the expected smoothness, anharmonicity and periodicity, is obtained from such a small data set. The 1000-data point and subsequent DD-Verlect solutions are ostensibly converged. In order to make the convergence analysis quantitative, we measure the discrepancy between the DD-Verlet and Morse solutions using the weighted L 2 -norm [32] (4.2)

The norm quantifies both amplitude and frequency errors in infinite wave-train solutions. The simple identities
exemplify those properties. Fig. 5 shows a log-log plot of the error incurred by the DD-Verlet solutions with respect to the Morse solutions as a function of the size of the force-field data set. As expected from Props. 3.2 and 3.3, the DD-Verlet solutions exhibit robust convergence towards the solution of the underlying Morse force field, whence the data is sampled. It bears emphasis that the underlying Morse force-field is assumed unknown for purposes of the DD-Verlet calculations. In addition, all DD-Verlet results stem directly from the sampled force-field data and no intermediate interatomic-potential modeling step is required at any stage of the calculations.

4.2.
Buckminster fullerene C 60 . We close with an example concerned with the radiation-induced fragmentation of C 60 Buckminster fullerene, or buckyball for short. The aim of the example is to demonstrate the ability of the DD-Verlet algorithm to navigate complex conformational transitions using relatively small ab initio data sets. Figure 6. a) Relaxed ground-state configuration of the buckminsterfullerene C 60 using the Stillinger-Weber (SW) [13] potential parametrized as in Table 2. b) Local environment of an atom in the SW ground-state configuration of C 60 comprising: i) 3 nearest neighbors (black), at an equilibrium distance r NN < r c ≡ SW cutoff radius, and ii) 6 connected next-to-nearest neighbors (white) at a distance r NNN > r c .
Buckminster fullerene C 60 , discovered in 1985 by Robert Curl, Harold Kroto, and Richard Smalley [33], is an allotrope of carbon consisting of 60 carbon atoms in a spatial arrangement of positions and bonds in the form of a dotriacontahedron (32- sided figure), Fig. 6a. The term fullerene more generally refers to compounds C n consisting of n carbon atoms, each of which is bonded to three other carbons through two single bonds and one double bond, forming a closed surface.
Interactions of intense ultrashort laser pulses with C 60 fullerene, and the ensuing competition between ionization and fragmentation, have attracted considerable attention [34,35,36] as a means of studying the mechanisms of molecular energy deposition. At low intensity irradiation, the vibrational modes of C 60 are excited, especially its fundamental A g , or breathing, mode [37,38]. Beyond a certain intensity, C 60 undergoes fragmentation and decomposes into C n fragments [39,36].
For purposes of demonstration, we carry out direct and DD-Verlet calculations based on the Stillinger-Weber (SW) potential [13] (4.4) where the two-body and three-body interaction energies are chosen of the form with (4.6) θ ijk = acos r ij · r ik |r ij | |r ik | representing the angle subtended by the vectors r ij and r ik . In calculations, we use the parametrization of [40], Table 2.
We note that the SW potential has finite range with cutoff radius r c = b σ = 2.686Å and an equilibrium bond length r e = 1.59169Å. From these values, an examination of the local interaction relations (3.37) for an atom i in the SW ground-state configuration of C 60 determines the local interaction neighborhood N i to be of the form shown in Fig. 6b. The local interaction neighborhood comprises: 3 nearest neighbors (black), at an equilibrium distance r NN < r c ≡ SW cutoff radius, and ii) 6 connected next-to-nearest neighbors (white) at a distance r NNN > r c , for a total #N i = 9. The actual calculation of the Wasserstein distance and the associated transformations is carried using the comparison algorithm of [41]. Fig. 7 shows the time evolution of the average radius of the C 60 molecule after irradiation as computed directly from the SW potential. Fig. 7a corresponds to an initial radial velocity v 0 =  Table 2. Stillinger-Weber (SW) potential parameters for carbon [40]. 0.00382Å/fs, below the fragmentation threshold, and Fig. 7b to an initial radial velocity v 0 = 0.0382Å/fs, above the fragmentation threshold. In the former case, the fullerene undergoes an anharmonic vibration in the breathing mode, whereas in the latter case the molecule initially dissociates into identical pentarings C 5 that scatter ballistically. A sequence of snapshots showing the scattering of the pentarings following fragmentation is shown inset in Fig. 7b. We wish to ascertain whether these behaviors can be reproduced by DD-Verlet calculations based on force-field data sets sampled from the SW potential. To this end, we generate synthetic force-field data sets that tabulate local interaction neighborhoods N i and corresponding forces on the central atom r i . Since the dimensionality is too large for local configurations to be sampled uniformly, we adopt a goal-oriented, or importance-sampling, sampling strategy, whereby the goal is to bias the sampling towards local configurations that are likely to arise under the specific conditions of the application of interest. Since the trajectories followed by the C 60 molecule after low-intensity irradiation are expected to consist predominantly of radial expansion/contraction, we generate local configurations by randomly expanding/contracting the ground-state configuration, relaxing it and adding random . DD-Verlet simulation of C 60 (r c = 2Å). Time evolution of the average radius after low-intensity irradiation. Comparison of trajectories obtained directly from the SW potential and from DD-Verlet for local force-field data sets of size 100 and 500.
noise to the positions of the atoms. Since the effective cutoff radius r c , if any, is presumed unknown a priori, we generate data sets over a range of assumed values of r c . Examples of local configuration sets generated by this procedure assuming a cutoff radius r c = 2Å are shown in Fig. 8. Fig. 9 shows results of low-intensity irradiation DD-Verlet calculations for data sets generated assuming r c = 2Å compared to the direct calculations. As may be seen from the figure, the DD-Verlet trajectories do closely match the results of the direct SW calculation even for the small data set of size 100.
Generating good force-field data sets for the high-intensity irradiation case a priori is challenging, as many different structural changes must be covered by the data in order to reproduce the fragmentation behavior of the fullerene. We may circumvent this difficulty by exploiting the variational structure of the DD-Verlet solver in order to adapt the force-field data sets on-the-fly, in the spirit of active learning [6]. Thus, we recall that DD-Verlet uses the Wasserstein distance dist(r, s) from a local configuration r to find similar local configurations s in the force-field data set. Therefore, situations in which the trajectory wanders into poorly sampled regions of phase space are immediately identified, since they result in large distances to the data. To correct for such sampling deficiencies, we perform an explicit Stillinger-Weber calculation whenever in the calculation we encounter a local configuration r such that no local configuration s in the force-field data set satisfies the condition dist(r, s) ≤ , for some prespecified tolerance , i. e., if the distance from r to the data set exceeds the tolerance.  Table 3. DD-Verlet simulation of C 60 . Response up to 0.6ps following high-intensity irradiation. Number of force-evaluations for different values of the distance tolerance .
As may be seen from Fig. 10a, for a neighborhood cutoff radius of r c = 2Å , DD-Verlet exhibits robust convergence with respect to the data-set tolerance (and thus the data-set size) but converges to a trajectory that is at variance with SW. Evidently, in hindsight this behavior is to be expected since the data is generated assuming a cutoff radius that underestimates the actual value r c = 2.686Å. This underestimation results in local force-field data sets with poor coverage of the actual SW force-field. By contrast, Fig. 10b collects similar curves obtained from force-field data sets generated assuming a cutoff radius r c = 4Å. As may be seem from the figure, in this case the DD-Verlet trajectories converge robustly to the correct SW trajectory as the sampling tolerance is decreased. Remarkably, convergence is achieved for modest force-field data set sizes, Table 3.

Concluding remarks
We have presented a paradigm that enables molecular dynamics calculations to be performed directly from sample data, e. g., obtained from ab initio calculations, thereby eschewing the conventional modeling of the data by empirical interatomic potentials entirely. The DD paradigm is lossless, i. e., it uses the available data in its entirety; unbiased in the sense of being free of modeling assumptions, ad hoc representational 'features', parametrization and regression choices, and other extraneous elements; and, by its direct and exclusive reliance on data, it is modeling-error free.
The data used by the DD paradigm specifically samples local configurations and corresponding atomic forces and is, therefore, fundamental, i. e., it is not beholden to any particular model. This strict reliance on fundamental data is potentially far-reaching as it facilitates data generation, management and sharing. Thus, fundamental data is by its very nature fungible, i. e., data sets of diverse provenance can be merged directly to obtain larger data sets. This property opens the way to the standing of publicably-editable, public-domain, force-field data sets that can be built collaboratively worldwide.
Considerations of sampling error, reliability and data provenance can easily be incorporated into the DD paradigm as quality-control measures [8].
The computation of distances between clusters of atoms and the need to execute fast searches in large data sets are among the main computational bottlenecks of the DD paradigm. For large local atomic clusters, it becomes impractical to compute the discrete Wasserstein distance by examining exhaustively all the rearrangements of the atoms, which is an operation of combinatorial complexity. When the sampling is sufficiently dense, the ideas in [41] can be used to limit the combinatorial complexity of the distance calculation to cubic order. In other cases, a Kantorovich relaxation, which reduces the operation to the solution of a linear programming problem, interior-point regularizations thereof, and generally considerations of efficiency in the computation of the distance become of the essence. In addition, fast search algorithms for purposes of Big Data analysis have been extensively investigated and enable the efficient utilization of billion-size data sets (cf., e. g., [12] and references therein).
Evidently, the local force-field data sets used in the DD calculations need to provide adequate coverage of the regions of phase space traversed by the trajectories of interest. However, data sets covering the entirety of phase space uniformly are not possible and the sampling must of necessity be goal-oriented, i. e. suited to a particular class of trajectories. For instance, in the application to irradiated C 60 fullerenes presented above, the expectation that the molecules undergo predominantly radial expansion/contraction at low irradiation intensities has been used in order to fashion a priori goal-oriented data sets delivering fast convergence.
In more general settings, it is not always obvious at the outset whether a given data set is likely to supply adequate coverage. Conveniently, the DD paradigm is self-correcting in that regard. Thus, a trajectory that wanders into a poorly sampled region of phase space inevitably results in large distances to the data set being recorded during the calculations, which immediately flags the data as insufficient for the intended purpose. Furthermore, the poorly sampled regions of phase space thus identified can then be sampled by carrying out additional ab initio calculations, with the process repeated until the distance between the computed trajectory and the data set is below a prespecified tolerance at all times. The resulting iteration, something referred to to as active learning [6], sets forth an approximating scheme with respect to the data itself. If convergent, the scheme results in trajectories that are close to those of the underlying-and unknown-force field sampled by the data, a remarkable feat indeed.
These and other questions and improvements of the DD paradigm presented here suggest themselves as worthwhile directions for further research. By the assumption on the density of sampling, we have |r h (t) − s h (t)| ≤ h , whereupon (A.2) gives which shows that the Lagrange multiplier w h (·) tends to zero uniformly over [0, T ]. Let r(·) be the exact solution of the initial-value problem (2.1). Subtracting equations of motion and initial conditions, we obtain m(r h (t) −r(t)) + g h (t) − f (r(t)) = κ −2 w h (t). (A.4a) r(0) − r h (0) = 0,ṙ(0) −ṙ h (0) = 0. (A.4b) Integrating with respect to time and using the second of (A.4b), Estimating with the aid of (A.3), (A.6) m|ṙ h (t) −ṙ(t)| ≤ t 0 g h (t ) − f (r(t )) dt + T 3 mκ 2 h . Furthermore, by the Lipschitz continuity and sampling assumptions, Inserting into (A.6), By Pointcaré's inequality and the first of (A.4b), there is a constant C > 0 such that Finally, by the integral form of Gronwall's inequality [25], Integrating with respect to time and using the second of (A.11b), Estimating, By the Lipschitz continuity and sampling assumptions, there is s h ∈ R 3N such that By Pointcaré's inequality and the first of (A.11b), there is a constant C > 0 such that Finally, by the integral form of Gronwall's inequality [25], whence we conclude that the velocity error |ṙ h (t) −ṙ(t)| and, together with the first of (A.11b), the position error |r h (t) − r(t)| converge to zero uniformly over [0, T ], as advertised.