Interatomic-Potential-Free, Data-Driven Molecular Dynamics

INTERATOMIC-POTENTIAL-FREE, DATA-DRIVEN MOLECULAR DYNAMICS J. BULIN1, J. HAMAEKERS1, M. P. ARIZA2AND M. ORTIZ3,4 Abstract. 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. 1. 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 Vhof approximate interatomic potentials that are guaranteed to converge to the underlying–and unknown–exact potential Vab 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 arXiv:2208.04937v1 [physics.chem-ph] 6 Aug 2022 2 J. BULIN, J. HAMAEKERS, M. P. ARIZA AND M. ORTIZ 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 `a 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 C60 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. 2. 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 Natoms adopting configurations described by coordinates ri= (rα i)3 α=1,i= 1, . . . , N, in three-dimensional Euclidean space, which we identify with R3equipped with the standard Euclidean metric. For shorthand, we denote by r≡ {ri}N i=1 ∈R3Nthe collection of all the position vectors of the atoms in the system and refer to R3Nas the configuration space. The motion of the atoms obeys Newton’s second law mi¨ri(t) + fi(r(t)) = fext i(t),(2.1a) r(0) = r0,˙r(0) = v0,(2.1b) where tdenotes time, miis the mass of atom i,fi(r) is the force on atom idue to atomic interactions, fext i(t) are applied forces, r0∈R3Nis the initial configuration, v0∈R3Nis the initial velocity field and a superimposed dot denotes time differentiation. The force field of the system is the function f(r)≡ {fi(r)}N i=1 from configuration space R3Nto force-field space R3N. We note that the configuration space, identified with R3N, is to be regarded as an affine space of points, whereas the force-field space, while also identified with R3N, 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.41. The interaction forces are conservative if there exists an interatomic potential V(r) such that (2.2) fi(r) = ∂V ∂ri (r). 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 Din phase space Z=R3N×R3N, namely, the space of all pairs (r, f) of system configurations and forces2. In this representation, the force field is regarded as a material-specific 3N- dimensional manifold, or graph, in 6N-dimensional phase space Zcharacterizing the entire range of possible atomic interactions. 1More 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. 2It 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. DATA-DRIVEN MOLECULAR DYNAMICS 3 We shall additionally denote by Zthe linear space of all trial trajectories z(·)≡(r(·), f(·)) of the system over a given time interval [0, T]. Thus, the elements of Zare curves z(·) in Zparameterized by time, not necessarily satisfying the equations of motion or compatible with the force field of the system. NB: In order to carefully differentiate between points and trajectories, henceforth we shall denote by r,fand z= (r, f)points in configuration space R3N, force space R3Nand phase space Z=R3N×R3N, respectively; and we shall denote by r(·),f(·)and z(·)=(r(·), f(·)) trajectories in the same spaces defined over a given time interval [0, T]. In particular, r(t),f(t)and z(t)=(r(t), f(t)) denote the values of trajectories r(·),f(·)and z(·)=(r(·), f(·)) at time t∈[0, T], respectively. For a specific material, the phase-space trajectories of the system must take values in the force-field graph Dat all times, i. e., they must be contained in the set (2.3) D=(r(·), f(·)) ∈ Z : (r(t), f(t)) ∈D, t ∈[0, T]. Thus, Dis 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)∈Dif and only if (Q(r−c), Qf)∈Dfor all translations c∈R3and orthogonal transformations Q∈O(3). Thus, Dmust 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., (ri, fi)N i=1 ∈Dif and only if (rσ(i), fσ(i))N i=1 ∈Dfor all permutations σof the index set {1, . . . , 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, E=z(·)≡(r(·), f(·)) ∈ Z : mi¨ri(t) + fi(t) = fext i(t), t ∈[0, T]; r(0) = r0,˙r(0) = v0. (2.4) We note that the admissible set Edepends parametrically on the applied force field fext(t) and the initial conditions (r0, v0). However, the trial trajectories (r(·), f(·)) in Eneed 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. 3. Data-Driven reformulation Suppose that, as is often the case in practice, the graph Dof the force field is not known in its entirety, but only through an approximating sequence of data sets Dh,h= 0,1, . . . . For instance, the sequence (Dh) may consist of increasing collections of points (r, f) in phase space Zobtained by means of ancillary ab initio calculations or by some other means. As noted before, Dhmust 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 Dhare clusters of points and, therefore, their representation must be invariant under relabeling of the atoms. Example 3.1 (Molecular oxygen O2).A simple notional example concerns molecular oxygen O2, also known as dioxygen or diatomic oxygen, vibrating in a bound configuration without rotation. The ground state of O2has a molecular weight of 31.9988, bond length of 1.21 ˚ A 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 #D1= 10, #D2= 34 and #D3= 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 Dhdefine approximate material trajectory sets of the form (3.1) Dh=(rh(·), fh(·)) ∈ Z : (rh(t), fh(t)) ∈Dh, t ∈[0, T]. 4 J. BULIN, J. HAMAEKERS, M. P. ARIZA AND M. ORTIZ (a) 0 . 6 0 . 8 1 . 0 1 . 2 1 . 4 1 . 6 1 . 8 0 100 200 300 400 500 600 B in d i n g f o r c e [ e V / A ] B o n d l e n g t h [ A ] (b) 0 . 6 0 . 8 1 . 0 1 . 2 1 . 4 1 . 6 1 . 8 0 100 200 300 400 500 600 B in d i n g f o r c e [ e V / A ] B o n d l e n g t h [ A ] (c) 0 . 6 0 . 8 1 . 0 1 . 2 1 . 4 1 . 6 1 . 8 0 100 200 300 400 500 600 B in d i n g f o r c e [ e V / A ] B o n d l e n g t h [ A ] Figure 1. Sequence of force-field data sets of molecular oxygen O2computed using the Vienna ab initio Simulation Package (VASP). a) #D1= 10 points; b) #D2= 34 points; c) #D3= 100 points. Thus, if Dhis a sample of points in phase space, the corresponding material trajectory set Dhconsists of piecewise trajectories that ’visit’ the orbits in Dhin turn. It is clear that, in general, the intersection between the admissible trajectory set Eand the approximating material trajectory sets Dhmay 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 Dh. In particular, we would like such Data-Driven solutions to converge to the exact trajectory as the data sets Dhsample the exact force-field graph Dwith increasing fidelity. 3.1. Optimal-control Data-Driven reformulation. In order to derive approximations directly from a sequence of data sets Dh, we relax the notion of solution in the spirit of [7, 8, 9, 21]to mean a pair (yh(·), zh(·)) of trajectories in phase space such that yh(·)∈ Dhis closest to E, and, conversely, zh(·)∈ E is closest to Dh, in some metric to be defined. Thus, yh(·) is the material trajectory that is closest to being admissible and zh(·) is the admissible trajectory that is closest to being material. Formally, (3.2) yh(·), zh(·)∈argminndisty(·), z(·):y(·)∈ Dh, z(·)∈ Eo, 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. DATA-DRIVEN MOLECULAR DYNAMICS 5 Evidently, if the material data set Dhcoincides with the entire force-field graph Dthe 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 Zis of the integral form (3.3) disty(·), z(·)=ZT 0 dist2(y(t), z(t)) dt1/2, 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 (3.4) Fh(z(·)) = inf y(·)∈Dh 1 2dist2y(·), z(·)=ZT 0 1 2dist2(z(t), Dh)dt, where (3.5) dist(z, Dh) = inf y∈Dh dist(y, z), is the distance from state z∈Zto the data set Dh. Evidently, the cost function Fh(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 Minimize: Fh(z(·)) in Z,(3.6a) subject to: Eqs. of motion (2.1).(3.6b) 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,˙w(0) = 0, resulting in the Lagrangian Ly(·), z(·), w(·)=1 2dist2y(·), z(·) +ZT 0m¨r(t) + f(t)−fext(t)w(t)dt, (3.8) to be rendered stationary. Assuming, for definiteness, an Euclidean distance of the form (3.9) dist2(y, z) = N X i=1 |ri−si|2+κ2|fi−gi|2, with y= (s, g), z= (r, f), the Euler-Lagrange equations corresponding to variations in z(·) and w(·) are rh(t)−sh(t) + m¨wh(t)=0,(3.10a) κ2(fh(t)−gh(t)) + wh(t)=0,(3.10b) m¨rh(t) + fh(t)−fext h(t)=0.(3.10c) where we write yh(·) = (sh(·), gh(·)) and zh(·) = (rh(·), fh(·)) for the resulting Data-Driven solutions. In (3.10), (sh(t), gh(t)) is to be chosen at all times t∈[0, T] as the point in Dhclosest to (rh(t), fh(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 rh(·) and wh(·). This duplicate structure is reminiscent of 6 J. BULIN, J. HAMAEKERS, M. P. ARIZA AND M. ORTIZ static Data-Driven problems in which the equilibrium constraint is enforced by means of Lagrange multipliers [7]. Suppose that Dh=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 rh(t) = sh(t) and fh(t) = gh(t), which, together with (3.7) give wh(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 wh(·) 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 Dhapproximate Dwith 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 −Fh(r(·), f(·)), eq. 3.4 with z(·) = (r(·), f(·)), as the payoff for the force player. Evidently, −Fh(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 Fh(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 (3.11) G(r(·), f(·)) = ZT 0 N X i=1 mi 2|˙ri(t)|2−fi(t)−fext i(t)·ri(t)dt, to be minimized with respect to r(·) at fixed f(·), with initial conditions (2.1b) replaced by the boundary conditions (3.12) r(0) = r0, r(T) = rT, with r0and rTgiven. 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 fh(·)∈argmin Fh(rh(·),·),(3.13a) rh(·)∈argmin G(·, fh(·)),(3.13b) 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 Fh(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 Fh(r(·), f(·)) is minimized with respect to f(·) at fixed r(·). Therefore, the Data-Driven approximations (rh(·), fh(·)) 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 (3.14) f(r(t); Dh) = f(t) : (s(t), f(t)) ∈Dh,dist(s(t), r(t)) →min!, which assigns to every system position r(t) the force f(t) such that (s(t), f(t)) ∈Dhand s(t) is nearest to r(t). The optimal position strategy rh(·) then follows by solving the initial-value problem mi¨ri(t) + fi(r(t); Dh) = fext i(t),(3.15a) r(0) = r0,˙r(0) = v0,(3.15b) formally identical to (2.1) with the exact, but unknown, force-field f(r) replaced by the ’learned’ Data-Driven force field f(r;Dh). It bears emphasis that the ’learned’ Data-Driven force field f(r;Dh) 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 (rh(·), fh(·)) obtained from (3.15) to converge to the exact solution as the data sets Dhapproximate DATA-DRIVEN MOLECULAR DYNAMICS 7 the graph Dwith 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 (rh(·), fh(·)) to the exact trajectory (r(·), f(·)) as the material data set Dhsamples the force-field graph Dwith 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 (rh(·), fh(·)) thereof, in keeping with the main focus of this work. We recall that a function f:Rm→Rnis Lipschitz continuous if there is a constant L > 0 such that (3.16) kf(r0)−f(r00)k ≤ Lkr0−r00k. for all r0,r00 ∈Rm. The space of Lipschitz-continuous functions of time over [0, T] with values in Rn is denoted W1,∞((0, T); Rn); it is a Banach space with norm (3.17) krk1,∞= max ess sup t∈(0,T) |r(t)|, T ess sup t∈(0,T) |˙r(t)|, where ess sup denotes the essential supremum (cf., e. g., [25]). Proposition 3.2 (Convergence of optimal-control Data-Driven problem).Suppose that f(r)is Lipschitz continuous and the exact initial-value problem (2.1) has solutions r(·)in W1,∞((0, T); R3N). Suppose that the data sets Dhare generated by sampling the exact force field f(r)at points s∈R3Nand that there is h↓0such that for every r∈R3Nthere is (s, g)∈Dhwith |r−s| ≤ h. Then, the optimalcontrol Data-Driven solutions rh(·)converge to the exact solution r(·)strongly in W1,∞((0, T),R3N). The assumption of Lipschitz continuity places restrictions on the force field. For instance, for a dimer such as molecular oxygen O2, 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 W1,∞means, in particular, that trajectories r(·) in configuration space are Lipschitz continuous and the velocities ˙r(·) exist almost everywhere in [0, T] and are essentially bounded. The convergence of the game-theoretical Data-Driven solutions can be verified likewise. Proposition 3.3 (Convergence of game-theory Data-Driven problem).Under the assumptions of Prop. 3.2, the game-theoretical Data-Driven solutions rh(·)converge to the exact solution r(·)strongly in W1,∞((0, T),R3N). The proofs of the preceding propositions are straightforward and are presented in Appendix A for completeness. Michael Ortiz EMMC18 2022 Data-Driven Molecular Dynamics transport map polar decomposition! (c, Q) Figure 2. Schematic representation of the Wasserstein distance between the E(n)- orbits of two clusters of points (yellow and red). 3.4. The distance between E(3)-orbits of point sets. Recall that the force-field graph Dand material data sets Dhthereof are invariant under the action of the Euclidean group E(3), therefore contain entire E(3)-orbits, and are additionally invariant under permutations of the atoms. However, in practice actual material data sets are likely to consist of representatives of the E(3)-orbits pegged 8 J. BULIN, J. HAMAEKERS, M. P. ARIZA AND M. ORTIZ to a particular numbering of the atoms. Details of the calculation of distances between E(3)-orbits of point sets from representatives thereof are presented in this section. 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={ri}N i=1 and s={si}N j=1 be two clusters of points in configuration space R3N(e. g., red and yellow in Fig. 2) and let f={fi}N i=1 and g={gi}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 ∪N i=1{yi}and ∪N i=1{zi} in R3×R3is (3.18) dist2(y, z) = min σ∈SN,(c,Q)∈E(3) N X i=1 |ri−Q(sσ(i)−c)|2+κ2|fi−Qgσ(i)|2, where the minimum is sought over the group SNof all permutations σof the index set {1, . . . , N}, translations c∈R3and orthogonal transformations Q∈O(3), Fig. 2. Fixing σ∈SNand taking variations of (3.18) with respect to (c, Q) gives the Euler-Lagrange equations N X i=1 Q−1ri−(sσ(i)−c)= 0,(3.19a) N X i=1 Q−1ri⊗(sσ(i)−c) + κ2fi⊗gσ(i)=U=UT.(3.19b) From (3.19a), we obtain (3.20) c=1 N N X i=1 Q−1ri−sσ(i), which, inserted into (3.19b), gives (3.21) N X i=1 Q−1ri⊗sσ(i)+ri⊗1 N N X j=1 sσ(j)+κ2fi⊗gσ(i)=U=UT, upon a redefinition of U. Evidently, the orthogonal transformation Qresults from the polar decomposition of a 3 ×3 matrix that is explicit in the data. Once, Qis determined, the translation cfollows from (3.20), and the corresponding trial distance as (3.22) dist2 σ(y, z) = N X i=1 |ri−Q(sσ(i)−c)|2+κ2|fi−Qgσ(i)|2. 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] dist2(r, s) = inf µ∈RN×N (c, Q)∈E(3) nN X i=1 N X j=1 µij|ri−Q(sj−c)|2+ κ2|fi−Qgj|2:µij ≥0, N X j=1 µij = 1o, (3.23) which reduces the computation to the solution linear programming problem for the weights µij. 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) DATA-DRIVEN MOLECULAR DYNAMICS 9 by dist2(r, s) = inf µ∈RN×N (c, Q)∈E(3) nN X i=1 N X j=1 µij|ri−Q(sj−c)|2+ κ2|fi−Qgj|2+1 βlog(µij): N X j=1 µij = 1o, (3.24) 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 t0= 0, t1, . . . , tn, tn+1 =tn+τ, . . . , T , where, for simplicity, we consider a constant time step τwith integer T/τ. The discrete trajectories are then sequences (r·, f·)≡(rn, fn)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, rn,fnand zn= (rn, fn)denote the values of trajectories r·,f·and z·= (r·, f·)at time tn=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 mi rn+1 i+rn−1 i−2rn i τ2+α−1fn−1 i−fext i(tn−1)+ α0fn i−fext i(tn)+α1fn+1 i−fext i(tn+1)= 0, (3.25) with (3.26) α−1+α0+α1= 1, and denote by Eτthe set of discrete trajectories satisfying the discrete equations of motion (3.25) and initial conditions, namely, (3.27) Eτ=z.= (r., f.)∈ Zτ: (3.25),(r0, r1) given. Likewise, we identify the set of discrete material trajectories as (3.28) Dτ=y·= (r·, f·)∈ Zτ: (rn, fn)∈D, n = 1, . . . , T/τ, and the set of discrete material-data trajectories as (3.29) Dh,τ =y·= (r·, f·)∈ Zτ: (rn, fn)∈Dh, n = 1, . . . , T/τ. 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/τ X n=0 dist2(yn, zn)τ1/2. Then, the time-discrete version of the optimal-control Data-Driven problem (3.2) is (3.31) yh,τ ·, zh,τ ·∈argminndistτy·, z·:y·∈ Dh,τ , z·∈ Eτo. 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·. 16 J. BULIN, J. HAMAEKERS, M. P. ARIZA AND M. ORTIZ A(eV) B b λ (eV) γ κ σ (˚ A) 5.3790 0.5082 1.8945 18.7079 1.2 -0.5 1.41800 Table 2. Stillinger-Weber (SW) potential parameters for carbon [40]. (a) 0 1 2 3 4 5 3 . 8 5 3 . 9 0 3 . 9 5 4 . 0 0 4 . 0 5 r a d i u s ( A ) t i m e ( p s ) (b) Figure 7. Buckminsterfullerene C60. Time evolution of the average radius predicted by the SW potential after: a) low-intensity irradiation; b) high-intensity irradiation. Inset: computed configurations at 0.5ps intervals (plotted to scale). 0.00382˚ A/fs, below the fragmentation threshold, and Fig. 7b to an initial radial velocity v0= 0.0382˚ A/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 C5that 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 Niand corresponding forces on the central atom ri. Since the dimensionality is too large for local configurations to be sampled uniformly, we adopt agoal-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 C60 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 DATA-DRIVEN MOLECULAR DYNAMICS 17 (a) (b) (c) Figure 8. Local configuration data sets for assumed cutoff radius rc= 2˚ A. Number of points: a) #Dloc = 100; b) #Dloc = 500; c) #Dloc = 1000. 0 . 0 0 . 1 0 . 2 0 . 3 0 . 4 0 . 5 0 . 6 0 . 7 0 . 8 3 . 8 5 3 . 9 0 3 . 9 5 4 . 0 0 4 . 0 5 r a d i u s ( A ) t i m e ( p s ) S t i l l i n g e r - W e b e r D D ( 1 0 0 p o i n t s ) D D ( 5 0 0 p o i n t s ) Figure 9. DD-Verlet simulation of C60 (rc= 2˚ A). 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 rc, if any, is presumed unknown a priori, we generate data sets over a range of assumed values of rc. Examples of local configuration sets generated by this procedure assuming a cutoff radius rc= 2˚ A are shown in Fig. 8. Fig. 9 shows results of low-intensity irradiation DD-Verlet calculations for data sets generated assuming rc= 2˚ A 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 rto find similar local configurations sin 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 rsuch that no local configuration sin the force-field data set satisfies the condition dist(r, s)≤, for some prespecified tolerance , i. e., if the distance from rto the data set exceeds the tolerance. 18 J. BULIN, J. HAMAEKERS, M. P. ARIZA AND M. ORTIZ (a) 0 . 0 0 . 1 0 . 2 0 . 3 0 . 4 0 . 5 0 . 6 0 3 6 9 1 2 1 5 1 8 2 1 r a d i u s ( A ) t i m e ( p s ) S t i l l in g e r - W e b e r D D ( ε = 0 . 0 1 A ) D D ( ε = 0 . 0 5 A ) D D ( ε = 0 . 1 A ) (b) 0 . 0 0 . 1 0 . 2 0 . 3 0 . 4 0 . 5 0 . 6 0 3 6 9 1 2 r a d i u s ( A ) t i m e ( p s ) S t i l l in g e r - W e b e r D D ( ε = 0 . 0 1 A ) D D ( ε = 0 . 0 5 A ) D D ( ε = 0 . 1 A ) Figure 10. DD-Verlet simulation of C60. Time evolution of the average radius after high-intensity irradiation. Comparison of trajectories obtained directly from the SW potential and from DD-Verlet with sampling tolerance = 0.1˚ A, 0.05˚ Aand 0.01˚ A. a) Local force-field data sampled using a cutoff radius rc= 2˚ A; b) Local force-field data sampled using a cutoff radius rc= 4˚ A. Stillinger-Weber DD (= 0.1˚ A) DD (= 0.05˚ A) DD (= 0.01˚ A) 6000 269 370 1460 Table 3. DD-Verlet simulation of C60. 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 rc= 2˚ A , 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 rc= 2.686˚ A. 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 rc= 4˚ A. 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. 5. 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. DATA-DRIVEN MOLECULAR DYNAMICS 19 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 C60 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. Acknowledgments This work has been funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) via project 211504053 - SFB 1060; project 441211072 - SPP 2256; and project 390685813 - GZ 2047/1 - HCM. Hausdorff Center for Mathematics (HCM). P. Ariza gratefully acknowledges financial support from the Consejer´ıa de Econom´ıa y Conocimiento of the Junta de Andaluc´ıa, Spain, under grant number P18-RT-1485; and from Ministerio de Ciencia, Innovaci´on y Universidades under grant number RTI2018-094325-B-I00. Appendix A. Proofs of Propositions 3.2 and 3.3 Proof. (of Prop. 3.2) For simplicity, we assume that all atoms have identical mass mi=m. We work directly with the Euler-Lagrange equations (3.10). Integrating (3.10a) twice using (3.7) we obtain (A.1) wh(t) = −Zt 0Zt0 0 1 m(rh(t00)−sh(t00)) dt00dt0. Estimating, (A.2) |wh(t)| ≤ Zt 0Zt0 0 1 m|rh(t00)−sh(t00)|dt00dt0. By the assumption on the density of sampling, we have |rh(t)−sh(t)| ≤ h, whereupon (A.2) gives (A.3) |wh(t)| ≤ T2 mh, 20 J. BULIN, J. HAMAEKERS, M. P. ARIZA AND M. ORTIZ which shows that the Lagrange multiplier wh(·) 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(¨rh(t)−¨r(t)) + gh(t)−f(r(t)) = κ−2wh(t).(A.4a) r(0) −rh(0) = 0,˙r(0) −˙rh(0) = 0.(A.4b) Integrating with respect to time and using the second of (A.4b), (A.5) m( ˙rh(t)−˙r(t)) = Zt 0κ−2wh(t0)−(gh(t0)−f(r(t0))dt0. Estimating with the aid of (A.3), (A.6) m|˙rh(t)−˙r(t)| ≤ Zt 0gh(t0)−f(r(t0))dt0+T3 mκ2h. Furthermore, by the Lipschitz continuity and sampling assumptions, gh(t)−f(r(t))=f(sh(t)) −f(r(t))≤Lsh(t)−r(t)≤ L|rh(t)−sh(t)|+|rh(t)−r(t)|≤Lh+L|rh(t)−r(t)|. (A.7) Inserting into (A.6), (A.8) m|˙rh(t)−˙r(t)| ≤ LZt 0 |rh(t0)−r(t0)|dt0+LT +T3 mκ2h. By Pointcar´e’s inequality and the first of (A.4b), there is a constant C > 0 such that (A.9) m|˙rh(t)−˙r(t)| ≤ CL Zt 0 |˙rh(t0)−˙r(t0)|dt0+LT +T3 mκ2h. Finally, by the integral form of Gronwall’s inequality [25], (A.10) |˙rh(t)−˙r(t)| ≤ h1 + CL mtexp CL mtiLT m+T3 m2κ2h, whence we conclude that the velocity error |˙rh(t)−˙r(t)|and, together with the first of (A.4b), the position error |rh(t)−r(t)|converge to zero uniformly over [0, T], as advertised.  Proof. (of Prop. 3.3) Again we work directly with the Euler-Lagrange equations (3.15). Subtracting equations of motion and initial conditions, we obtain m(¨rh(t)−¨r(t)) + f(rh(t), Dh)−f(r(t)) = 0.(A.11a) r(0) −rh(0) = 0,˙r(0) −˙rh(0) = 0.(A.11b) Integrating with respect to time and using the second of (A.11b), (A.12) m( ˙rh(t)−˙r(t)) = −Zt 0f(rh(t0), Dh)−f(r(t0)dt0. Estimating, (A.13) m|˙rh(t)−˙r(t)| ≤ Zt 0f(rh(t0), Dh)−f(r(t0))dt0. By the Lipschitz continuity and sampling assumptions, there is sh∈R3Nsuch that f(rh(t), Dh)−f(r(t))=f(sh(t)) −f(r(t))≤Lsh(t)−r(t)≤ L|rh(t)−sh(t)|+|rh(t)−r(t)|≤Lh+L|rh(t)−r(t)|. (A.14) Inserting into (A.13), (A.15) m|˙rh(t)−˙r(t)| ≤ LZt 0 |rh(t0)−r(t0)|dt0+LTh. By Pointcar´e’s inequality and the first of (A.11b), there is a constant C > 0 such that (A.16) m|˙rh(t)−˙r(t)| ≤ CL Zt 0 |˙rh(t0)−˙r(t0)|dt0+LTh. DATA-DRIVEN MOLECULAR DYNAMICS 21 Finally, by the integral form of Gronwall’s inequality [25], (A.17) |˙rh(t)−˙r(t)| ≤ h1 + CL mtexp CL mtiLT mh, whence we conclude that the velocity error |˙rh(t)−˙r(t)|and, together with the first of (A.11b), the position error |rh(t)−r(t)|converge to zero uniformly over [0, T], as advertised.  References [1] P. Xu, E. B. Guidez, C. Bertoni, and M. S. Gordon. Perspective: ab initio force field methods derived from quantum mechanics. Journal of Chemical Physics, 148:090901, 2018. [2] The Knowledgebase of Interatomic Models. https://openkim.org/. [3] The NoMaD Repository. http://nomad-repository.eu/cms/. [4] The Materials Project. https://materialsproject.org/. [5] The NIST Materials Genome Initiative. https://mgi.nist.gov/materials-data-repository/. [6] B. Settles. Active Learning, volume 18 of Synthesis Lectures on Artificial Intelligence and Machine Learning. Morgan & Claypool, 2012. [7] T. Kirchdoerfer and M. Ortiz. Data-Driven computational mechanics. Computer Methods in Applied Mechanics and Engineering, 304:81–101, 2016. [8] T. Kirchdoerfer and M. Ortiz. Data-Driven computing with noisy material data sets. Computer Methods in Applied Mechanics and Engineering, 326:622–641, 2017. [9] T. Kirchdoerfer and M. Ortiz. Data-Driven computing in dynamics. International Journal for Numerical Methods in Engineering, 113(11):1697–1710, 2018. [10] L. C. Evans. Partial Differential Equations and the Monge-Kantorovich Mass Transfer. Current Developments in Mathematics. International Press of Boston, Somerville, MA, 1997. [11] F. Shu-Cherng and J. H. S. Tsao. Entropy optimization: interior point methods, volume 2, pages 544–548. Springer, Boston, MA, 2001. [12] R. Eggersmann, L. Stainier, M. Ortiz, and S. Reese. Efficient data structures for model-free Data-Driven computational mechanics. Computer Methods in Applied Mechanics and Engineering, 382:113855, 08 2021. [13] F. H. Stillinger and T. A. Weber. Computer simulation of local order in condensed phases of silicon. Physical Review B, 31:5262–5271, 1985. [14] R. Abraham, J. E. Marsden, and T. S. Ratiu. Manifolds, Tensor Analysis, and Applications, volume 75 of Applied Mathematical Sciences. Springer–Verlag, New York, NY, second edition, 1988. [15] M. Finnis. Interatomic forces in condensed matter, volume 1 of Oxford Series on Materials Modelling. Oxford University Press, New York, NY, 2003. [16] Mingjian Wen, Yaser Afshar, Ryan S. Elliott, and Ellad B. Tadmor. Kliff: A framework to develop physics-based and machine learning interatomic potentials. Computer Physics Communications, 272:108218, 2022. [17] David R. Lide, editor. CRC Handbook of Chemistry and Physics: A Ready-Reference of Chemical and Physical Data. CRC Press LLC, Boca Raton, FL, 85th edition, 2005. [18] P. E. Bl¨ochl. Projector augmented-wave method. Physical Review B, 50:17953–17979, Dec 1994. [19] John P. Perdew, Kieron Burke, and Matthias Ernzerhof. Generalized gradient approximation made simple. Physical Review Letters, 77:3865–3868, Oct 1996. [20] Hendrik J. Monkhorst and James D. Pack. Special points for Brillouin-zone integrations. Phys. Rev. B, 13:5188–5192, Jun 1976. [21] Sergio Conti, Stefan M¨uller, and Michael Ortiz. Data-driven problems in elasticity. Archive for Rational Mechanics and Analysis, 229(1):79–123, 2018. [22] D. Bucur and G. Buttazzo. Variational methods in some shape optimization problems. Birkh¨auser, Boston, MA, 2005. [23] T. Roubicek. Relaxation in Optimization Theory and Variational Calculus, volume 4 of Nonlinear Analysis and Applications. De Gruyter, Berlin/Boston, second edition, 2020. [24] M. Gonzalez, B. Schmidt, and M. Ortiz. Force-stepping integrators in Lagrangian mechanics. International Journal for Numerical Methods in Engineering, 84(12):1407–1450, 2010. [25] L. C. Evans. Partial Differential Equations, volume 19 of Graduate Studies in Mathematics. American Mathematical Society, Providence, RI, 1998. [26] M. Arroyo and M. Ortiz. Local maximum-entropy approximation schemes: a seamless bridge between finite elements and meshfree methods. International Journal for Numerical Methods in Engineering, 65(13), 2006. [27] P. M. Morse. Diatomic molecules according to the wave mechanics. ii. Vibrational levels. Physical Review, 34(1):57– 64, 1929. [28] T. J. Zielinski. Exploring the Morse potential. Journal of Chemical Education, 75(9):1191, September 1998. [29] K. P. Huber and G. H. Herzberg. NIST Chemistry WebBook, chapter Constants of Diatomic Molecules. Number 69 in NIST Standard Reference Database. National Institute of Standards and Technology, Gaithersburg MD, 20899, May 2022. data prepared by J. W. Gallagher and R. D. Johnson, III. [30] G. Herzberg. Molecular spectra and molecular structure, volume 1. van Nostrand, 1950. [31] O. B. Gadzhiev, S. K. Ignatov, M. Y. Kulikov, A. M. Feigin, A. G. Razuvaev, P. G. Sennikov, and O. Schrems. Structure, energy, and vibrational frequencies of oxygen allotropes on (n≤6) in the covalently bound and van der Waals forms: ab initio study at the CCSD(T) level. Journal of Chemical Theory and Computation, 9:247–262, 2013. 22 J. BULIN, J. HAMAEKERS, M. P. ARIZA AND M. ORTIZ [32] Ortiz M., E. D. Sotelino, and B. Nour-Omid. Efficiency of group implicit concurrent algorithms for transient finite element analysis. International Journal for Numerical Methods in Engineering, 28(12):2761–2776, 1989. [33] H. W. Kroto, J. R. Heath, S. C. O’Brien, R. F. Curl, and R. E. Smalley. C60: Buckminsterfullerene. Nature, 318(6042):162–163, 1985. [34] I. V. Hertel, T. Laarmann, and C. P. Schulz. Ultrafast excitation, ionization, and fragmentation of C60. Advances In Atomic, Molecular, and Optical Physics, 50:219–286, 2005. [35] J. Kou, V. Zhakhovskii, S. Sakabe, K. Nishihara, S. Shimizu, S. Kawato, M. Hashida, K. Shimizu, S. Bulanov, Ya. Izawa, et al. Anisotropic coulomb explosion of C60 irradiated with a high-intensity femtosecond laser pulse. The Journal of Chemical Physics, 112(11):5012–5020, 2000. [36] Z. Z. Lin and X. Chen. Ultrafast dynamics and fragmentation of C60 in intense laser pulses. Physics Letters A, 377:797–800, 2013. [37] W. Kr¨atschmer, L. D. Lamb, K. Fostiropoulos, and D. R. Huffman. Solid C60: a new form of carbon. Nature, 347(6291):354–358, 1990. [38] R. Meilunas, R. P. H. Chang, S. Z. Liu, M. Jensen, and M. Kappes. Infrared and Raman spectra of C60 and C70 solid films at room temperature. Journal of Applied Physics, 70(9):5128–5130, 1991. [39] S. C. O’Brien, J. R. Heath, R. F. Curl, and R. E. Smalley. Photophysics of buckminsterfullerene and other carbon cluster ions. Journal of Chemical Physics, 88(1):220–230, 1988. [40] F. F. Abraham and I. P. Batra. Theoretical interpretation of atomic-force microscope images of graphite. Surface Science Letters, 209:L125–L129, 1989. [41] Johannes Bulin and Jan Hamaekers. Similarity of particle systems using an invariant root mean square deviation measure, 2021. 1Fraunhofer Institute for Algorithms and Scientific Computing, Schloss Birlinghoven, 53757 Sankt Augustin, Germany, 2Escuela T´ ecnica Superior de Ingenier´ ıa, Universidad de Sevilla, Sevilla 41092, Spain, 3Hausdorff Center for Mathematics, Universit¨ at Bonn, Endenicher Allee 60, 53115 Bonn, Germany, 4Division of Engineering and Applied Science, California Institute of Technology, 1200 E. California Blvd., Pasadena, CA 91125.