FlashMD: long-stride, universal prediction of molecular dynamics
Authors: Filippo Bigi
Sanggyu Chong
Agustinus Kristiadi
Michele Ceriotti
Venue: NeurIPS 2025 spotlight ()
Keywords: molecular dynamics
atomistic modeling
geometric machine learning
statistical mechanics
Primary area: machine_learning_for_sciences
TL;DR: A method for the prediction of molecular dynamics trajectories using long time steps

Abstract: Molecular dynamics (MD) provides insights into atomic-scale processes by integrating over time the equations that describe the motion of atoms under the action of interatomic forces. Machine learning models have substantially accelerated MD by providing inexpensive predictions of the forces, but they remain constrained to minuscule time integration steps, which are required by the fast time scale of atomic motion. In this work, we propose FlashMD, a method to predict the evolution of positions and momenta over strides that are between one and two orders of magnitude longer than typical MD time steps. We incorporate considerations on the mathematical and physical properties of Hamiltonian dynamics in the architecture, generalize the approach to allow the simulation of any thermodynamic ensemble, and carefully assess the possible failure modes of a direct MD approach. We validate FlashMD’s accuracy in reproducing equilibrium and time‐dependent properties, using both system‐specific and general-purpose models, extending the ability of MD simulation to reach the long time scales needed to model microscopic processes of high scientific and technological relevance.

FlashMD: long-stride, universal prediction of molecular dynamics
Filippo Bigi"
Sanggyu Chong"
Institute of Materials Ecole Polytechnique Fédérale de Lausanne Lausanne 1015, Switzerland filippo. maag@eepfl.cmm ch
Institute of Materials
Ecole Polytechnique Fédérale de Lausanne
Lausanne 1015, Switzerland
sanggyu. wwwnggepfl.....mmm
Agustinus Kristiadi
Michele Ceriotti
Department of Computer Science Western University Vector Institute London, ON N6A 3K7, Canada wwwiistt@ww.... ca.
Institute of Materials
Ecole Polytechnique Fédérale de Lausanne
Lausanne 1015, Switzerland
michele. cew.ollieeriottteep...... ch
Abstract
Molecular dynamics (MD) provides insights into atomic-scale processes by inte- grating over time the equations that describe the motion of atoms under the action interatomic forces. Machine learning models have substantially accelerated MD by providing inexpensive predictions of the forces, but they remain constrained to minuscule time integration steps, which are required by the fast time scale of atomic motion. In this work, propose FlashMD, method to predict the evolu- tion of positions and momenta over strides that are between one and two orders magnitude longer than typical MD time steps. We incorporate considerations on the mathematical and physical properties of Hamiltonian dynamics in the architecture, generalize the approach to allow the simulation of any thermodynamic ensemble, and carefully assess the possible failure modes of such a long-stride MD approach. We validate FlashMD's accuracy reproducing equilibrium and time-dependent properties, using both system-specifie and general purpose models, extending the ability of MD simulation to reach the long time scales needed to model microscopic processes of high scientific and technological relevance.
general-purpose
Introduction
Simulations of atomic-scale systems are at the core of computational physics, chemistry, biology, and materials science [1]. Molecular dynamics (MD) in particular is a powerful tool for investigating the behavior of microscopic systems, from proteins [2-5] chemical reactions [6-8] and materials [9- 11]. MD elucidates atomic-scale structure and mechanisms by numerically solving the equations of atomic motion, driven by forces that can be estimated by quantum ground state electronic- structure calculations [12]. This has allowed simulating the time-evolution of various systems at the atomic scale, and probing thermodynamic and kinetic properties that are often difficult measure experimentally.
Despite its utility, MD has long been hindered by the trade-off between computational cost and accuracy of employed methods. Machine learning interatomic potentials [13-19] (MLIPs) have reme- died this in part by allowing the cheap approximation of otherwise expensive quantum mechanical
These two authors contributed equally to this work.
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

MD conti..
FlashMD
Embed.
MLP
Pilt ar)
COM
PiC
GNN (PET)
pP())
MLP
MLP
30.(1 +
COM
-q Ar)
Hal Sie
ction
Optirall
(g)
Enarg race Thersstane
Figure 1: Schematic overview of FlashMD, Atoms of the system at time step T are taken as inputs, with atomic numbers Z, and momenta (T( embedded into the node features h,, and relative coordinates 4ıy (7) embedded into the edge features Eiy of a GNN for the system. The node outputs are used to predict the new configuration p; ÂT) and Aq, (T + År)
multi-head manner. Center-of-mas mass constraints also enforced. Uncertainty quantification can be enabled shown in the navy inset. Optional filters for energy conservation enforcement, thermodynamic ensemble control, and random rotation are provided, as discussed below. Conventional MD would require Ar explicit numerical integrations to reach the final configuration as opposed to pass of FlashMD.
energies and forces. Nonetheless, MLIPs exhibit their constraints of having obey the physical symmetries and requiring expensive gradient computation to obtain the forces for MD propagation. Furthermore, stable and theoretically meaningful MD simulations require mathematical integration of the equations motion with sufficiently small time steps (~ 1 fs), limiting the simulations to a regime far removed from experimentally relevant time scales.
Motivated by the latest developments in MLIPs involving symmetry breaking and direct force learn- ing [20-22], well by generative approaches that construct representative atomic configurations without following the physically motivated equations of motion [23], this work will focus on the direct prediction of MD trajectories (sce Fig. 1). This approach avoids both the explicit calculation of interatomic forces and the numerical integration of the equations of motion, allowing one to use much larger strides compared to traditional MD integrators with corresponding, dramatic extension of the time scales accessible via atomistic modeling.
Our novel contributions are summarized as follows: (i) We provide thorough theoretical analysis of the problem of directly learning MD trajectories, discussing potential pitfalls. (ii) We introduce techniques for higher accuracy and larger time steps, e.g. by enforcing exact conservation of energy at inference time, proving its importance in stabilizing trajectories and reproducing physically correct thermodynamic
behavior. (iii) We generalize the approach to MD in arbitrary, experimentalll-reeeewee ensembles. (iv) We present universal models for direct, large-stride MD, capable of predicting trajectories across wide range of chemical systems with diverse structure and composition.
Background and related work
We aim to predict the sequence of position g, (T) and momentum p; (T) for each atom at the discrete time step time step t At). Conventional MD evolves the dynamics using positions and momenta of all atoms at time step to predict the new positions and momenta at time 7 + 1. Our goal here is to be able to take longer strides AT, skipping all intermediate steps, thereby achieving corresponding reduction in computational cost and eventually allowing practitioners to access much longer time scales in MD.
MD trajectory integrated with a small time step At (so the actual simulation time at
2.1 Molecular dynamics
In its simplest form, MD is the numerical solution of Hamilton's 's equations
dg:
0H
dp:
0H
(1)
dt
0p,-
dt
0q,

---

for an atomistic system with N atoms at positions /m and momenta {p. The Hamilto- nian function H describing the dynamics, in the absence of external perturbations, takes the form H ((p,. I = Ei= p /2m; + V (lq, (,1.)))))) where m, are the atomic masses and (( )) is the potential energy of the system. In practice, Eq. 1 is discretized using time step At. Among the many algorithms that could be used for numerical integration, the velocity Verlet (VV) algorithm [24] has become the standard due to its simplicity and the fact it preserves exactly some of the key properties of the underlying continuous Hamiltonian dynamics (Sec. 2.2). A single VV step reads
1 av At, P P; 2 dqq
q: 7 q, +
PiAt,
Ps mi
P; - P; .qq At.
1V
(2)
If integrated with sufficiently small At, the VV algorithm approximately conserves the energy of the system, making it suitable to sample the NVE thermodynamic ensemble (where the number of particles N, volume V and total energy E are fixed).
The NVE ensemble rarely corresponds to realistic experimental conditions. For this reason, variants of MD have been developed to target other types of ensembles [25, 26], accelerate their statistical sampling [27], account for nuclear quantum effects [28], etc. Practically, NVE MD can be modified via the inclusion of thermostats (e.g.. [29, 30]), which allows sampling the constant-temperrauue (NVT) ensemble, as well as the addition of barostats to sample constant-presstin Specialized Monte-Carlo moves are also used to access ensembles with varying number of particles and constant chemical potential (μVT). These different variants of MD are discussed further in Appendix D, and are all built around the NVE integrator, whose accuracy is therefore of central importance to achieve sampling of configurations with the correct probabilities.
(Np') ensemble.
Molecular dynamics with machine learning interatomic potentials The integrator in (2) is simple and computationally inexpensive, even though requires a small time step. The bottleneck is typically the evaluation of the force F; = -V; at every step, which is traditionally done using affordable but inaccurate empirical potentials, or accurate but very demanding quantum mechanical calculations. Over the last two decades, most efforts of applying machine learning (ML. to accelerate MD have revolved around MI.IPs that approximate the potential energy surface V( ) from quantum mechanical calculations at much reduced cost. Though MLIPs were traditionally focused on describing a specific chemical system, the last few years have seen the development of "universal" MLIPs (3 1-440, which aim to provide good accuracy across the whole periodic table, in principle allowing users to simply deploy the model for the desired system without further, dedicated training.
2.2 Physical and mathematical considerations
The potential (((g,)=2 obeys two fundamental physical symmetries: (i) Sy-invariance: (( is invariant with respect to permutations of atom indices; (ii) B(3)-invariance )
Symmetries of the potential energy function
is invariant with respect to translations, rotations and inversions of the atomic structure. Although traditional MLIPs incorporate all these symmetries through symmetry-constrained architectures, some recent models do not directly enforce rotational (and inversion) symmetry and use instead data augmentation strategies to encourage the model to approximately capture it at training-time [18, 22]. This allows the models to avoid expensive cquivariant operations and make inference more computationally efficient, without significantly affecting the physical observables in MD [20]. It should also be noted that preserving the two remaining symmetries (those with respect to translations and permutations), which are much harder to augment, naturally leads to the choice of graph neural networks (GNNs) for learning interatomic potentials, explaining their widespread adoption in last-generation MLIPs
IN
Conservative and non-conservative forces
Using forces that are the derivatives of V in the propagation of Hamiltonian dynamics conserves the total energy H. Even though this is a requirement to sample the NVE ensemble, some recent models have implemented the direct prediction of F; as an arbitrary vector field, leading to conceessevvetve dynamics [22, 32, 41]. These "direct-force" models do not even apply data augmentation to encourage energy conservation, as doing so involves computing the Jacobian of the forces, which is computationally impractical. Nevertheless, high model accuracy and the use of hybrid integration strategies [21, 42] allow non-conservative forces to be used in MD without generating noticeable unphysical artifacts, and with the efficiency benefits given by skipping the differentiation step.
3

---

Symmetries in molecular dynamics The possibility of performing stable MD with direct force predictions, and therefore without an underlying potential energy surface, suggests that an explicit potential energy function might not be necessary to predict the system evolution over time. It is therefore natural to consider the underlying symmetries in MD, and understand whether they can be incorporated efficiently into a notential-free model, or encouraged at training-time.
The VV algorithm, as presented in Eq. displays two fundamental symmetries, both of which are properties of the underlying continuous solution of Hamilton's ' equations: (i) MD is symplectic. That is, if (p, Ji {q are the momenta and positions a time step T1 and {p {q.)- are those at a different time step 72, then
Op ddh Opi,a dq.a
1,
(3)
Ppt.o dq:a dqr.a Opi,o
for all i= N...,N (ii) MD is rive--rversiblle This means that, considering positions and momenta evolved from T3 to 72 T1, then evolving {-p} } {q()== for т2 time steps will result in the state {-p,)- (q)) Furthermore, (iii) MID is (approximately) energy-conserving While energy conservation is exact for the solution of Eq. I, its discretization in the form of Eq. 2 is only approximately energy-conserving due to the finite integration step At. Since smaller At values afford better energy conservation, the latter is an important metric of integrator quality and sampling accuracy, which is why conventional MD is limited to short time steps and therefore simulation times.
N and a 2. y, 2. This crucially implies conservation of volume in phase space [43].
MD as time series The sequence of configurations generated by a MD trajectory obeys a few additional mathematical properties. (i) MD is Markovian. It is trivial to see that Eq. 2 defines a Markovian process, i.c. that the present time step 7 contains all the necessary information predict any future time step + Ar. (ii) MD is deterministic. Even though it is possible to describe MD in probabilistic framework (see e.g. Appendix D), barring effects due to numerical error and parallel programming, the initial conditions determine unequivocally the trajectory, both for infinitesimal and finite time steps. (iii) MD chaotic. Very often, even moderately complex systems simulated with MD exhibit a positive Lyapunov exponent [44] (i.e., the speed at which trajectories diverge), meaning that the phase-space distance of two very close initial states increases exponentially fast arithmetic, this implies that the targets of the learning exercise will become excessively noisy (and therefore difficult or impossible to learn) for large strides Ar. Furthermore, the Lyapunov exponent is highly dependent on the physical system under consideration.
as the simulation time progresses. Since MD is executed in
2.3 ML modeling of molecular dynamics trajectories
few previous works have applied ML techniques to model MD trajectories or their target distribu- tions. In the following, we provide a non-exhaustive set of related works broadly categorized by their conceptual approach. We focus on methods that aim to avoid VV integration entirely, rather than on methods based on multiple time-stepping [42] that reduce the computational cost by combining the evaluation of cheap-but-ig -utinccutti ML models and more expensive physics [45, 46] or ML-based approximations of V along the trajectory [47]. We also discuss how the nature and the practical implementation of previous approaches compare with the fundamental properties of MD.
Thermodynande ensemble generators Generative ML techniques have been applied to directly sample the thermodynamically accessible system configurations [48-51], which is especially useful for biomolecular systems with slow conformational transitions. Such models allow cheap predic- tion of the thermodynamic properties that conventionally require long MD simulations, but they disregard the time-dependent behavior of the systems and hence are unsuitable for investigating the physicochemical phenomena that can only be explained via the system's 's dynamics. More recently, the implicit transfer operator (ITO) has been proposed to extend Boltzmann generators to also re- cover, with stochastic trajectory, a measure of long-time dynamical processes [52]. Another relevant approach is Timewarp [53], which employs a normalizing flow as a proposal distribution within Markov chain Monte Carlo scheme targeting the Boltzmann distribution, achieving effective time steps on the order of 10s-106 fs for molecular systems.
Time-series approaches Several works [54-57] have interpreted the MD trajectory as a time series and adopted recurrent neural network (RNN)-type architectures, particularly the long short-

---

term memory [58] (LSTM), for the learning task. These models take time series of past system configurations as imputs make predictions of the future trajectory. Some have taken a stochastic approach to predict probable distribution of system states at a future time, and this has been successfully demonstrated in both all-atomic [54] and grarse- [56] contexts. However, in light of the Markovian nature of MD, sequence models such as LSTTM use redundant information and are not necessary to learn MD trajectories. Furthermore, the deterministic nature of MD as presented in Eq. 2 makes it superfluous to use probabilistic models (VAEs [59], normalizing flows [60], diffusion models [61], etc.).
Direct MD propagators This is the class of methods that most closely resembles our approach. It distinet from the former two approaches in that no generative approaches or multiple time step information is used. In this case, the MI. models take only the positions and momenta of atoms of the system at time step 7 as inputs and deterministically predict their changes at + Âr, hence "directly propagating" the dynamics. In MDnet by Zheng et al. [62], the chemical system is described as graph, with the edge features incorporating both the relative positions and momenta between the atoms. The model then predicts the changes in the positions and momenta for a fixed large time step AL. Very recently, Thiemann et al. have demonstrated Traj [63], an autoregressive equivariant network for direct MD prediction. Their framework has been shown to achieve good accuracy in reproducing NTT trajectories for individual molecular or bulk systems at a specific thermodynamic state point. It is also worth mentioning GICnet [64] and its transferable, transformer- beabooeeerrr-eeddddd variant MDtrajNet-1 1 [65], which learns a function that takes as inputs the initial positions, velocity, and At to return the positions at time f + At, the Equivariant Graph Neural Operator [66] (EGNO) approach, which predicts the evolution of the system at multiple times using equivariant temporal convolution in Fourier space, and the Graph Network-based Simulators [67, 68] (GNS), which have been developed for arbitrary particle-based systems without the chemical context.
3 The FlashMD framework
In this section, we explain the design choices made for "FlashMD", our proposed approach for the direct learning and prediction of MD trajectories.
3.1 Learning molecular dynamics trajectories with graph neural networks
Given that MD shares with interatomic potentials the properties of E(3) ((3-i)-qivvriancee and the use of atomic geometries inputs, we propose that FlashMD should be built on top of similar GNN architectures to those that have successfully been used to model machine-learnin potentials (MLIPs). In this work, we choose the Point-Edge Transformer [18] (PET), although any GNN architecture could be used. Compared to the original architecture in Ref. [18], we make two physically motivated modifications adapt it to the prediction of trajectories: (i) Each node state is also initialized using the particle momentum P.. encoded via multi-layer perceptron, in addition to the chemical species of the atom under consideration. The initialization of the edge states remains unchanged, and it includes the interatomic vectors q, q,- (ii) The outputs (1 + Ar) and q, (T Ar) are node properties and are therefore predicted via two distinct multi-layer perceptron heads starting from the node representation of PET.
interatomic
It should be noted that the "raw" FlashMD predictions are chosen to be mass-scaled, i.c. p //m and (g:- 9,) / This ensures treatment of displacements and momenta on equal footing for atoms of different mass, although a data-driven approach is also possible (App. A). Further details on the architecture and the training procedure are available in Appendices and B, respectively.
3.2 Addressing the many pitfalls of direct molecular dynamies predictions
Despite the fact that we have identified and justified graph neural networks as a suitable model to predict MD trajectories, there still exist many problematic aspects of this exercise which, if ignored, could make the resulting models practically useless. We note that these, as well as many of the theoretical considerations above, have been almost entirely ignored in previous related works.
ditt---distribution predictions Robust epistemic uncertainty schemes, capable of predicting errors associated with limited data sampling, are generally highly recommended when sampling
5

---

configurations using MLIPs [69]. They become essential for models that predict MD trajectories directly, that are less physically grounded and more susceptible to exhibiting pathological and unphysical behavior when queried outside of the domain of their training data.
Chaoticity The chaoticity of MD limits the time scale that can be reached with deterministic predictions. It also introduces aleatoric component to the model error, which varies in intensity depending on the system (see Sec. 2.2) and should be accounted for when building an uncertainty quantification scheme. For more information, see App. M.
Timee-eevessibility Time-reversibilityy one consequence of the symmetries of MD, can easily be incorporated by data augmentation [62]. This follows trivially from the discussion in Sec. 2.2.
Conservation of energy Conservation of energy is another consequence of the fundamental symmetries of MD, namely the translational symmetry of time. Given the radical importance of translational symmetries in 3D space, which make GÑNs so effective for MLIPs, it is clear that conservation of energy should be considered centerpiece of MD trajectory modeling, it encodes the symmetry in the time dimension. Previous works have not carefully monitored conservation of energy in their predicted MD trajectories. In FlashMD,
implement two approaches improve energy conservation: (i) we utilize errors in energy conservation during training, in addition to the terms corresponding to the position and momenta (see App. B), (ii) enforce energy conservation at inference time by rescaling momenta after each FlashMD step (see App. C). The latter adjustment makes it possible run long trajectories targeting the NVE ensemble, these would otherwise be affected by large energy drift (see App. F).
Symplecticity Symplectic behavior is together with energy conservation a necessary and sufficient condition for correct thermodynamic sampling. Unfortunately, penalizing non- somppllletttitty in the loss function is impractical, as evaluating Eq. (3) involves the computation of the full 3N x 3N Jacobian similar to energy conservation in direct force prediction [21]. We will discuss some numerical results on the violation of (3) by FlashMD, but mainly focus on the empirical measures of the accuracy of dynamics and sampling, through comparison with conventional MD simulations.
Symmetry breaking Although equivariant GNNs include strict rotational symmetries in the model, many GNN architectures do not enforce rotational equivariance explicitly [18. 36, 38]. Given that directly learning MD trajectories a fundamentally more challenging problem than learning potential energy surface, symmetry breaking might affect models for MD more than MLIPs. Therefore, if using rotationally unrestricted GNNs, we recommend correcting for rotational (and inversion) symmetry breaking at inference time, using similar techniques as those proposed for MLIPs [20]. Given that the PET architccture [18] also does not enforce rotational equivariance, we rotational and inversion augmentation at training time, and optionally perform random rotation(s) of the system at each step of FlashMD simulations (see App. A).
3.3 Generalization to arbitrary Bermodyyamii ensembles
FlashMD is trained reproduce, with longer stride, NVE trajectories obtained with VV integrator. However, nearly all other MD variants can be discretized (and are often implemented) using a split- operator formalism where VV is one of the components of the algorithm for a single time-step. This construction, which further discussed in App. D, allows using FlashMD to accelerate the majority of MD variants and ensembles.
4 Results
To demonstrate the capabilities of FlashMD, we trained two types of models: water-specific models trained on dataset of MD trajectories for liquid water, and general- on MD trajectories of structures sampled from the MAD dataset (see Ref. [39]). All reference MD simulations were performed with the PET-MAD universal MLIP [39]. For both the water-specific and universal cases, we trained separate models targeting different time strides. Full training details are available in App. B. The following subsections will focus on the testing of such models predicting meaningful physical observables for the corresponding systems. Ablation studies are discussed in
universal models trained

---

App. F, and timings are provided in App. G. Test set accuracy benchmarks against MDNet [62] and TrajCast [63] are reported in App. H and App. respectively.
4.1 Liquid water
Liquid water is central to many physical, chemical, biological and environmental processes, in great part thanks to its microscopic hydrogen- beggennbonding network and resulting physical and chemical properties. The study of liquid water at the microscopic level with MD is therefore a very active area of research [70-73], and we use it here as an example of how FlashMD models can accurately predict physical observables at the molecular level. For consistency with the universal model, we train the water-specific model on trajectories obtained with PET-MAD, even though its reference electronic-structure method (PBEsol [74]) is known to grossly overestimate the melting point of water. For this reason, we perform simulations with target temperature of 450 K, above the melting temperature of this model. Results for the q- T-I44PP/ model [75] at room temperature are discussed in Appendix J.
We focus on the evaluation of static observables, i.c. the equilibrium, time-independent properties of physical system which can be estimated as averages over MD) trajectories. As easy-to-compute diagnostics, we investigate the mean kinetic energy and atomic radial distribution functions NVT simulations, as well as the equilibrium density predicted by Np7 simulations. The mean kinetic energies (expressed as effective temperatures) are shown in Table It can be seen that, while the models without energy conservation can show pathological deviations from the target temperatures, the energy conservation enforcement approach described in App. C always recovers the correct temperatures for the overall simulation. However, significant deviations can still be observed for the global stochastic velocity rescaling [30] (SVR) thermostat in the kinetic energies resolved by atom type. This is spurious effect (classically, each degree of freedom should have mean kinetic energy equal to 22)) which is also observed in non-conservative force models [21]. The link to direct force prediction (that can be seen as the At - 0 of trajectory prediction) is also suggested by the fact that lack of equipartition is observed also for short-time FlashMD models, that have excellent validation accuracy. As a consequence, those based on Langevin dynamics, similar to what was done in Bigi et al. [21]. With a judicious choice, one can achieve accurate sampling of equilibrium properties, without reducing substantially sampling efficiency. However, local thermostats disrupt dynamical properties, to an extent that depends on the strength of the thermostat coupling, and a thorough quantitative analysis of dynamical properties require disentangling the effect of the long-stride sampling and that of the thermostats needed obtain accurate equilibrium sampling (see App. J). For these reasons, we primarily focus our quantitative analysis on independeenn equilibrium properties, and discuss examples where FlashMD qualitatively captures time-dependent behavior.
needs to employ local thermostats, such as
The atomic radial distribution functions (for the oxygen and hydrogen atoms, respectively), are shown in the left and center panels of Fig. 2. Here, it can be seen that FlashMD is able to correctly reproduce the distribution functions from the reference MD simulations, for both water-specific and universal models. To demonstrate the behavior when simulating the constant-pressure NpT ensemble, we also compute densities at ambient pressure (Fig. 2, right panel). The water-specific FlashMD models are
Table 1: Difference between effective and target temperatures in NVT simulations of liquid water using FlashMD, using different models and thermostats, and comparing results with and without enforcement of energy conservation. Characteristic times of 100 fs and 10 fs were used for the Langevin and SVR thermostats, respectively. The "all", "O", "H" labels refer to the subset of atoms under consideration. Subscripts on the results refer to statistical sampling errors. All units are in Kelvin; numbers close zero are better.
Without energy conservation
With enengy conservation
Model
I Langevin ATo
SVR ATo
Langevin ATo
SVR
ATa
ATn
4Ta
ATa
ATa
ATh
ATa
ATo
A7a
Water, fs
-1,3gg9
-1.1m
-1.4ps;
.04am
13.8
-74p%
-0.3mm
-1.6
0.3mm
-033m
-54,6%)
2.2
Water,
1.4m9)
-1.4az:
288
0.1m3)
-21429
10.8
0.4m%
1.2
0.1m
20.4m
10.0LD
Water, 16 fs
-0.2ga.)
-1.1g
0.2mm
.0.1am -16.009
7.8
1.309)
1.2
14rn
0.4
-10.7a
5.4
Universal, fs
33.8
36.2,19)
32.8.
8.0mo
-57.5
40.4,.
0.2oa:
-144n
0.9m
-0.5mm
-58.628
28.2a.
Universal, 4 is
10.7pm
7.3
12.5am
9.2
79.2gmm
-25.3am
40.7m%
-1.9am
-0.1p
0.4g4:
-31.am
16.0pm
Universal. 16
-22.5mm
-20.9nm
-23
2m.
-4.0pm
7.6g.3:
-9.8
0.409)
2.9
0.8mm
0.1
8.6
-4.4a
7

---

able to reproduce densities similar (although not statistically equivalent) to the reference MD. The universal models show significant deviations from the reference calculation, although smaller than the typical discrepancy expected when varying the details of the electronic structure calculation.
4.2 Universal long-stride sampling
Having compared and validated both system-specific and universal FlashMD models on the liquid water system, we now proceed to demonstrate the accuracy of the universal FlashMID models for more complex and chemically-divy the relevance for different classes of chemical and materials science problems: (i) We estimate the distribution Ramachandran angles for alanine dipeptide, a system that exhibits the basic features of protein dynamics and that is often used to benchmark sampling methodologies in biomolecular simulations; (ii) We model the meie-ttemperature dynamics of the AI(110) surface, a deceptively simple system that exhibits spontaneous formation of defects and surface pre-melting [76]; (iii) We compute the temperature- dependent diffusion coefficient Li atoms in lithium thiophosphate (LiPS) material that being investigated an electrolyte for solid-state batteries [77, 78), and that allows us to demonstrate the accuracy in capturing time-dependent properties.
systems. We consider three archetypal examples that showcase
Solvated alanine dipeptide We extend upon the previous liquid water simulations, generating constant-pressure trajectories of single alanine dipetide molecule solution, following closely the setup described in Morrone et al. [79]. The energy landscape is probed in three different simulations: MD using the PET-MAD potential with 0.5 fs stride, and universal FlashMD with 8 fs and 16 fs strides. Note that for this system 0.5 fs is a limit above which MD simulations exhibit severe instabilities. The Ramachandran plots (Figure 3a) show the characteristic distribution of molecular conformers in terms of the backbone dihedral angles. This is model for the backbone flexibility of proteins and demonstrates the ability of FlashMD to recover major features of the Ramachandran plot (particularly the low-energy basins in the S Φ < O range) with strides to 32 times larger.
Metal surface As discussed in previous work (76), the premelting behavior in AI(110) is charac- terized by two soft vibrational modes of the surface atoms: the top layer atoms show softening in the [001], or x component, and the second layer atoms show softening in the [110], or component. Figure 3b shows that the universal Flash MD models correctly describe this trend, with the mean square displacement (MSD) larger for the corresponding surface atoms and their associated softer vibrational modes at 500 K. The dynamical consequences of premelting is presented the trajectory traces of 3b, which is generated from the FlashMD simulation with 64 fs strides at 600 K. The trajectory traces provide a visual representation of the anisotropic softening of vibration for first and second layer atoms, well as the dynamic adatom formation pathways involving cooperative migration of both the surface and second layer atoms and effective adatom diffusion via exchange with nearby surface atoms. This demonstrates the ability of FlashMD not only capture material specific trends, but also describe meaningful dynamical behavior, despite the much larger strides.
4-
Density [kg/m²]
1000
3
goo(n)
universalie
universal
2-
950
925
(6.25 8o
water
untvenal
900
3
4
5
1
2
5
16
[A]
[A]
Time step [fs]
Figure 2: Comparison of physical observables obtained from MD (black) and FlashMD (other colors). Left and Center: radial distribution functions for oxygen and hydrogen atoms, respectively. from simulations in the NVT ensemble using the Langevin thermostat. Right: densities from simulations the NpT ensemble.

---

())
-1
[[K]
MD
FlackMD
FashMD
380
109
MSD [A?]
10-a
Layer
1000/7 [K-2]
Figure 3: Results of case studies conducted for the universal FlashMD models. (a) Ramachandran plots of the main backbone dihedrals for a simulation of solvated alanine dipeptide at 450 K. (b) Mean square displacement (MSD) of the Al (110) surface atoms 500 at different layers from the surface (B indicates the limiting value for the bulk). The premelting and defect formation phenomena are also visualized as traces of atomic positions from FlashMD simulation at 600 K, run with Ar 64 fs. The ideal atomic positions are also shown for reference. (c) Li conductivities of LigL, at varying T, along with the initial system configuration overlaid with traces of the Li atom positions from the FlashMD trajectory at 700 K. In both (b) and (c), traces are obtained with a moving average time to remove thermal fluctuations and visualize more clearly the diffusive behavior.
Solid-state electrolyte At high temperatures, the phase of lithium thiophosphate undergoes a phase transition to a superionic state that exhibits much higher conductivity. We reproduce simulations analogous those in Ref. [80], computing the Nernst- cell as a function of temperature. Results Figure 3c show that the universal FlashMD model successfully describes the superionic transition of y LiyPS, and predicts for it to take place 675 K. within the established transition temperature range determined in previous simulations using PET-MAD [39]. Li conductivities reasonably matched with the reference MD simulations, albeit with systematic over- and under-estimations in the low and high T regimes, respectively.
conductivity of LiPS
5 Discussion
Machine learning has been quietly revolutionizing the atomistic modeling of matter, accelerating the most time-consuming parts of physics-based calculations while striving to retain as much as possible of the underlying physical symmetries and constraints. As datasets and models grow in scale, there increasing interest in more radical approaches that trade the physical grounding of established practices for computational efficiency. Our work demonstrates that there is enormous potential in constructing GNN models that predict directly the evolution of atomic coordinates and momenta, allowing MD simulations to propagate with long strides, each replacing tens of costly force evaluations with miniscule time steps. Contrary to the very few previous works in this direction, which were restricted to reproducing MD trajectories for specific system in prescribed thermodynamic conditions, we show that our FlashMD architecture allows one to obtain a universal direct MD model that be applicable to different thermodynamic conditions and ensembles, and to wildly diverse atomic structures and compositions.
This is not to say that circumventing Hamiltonian dynamics is without problems. We highlight several ways an architecture similar to FlashMD could fail, by breaking some of the fundamental symmetries and conservation laws that are obeyed (at least approximately) by conventional integrators. We show how these shortcomings can be mitigated, e.g., by performing energy rescaling at inference time, or by including thermostats to control systematic drifts in the constant-energy trajectories. While the design choices FlashMD deviate from Hamiltonian dynamics, an alternative framework that preserves symplecticity (and therefore Hamiltonian dynamics) in the direct learning of trajectories is

---

proposed elsewhere [81]. Many of the shortcomings of non-symplectic FlashMD are shared by non- conservative force models, which have also become fashionable as a tool accelerate MD. We argue that the transformative speed-up afforded by FlashMD makes direct MD trajectory prediction more promising approach, at least when performing exploratory studies that require simulating long time scales. As shown concretely by challenging examples that simulate with bemi-quantitative accuracy three archetypal systems for biochemistry, surface science and energy technologies, FlashMD can already be applied to realistic simulation problems, capturing the essential equilibrium and dynamical processes while accelerating sampling significantly all cases (App. G).
In considering potential directions for further development, one should keep in mind that, contrary to the case of ML interatomic potentials that have been studied in great detail and brought to scale over the past decade, there is very little existing research on direct MD prediction. We recognize the possibility of incorporating further constraints the model architecture, or refinements to the training details, that can better enforce the conservation laws obeyed by the fine-grained VV integrator. One could also investigate scaling up the FlashMD universal model to more parameters and a larger trajectory datasets, implement modified architecture that targets multiple time strides with single model. Given that training relies on short MD trajectories built from a universal MLIP, it is relatively affordable to increase the dataset size by at least one order of magnitude. Moving forward, we envisage a future in which every MLIP would come with its own FlashMD-styl long-stride MD companion models, increasing even further the time scales within reach of ML-driven atomic-scale simulations.
Software and data
The FlashMD models presented in this work were trained with the metatrain package [82], and they support inference multiple simulation engines (including ASE [83], i-PI [84], LAMMPS [85]) through metatomic [82].
Helper functions to download universal FlashMD models and to prepare simulations are distributed with the flashmd package available on PyPI. Further information and instructions can be found at https: /t//f/aaad....m/////////////...m org, including links to the training datasets and scripts to reproduce the reported results on HuggingFace and Materials Cloud [86]. An example of the use of FlashMD is also /tgg/exmplees//ee wwtps////aaaa////aa/ee/aa/ haaamddfllasmmd demo html.
available at https: / /ttp::://ttmmsssiieiieee stannsssiiiiii stw.oonissccccooooooo....o.o
Besides the universal models presented in this work, which were trained to reproduce molecular dynamics at the PBEsol level of theory (through the PET-MAD universal potential), also make available FlashMD models based on the more accurate PSCAN functional which we recommend for scientific use. A separate set of case study results for these universal FlashMD models presented in App. I..
Acknowledgments and Disclosure of Funding
The authors would like to thank Davide Tisi and Federico Grasselli for help in processing LiPS trajectories. FB and MC acknowledge support from the NCCR MARVEL, funded by the Swiss Na- tional Science Foundation (SNSF grant number 182892) and from the Swiss Platform for Advanced Scientific Computing (PASC). MC and SC acknowledge funding from the Swiss National Science Foundation (Project 200020_214879) 214879).
References
[1] Daan Frenkel and Berend Smit. Understanding Molecular Simulation. Academic Press, London, second edition, 2002.
[2] David E. Shaw, Kevin J. Bowers, Edmond Chow, Michael P. Eastwood, Douglas J. Ierardi, John L. Klepeis, Jeffrey S. Kuskin, Richard H. Larson, Kresten Lindorff-Larsen, Paul Mara- gakis, Mark A. Moraes, Ron O. Dror, Stefano Piana, Yibing Shan, Brian Towles, John K. Salmon, J. P. Grossman, Kenneth M. Mackenzie, Joseph A. Bank, Cliff Young, Martin M. Deneroff, and Brannon Batson. Millisecond-scall molecular dynamics simulations on Anton. In Proc. Conf. High Perform. Comput. Netw. Storgee Anal. - SC 09, page 1, New York, New York, USA, 2009, ACM Press. ISBN 1-081----588--4--88 doi: 10. 4544///115499.55499999
10

---

[3] Kresten Lindorff-Larsen, Stefano Piana, Ron O. Dror, and David E. Shaw. How fast-folding proteins fold. Science, 334(6055):517-522, 2011. doi: 10.1126/sciemee.1200 1208351
[4] Paul Robustelli, Stefano Piana, and David E. Shaw. Mechanism of coupled folding-upont binding of an intrinsically disordered protein. Journal of the American Chemical Society, 142
(25): 1091-11101, 2020, doi: 10. w..1021/jac..0c02220
[5] Yohan Lee, Sanggyu Chong, Chiwoo Lee, Jihan Kim, and Siyoung Q. Choi. Structural determinants of chirally selective transport of amino acids through the a-hemolysin protein nanopores free-standing planar lipid membranes. Nano Letters, (H(2):881667) 2024, doi: 10. (0021/acs.aanolet..00..
[6] Yun Kyung Shin, Hyunwook Kwak, Alex V. Vasenkov, Debasis Sengupta, and Adri C.T. van Duin. Development of a reaxff reactive force field for fe/cr/o/s and application to oxidation of butane over pyrite-covered cr2o3 catalyst. ACS Catalysis, 5(12):7226-7233) 2015. doi:
10. 1021/acccaaal
[7] Pictro Vidossich, Agustí Llcdós, and Gregori Ujaque. First-principles molccular dynamics studies of organometallic complexes and homogeneous catalytic processes. Accounts of Chemical Research, 49(6): 2016. doi: 10. wll.0///ac..cccouuet 6600054.
[8] Christoph K. Jung, Laura Braunwarth, and Timo Jacob. Grand canonical reaxff molecular dynamics simulations for catalytic reactions. Journal of Chemical Theory and Computation, (5(11))580-581 6,, 2019. doi: 10. 1021 htt////...cct..99.
[9] Norman J. Wagner, Brad Lee Holian, and Arthur F. Voter. Molecular- simulations of two-dimensional materials at high strain rates. Phys. Rev. A, 45:8457-8477 Jun 1992. doi: 10. 100/PhyyRev
[10] Thomas E. III Gartner and Arthi Jayaraman. Modeling and simulations of polymers: A roadmap. Macromolecules, 52(3):755-778) 2019. doi:
[11] Sanggyu Chong, Sven M. Rogge, and Jihan Kim. Tunable electrical conductivity of flexible metal-organic frameworks. Chemistry of Materials, 34(1):254 265, 2022. doi:
10. (002laacs.chemmatrrr1003333
[12] R Car and M Parrinello. Unified Aproach for Molecular Dynamics and Density-Functtinal theory. R. Car and M. Parrinello.pda .df. Phys. Rev. Lett., 55(22)24711-244444 1985.
[13] Jörg Behler and Michele Parrinello. Generalized Neural-Nerwork Representation of High- Dimensional Potential-Energe Surfaces. Phys. Rev. Lett. 98(14);146401, April 2007. ISSN
0031-9007. doi: 10.1 103/PhysRewLeen.98. 146401.
[14] Albert P. Bartók, Mike C. Payne, Risi Kondor, and Gábor Csányi. Gaussian Approximation Potentials: The Accuracy of Quantum Mechanics, without the Electrons. Phys. Rev. Lett., 104 (13):166403, April 2010. ISSN 0031-9007. doi: 10.1 Ro.110//yysReeLett.... 104.136403.
[15] Han Wang, Linfeng Zhang, Jiequn Han, and Weinan E. DecPMD-kit: A deep learning package for many-body potential energy representation and molecular dynamies. Computer Physics Communications, 128:178-11 184, July 2018. ISSN 00104655.
doi: 10. 1...01//.cpc.211..000
[16] Ilyes Batatia, David Peter Kovacs, Gregor N. C. Simm, Christoph Ortner, and Gabor Csanyi. MACE: Higher order equivariant message passing neural networks for fast and accurate force fields. In Alice H. Oh, Alekh Agarwal, Danielle Belgrave, and Kyunghyun Cho, editors, Adv. Neural Inf. Process. Syst., 2022.
[17] Albert Musaelian, Simon Batzner, Anders Johansson, Lixin Sun, Cameron J. Owen, Mordechai Kornbluth, and Boris Kozinsky. Learning local equivariant representations for large-scale atomistic dynamics. Nat Commun, 14(1):579, February 2023. ISSN 2041-1723. doi: 10.1034 1038/ $41467-023-363399y
[18] Sergey Pozdnyakov and Michele Ceriotti, Smooth, exact rotational symmetrization for deep learning on point clouds. In Adv. Neural Inf. Process. Syst., volume 36, pages 79469-79501. Curran Associates, Inc., 2023.

---

[19] Ryan Jacobs, Dane Morgan, Siamak Attarian, Jun Meng, Chen Shen, Zhenghao Wu, Clare Yijia Xie, Julia H. Yang, Nongnuch Artrith, Ben Blaiszik, Gerbrand Ceder, Kamal Choudhary, Gabor Csanyi, Ekin Dogus Cubuk, Bowen Deng, Ralf Drautz, Xiang Fu, Jonathan Godwin, Vasant Honavar, Olexandr Isayev, Anders Johansson, Boris Kozinsky, Stefano Martiniani, Shyue Ping Ong, Igor Poltavsky, KJ Schmidt, So Takamoto, Aidan P. Thompson, Julia Westermayr, and Brandon M. Wood. A practical guide to machine learning interatomic potentials status and future. Current Opinion in Solid State and Materials Science, 35:101214, 2025. doi: https:/ddo..org/10.. 10wp///0.0011..ooooooo cos.is..ssssmmmmm 2025. 101214.
[20] Marcel F Langer, Sergey N Pozdnyakov, and Michele Ceriotti. Probing the effects of broken symmetries in machine learning. Mach. Learn.: Sci. Technol., 5(4)44L011, December 2024. ISSN 2632-2153, doi: 10. 1...08/2233221133/48660
[21] Filippo Bigi, Marcel Langer, and Michele Ceriotti. The dark side of the forces: assessing non- conservative force models for atomistic machine learning. arXiv preprint arXiv:2412.11569, 2024.
[22] Benjamin Rhodes, Sander Vandenhaute, Vaidotas Šimkus, James Gin, Jonathan Godwin, Tim Duignan, and Mark Neumann. Orb-v3: atomistic simulation at scale. arXiv preprint arXiv:2504. 062311, 2025.
[23] Frank Noé, Simon Olsson, Jonas Köhler, and Hao Wu. Boltzmann generators: Sampling equilibrium states of many-body systems with deep learning. Science, 365(6457) September 2019. ISSN 0036-8075, 1095-9203. doi: 10. 1266/ccienceea 147.
[24] Loup Verlet. Computer "Experiments" on Classical Fluids. I. Thermodynamical Properties Lennard-Jones Molecules. Phys. Rev, 159(1):98-103, July 1967. ISSN 0031-899X. doi: 10. 103/PhysRev. 159.98.
[25] Hans C. Andersen. Molecular dynamics simulations at constant pressure and/or temperature. 1980. ISSN 00219606. doi: $dU/dVS.
J. Chem. Phys.,
[26] M. Parrinello and A. Rahman. Polymorphic transitions in single crystals: A new molecular dynamics method. Appl. Phys., 52(12);7182-7190, 1981. ISSN 00218979, doi: 10.1063/1. 328693.
[27] A Laio and M Parrinello. Escaping free-energy minima. Proc. Natl. Acad. Sci., 99(20):
12562-12566, 2002.
[28] M Parrinello and A Rahman. Study of an F center in molten KCI. J. Chem. Phys., 80:860,
1984.
[29] G Bussi and M Parrinello. Accurate sampling using Langevin dynamics. Phys. Rev. E, 75(5):
56707, 2007. doi: 10.1 10
[30] G Bussi, Donadio, and M Parrinello. Canonical sampling through velocity rescaling. Chem. Phys., 126(1): (66(1):4401,, 2007.
[31] Chi Chen, Weike Ye, Yunxing Zuo, Chen Zheng, and Shyue Ping Ong- Graph Networks Universal Machine Learning Framework for Molecules and Crystals. Chem. Mater, 31(9): 3564-3572, May 2019. ISSN 0897-4756. 1520-5002. doi:
[32] Johannes Klicpera, Florian Becker, and Stephan Günnemann. GemNet: Universal directional graph neural networks for molecules. arxiv:2
2rxiv:2066.08903 002.000008888 2021.
[33] Chi Chen and Shyue Ping Ong. A universal graph deep learning interatomic potential for the periodic table. Nature Computational Science, 2(11):718-728 Nov 2022. doi: 10.
[34] Bowen Deng, Peichen Zhong, KyuJung Jun, Janosh Riebesell, Kevin Han, Christopher J. Bartel, and Gerbrand Ceder. Chgnet as pretrained universal neural network potential for charge-informed atomistic modelling. Nature Machine Intelligence, 103,103,,, Sep 2023. ISSN 2522-5839. doi: 10.1038/42256-03-00777
12

---

[35] Ilyes Batatia, Philipp Benner, Yuan Chiang, Alin M. Elena, Dávid P. Kovács, Janosh Riebesell, Xavier R. Advincula, Mark Asta, Matthew Avaylon, William Baldwin, Fabian Berger, Noam Bernstein, Arghya Bhowmik, Samuel M. Blau, Vlad Cärare, James P. Darby, Sandip De, Flaviano Della Pia, Volker I.. Deringer, Rokas Elijošius, Zakariya El-Machachi, Fabio Falcioni, Edvin Fako, Andrea C. Ferrari, Annalena Genneith--chhrieeer, Janine George, Rhys E. A. Goodall, Clare P. Grey, Petr Grigorev, Shuang Han, Will Handley, Hendrik H. Heenen, Kersti Hermansson, Christian Holm, Jad Jaafar, Stephan Hofmann, Konstantin S. Jakob, Hyunwook Jung, Venkat Kapil, Aaron D. Kaplan, Nima Karimitari, James R. Kermode, Namu Kroupa, Jolla Kullgren, Matthew C. Kuner, Domantas Kuryla, Guoda Liepuoniute, Johannes T. Margraf, Bean-Bogdan Magdäu, Angelos Michaelides, J. Harry Moore, Aakash A. Naik, Samuel P. Niblett, Sam Walton Norwood, Niamh O'Neill, Christoph Ortner, Kristin A. Persson, Karsten Reuter, Andrew S. Rosen, Lars I.. Schaaf, Christoph Schran, Benjamin X. Shi, Eric Sivonxay, Tamás K. Stenczel, Viktor Svahn, Christopher Sutton, Thomas D. Swinburne, Jules Tilly, Cas van der Oord, Eszter Varga-Umbrisn William C. Witt, Fabian Zills, and Gábor Csányi A foundation model for atomistic materials chemistry. arXiv preprint arXiv:2401 000096, 2024.
Tejs Vegge, Martin Vondrák, Yangshuai Wang,
[36] Mark Neumann, James Gin, Benjamin Rhodes. Steven Bennett, Zhiyi Li, Hitarth Choubisa, Arthur Hussey, and Jonathan Godwin. Orb: A fast, scalable neural network potential. arXiv preprint arXiv:2410.22557 2024.
[37] Han Yang, Chenxi Hu, Yichi Zhou, Xixian Liu, Yu Shi, Jielan Li, Guanzhi Li, Zekun Chen, Shuizhou Chen, Claudio Zeni, Matthew Horton, Robert Pinsler, Andrew Fowler, Daniel Zügner, Tian Xie, Jake Smith, Lixin Sun, Qian Wang, Lingyu Kong, Chang Liu, Hongxia Hao, and Ziheng Lu. Mattersim: A deep learning atomistic model across elements, temperatures and pressures. arXiv preprint arXiv:2405.04967 2024.
[38] Yutack Park, Jaesun Kim, Seungwoo Hwang, and Seungwu Han. Scalable parallel algorithm for graph neural network interatomic potentials in molecular dynamics simulations. Journal of Chemical Theory and Computation, 20(1 018445474888,, 2024. doi: 10. Ac0000000..
[39] Arslan Mazitov, Filippo Bigi, Matthias Kellner, Paolo Pegolo, Davide Tisi, Guillaume Fraux, Sergey Pozdnyakov, Philip Loche, and Michele Ceriotti. Pet-mad, a universal interatomic potential for advanced materials modeling. arXiv preprint 1rXii:2503.141136 2025.
[40] Haochen Yu, Mattco Giantomassi, Giuliana Materzanini, Junjie Wang, and Gian-Marew Rignanese. Systematic assessment of various universal machine-learning interatomic potentials. arXiv preprint arXiv:2403.055: 2024.
[41] Yi-Lun Liao, Brandon Wood, Abhishek Das, and Tess Smidt. Equiformerv2: Im- proved equivariant transformer for scaling to higher-degree representations. arXiv preprint arXiv:2306. 12059, 2023.
[42] M. Tuckerman, B. J. Bemne, and G. J. Martyna. Reversible multiple time scale molecular dynamics. J. Chem. Phys., 97(3):1990, 1992. ISSN 00219606. doi: 10.10377/.463137.
[43] Mark Tuckerman. Statistical Mechanics and Molecular Simulations. Oxford University Press, 2008.
[44] Giancarlo Benettin, Luigi Galgani, Antonio Giorgilli, and Jean-Marie Strelcyn. Lyapunov characteristic exponents for smooth dynamical systems and for hamiltonian systems; method for computing all of them. part I: Theory. Meccanica, 15(I):9-20, Mar 1980. ISSN 1572-9648. doi: 10. 00.1077/FF02228236 URL https: http:///dooooooo/////// 10. 10010007/BBBR /F000/22222222.6
[45] Venkat Kapil, Jörg Behler, and Michele Ceriotti. High order path integrals made easy. J. Chem. Phys., 145(23):234103, December 2016. ISSN 00219606. doi: 11111//1199144333
[46] Kevin Rossi, Veronika Jurásková, Raphael Wischert, Laurent Garel, Clémence Corminboeuf, and Michele Ceriotti. Simulating Solvation and Acidity in Complex Mixtures with First- Principles Accuracy: The Case of CH; SO H and 2 O in Phenol. J. Chem. Theory Comput., 16(8):5133-5144 August 2020. ISSN 1549-9618, 1549-9626. doi: 10. 1/221/ac.jccc /ct.ss/cccc..000033..
13

---

[47] Lars L. Schaaf, Ilyes Batatia, Christoph Brunken, Thomas D. Barrett, and Jules Tilly. Boostmd: Accelerating molecular sampling by leveraging ml force field features from previous time-steps. arXiv preprint arXiv:2412. 18633, 2024.
[48] Frank Noé, Simon Olsson, Jonas Köhler, and Hao Wu. Boltzmann generators: Sampling equilibrium states of many-body systems with deep learning. Science, 365(6457):eaaw 147, 2019.
[49] Giacomo Janson, Gilberto Maldes-Garcia Lim Heo, and Michael Feig. Direct generation of protein conformational ensembles via machine learning. Nature Communications, 14(1):774, Feb 2023.
[50] Leon Klein and Frank Noé. arXiv:2406. 14426, 2024.
Transferable
boltzmann generators.
arXiv preprint
[51] Selma Moqvist, Weilong Chen, Mathias Schreiner, Feliks Nüske, and Simon Olsson. Ther- modynamic interpolation: A generative approach to molecular thermodynamics and kinetics. Journal of Chemical Theory and Computation, 21(5):2535- 2545, Mar 2025.
[52] Mathias Schreiner, Ole Winther, and Simon Olsson. Implicit Transfer Operator Learning: Multiple Resolutsn Models for Molecular Dynamies. In Oh, Naumann, A. Glober- son, K. Saenko, M. Hardt, and S. Levine, editors, Advances in Neural Information Processing Systems, volume 36, pages
Curran Associates, Inc., 2023.
[53]
Leon Klein, Andrew Foong, Tor Fjelde, Bruno Mlodozeniec, Marc Brockschmidt, Se- bastian Nowozin, Frank Noe, and Ryota Tomioka. Timewarp: Transferable accelera- tion of molecular dynamics by learning time-coarsened dynamics. In A. Oh, T. Nau- mann, A. Globerson, K. Saenko, M. Hardt, and S. Levine, editors, Advances in Neural Information Processing Systems, volume 36, pages 32863-52883 Curran Associates, Inc., files/ /iles/paper//00/ /appr//2233/flle Paper- Conf erence pdf.
2023. URL https: https: //ppoe/ee///ccoc/eeeeeeee wtw.rrreeeedinneeeeeneeirrpepppe neurips cc cttpccc//pppaa cccpaper_.... 3558066788800044444fccc777777668
[54] Sun-Ting Tsai, En-Jui Kuo, and Pratyush Tiwary. Learning molecular dynamics with simple language model built upon long short-term memory neural network. Nature Communications, 11(1):5115, Oct 2020.
[55] J C S Kadupitiya, Geoffrey C Fox, and Vikram Jadhao. Solving newton's equations of motion with large timesteps using recurrent neural networks based operators. Machine Learning: Science and Technology, 3(2)022500, apr 2022.
[56] Xiang Fu, Tian Xie, Nathan J. Rebello, Bradley D. Olsen, and Tommi Jaakkola. Simulate time-integrated coarse-grainee molecular dynamics with multi-scale graph networks. arXiv preprint 2tti:22044.13443 2022.
[57] Guy Dayhoff and Sameer Varma. MImd: Machine learning velocities to propagate molecular dynamics simulations. ChemRxiv, 2025. doi: 10226344hhemrii22033--2222
[58] Sepp Hochreiter and Jürgen Schmidhuber. Long short-term memory, Neural Computation, (8):1735- 1997. ISSN 0899-7667. doi: 10.1 10..162/ccoo.1197998.7777
[59] Diederik P Kingma and Max Welling. Auto-encoding variational bayes. arXiv preprint arXiv:1312.6114.22 2013.
[60] Danilo Jimenez Rezende and Shakir Mohamed. Variational inference with normalizing flows. arXiv preprint arXiv:1505.05770, 2015.
[61] Jascha Sohl-Dickstein, Eric A. Weiss, Niru Maheswaranathan, and Surya Ganguli. Deep unsu- pervised learning using nonequilibrium thermodynamics. arXiv preprint arXiv: /50.003585 2015.
[62] Tianze Zheng, Weihao Gao, and Chong Wang. Learning large-time-step molecular dynamics with graph neural networks. arXiv preprint arXiv:2111.15776 2021.
14

---

[73] Bingqing Cheng, Mandy Bethkenhagen, Chris J. Pickard, and Sebastien Hamel. Phase be- haviours of superionic water at planetary conditions. Nat. Phys., 17(11):1228 12221 November 2021. ISSN 1745-2473, 1745-2481. doi: 10. 100000//115670001-01333444- [63] Fabian L. Thiemann, Thiago Reschützzeger, Massimiliano Esposito, Tseden Taddese, Juan D. Olarte-Plata, and Fausto Martelli. Force-free molecular dynamics through autoregressive equivariant networks. arXiv preprint 2rXiv:2503323374 2025.
[64] Fuchun Ge, Lina Zhang. Yi-Fan Hou, Yuxinxin Chen, Arif Ullah, and Pavlo O. Dral. Four- dimeeniionllspacctiime atomistic artificial intelligence models. The Journal of Physical Chemistry Letters, 14(34):732-7744,, 2023. doi: 10.1021/acs,... 10.1022//cs..ccceeet30.....
[65] Fuchun Ge and Pavlo O. Dral. Artificial intelligence for direct prediction of molecular dynamics across chemical space. arXiv preprint arXiv:2505. 2025.
[66] Minkai Xu, Jiaqi Han, Aaron Lou, Jean Kossaifi, Arvind Ramanathan, Kamyar Azizzade- nesheli, Jure Leskovec, Stefano Ermon, and Anima Anandkumar. Equivariant graph neural operator for modeling 3d dynamics. arXiv preprint arXiv:2401.11053 2024.
[67] Alvaro Sanchez-Gonzalez Jonathan Godwin, Tobias Pfaff, Rex Ying, Jure Leskovec, and Peter W. Battaglia. Learning to simulate complex physics with graph networks. arXiv preprint arXiv:2002.09400 2020.
[68] Yulia Rubanova, Alvaro Sanchez-Gonziie Tobias Pfaff, and Peter Battaglia. Constraint-basee graph network simulator. arXiv preprint arXiv:21/2.09) 2021.
[69] Federico Grasselli, Sanggyu Chong, Venkat Kapil, Silvia Bonfanti, and Kevin Rossi. Uncer- tainty in the era of machine learning for atomistic modeling. arXiv preprint ar arXiv:2503.0919, 2025.
[70] Michele Ceriotti, Jér' me Cuny, Michele Parrinello, and David E. Manolopoulos. Nuclear quantum effects and hydrogen bond fluctuations in water. Proc. Natl. Acad. Sei. U. S. A., 110 (39):15511-15396 September 2013, ISSN 00278424. doi: 10.1073/pnas. 1308560110.
[71] Gregory R. Medders, Andreas W. Götz, Miguel A. Morales, Pushp Bajaj, and Francesco Paesani. On the representation of many-body interactions in water. J. Chem. Phys., 143(10): 104102, September 2015. ISSN 0021-9606. doi: 10. 1933/444449
[72] Mariana Rossi, Michele Ceriotti, and David E. Manolopoulos. Nuclear Quantum Effects in H+ and OH- Diffusion along Confined Water Wires. J. Phys. Chem. Lett., 7(15):0011-3007, August 2016. ISSN 19487185. doi: 10. 1..1021/cc.jjclett.ooo
[74] John P. Perdew, Adrienn Ruzsinszky, Gábor I. Csonka, Oleg A. Vydrov, Gustavo E. Scuseria, Lucian A. Constantin, Xiaolan Zhou, and Kieron Burke. Restoring the Density-Gradien Expansion for Exchange in Solids and Surfaces. Phys. Rev. Lett. 100(13): April 2008.
ISSN 0031-9007, 1079-7114, doi: 10.1 1ee/PhysRevl.ett.. 100.11640.36
[75] Scott Habershon, Thomas E Markland, and David E Manolopoulos. Competing quantum effects in the dynamics of flexible water model. J. Chem. Phys., 131(2):44501, July 2009. ISSN 1089-7690. doi: 10. 10.11/00311677777
[76] Nicola Marzari, David Vanderbilt, Alessandro De Vita, and M. C. Payne. Thermal Contraction and Disordering the AI(110) Surface. Phys. Rev. Lett. 82(16):329-33299, April 1999. ISSN 0031-9007 1079-7114, doi: 10. 20.11033Phyyseeeeen.2...22666
[77] Li-Peng Hou, Xue- Que-QQang Zhang, Nan Yao, Xiang Chen, Bo-Quan Li, Peng Shi, Cheng- -iin Jin, Jia-QU Huang, and Qiang Zhang. An encapsulating lithium polyuude electrolyte for practical
Apr 2022. ISSN 2451-9294, doi: 10. lithium&#x201 3:uuuuu batteries. Chem, /ttpt://////ii.o//////////////00 1016/j chempr 2021 023. 10166j.chemp:2221.1 20.02.2..3..... URL https: //tttts:///ii..ee
[78] Xin Gao, Zhiao Yu, Jingyang Wang, Xueli Zheng, Yusheng Ye, Huaxin Gong, Xin Xiao, Yufei Yang, Yuelang Chen, Sharon E. Bone, Louisa C. Greenburg, Pu Zhang, Hance Su, Jordan Affeld, Zhenan Bao, and Yi Cui. Electrolytes with moderate lithium polysulfide solubility for high-performance long-calendar-lite lithium-sulfur batteries. Proceedings of the National 2023. doi: 10. 10733//pp 2300mma.....020.31..
Academy of Sciences,
15

---

[79] Joseph A. Morrone, Thomas E. Markland, Michele Ceriotti, and B. J. Berne. Efficient multiple time scale molecular dynamics: Using colored noise thermostats to stabilize resonances. Chem. Phys., 134(1):14100) January 201 1. ISSN 00219606. doi:
[80] Lorenzo Gigli, Davide Tisi, Federico Grasselli, and Michele Ceriotti. Mechanism of Charge Transport in Lithium Thiophosphate. Chem. Mater, 36(3):1482-1499) Febrarr 2024, ISSN 0897-4756, 1520-5002. doi: 10..122lacs..hhmmmatr.... chemmater..30022...
[81] Filippo Bigi and Michele Ceriotti. Learning the action for long-time-stey simulations of molecular dynamics. arXiv preprint 20000:0000000068,, 2025.
[82] Filippo Bigi, Joseph Abbott, Philip Loche, Arslan Mazitov, Davide Tisi, Marcel Langer, Alexander Goscinski, Paolo Pegolo, Sanggyu Chong, Rohit Goswami, et al. Metatensor and metatomic: foundational libraries for interoperable atomistic machine learning. arXiv preprint arXiv:2508. 15704, 2025.
[83] Ask Hjorth Larsen, Jens Jørgen Mortensen, Jakob Blomqvist, Ivano Castelli, Rune Chris- tensen, Marcin Dulak, Jesper Friis, Michael N Groves, Bjørk Hammer, Cory Hargus, et al. The atomic simulation environment a python library for working with atoms. Journal of Physics:
Condensed Matter, 29(27):273002 2017.
[84) Yair Litman, Venkat Kapil, Yotam M. Y. Feldman, Davide Tisi, Tomislav Begušić, Karen Fidanyan, Guillaume Fraux, Jacob Higer, Matthias Kellner, Tao E. Li, Eszter S. P6s, Elia Stocco, George Trenins, Barak Hirshberg, Mariana Rossi, and Michele Ceriotti. I-PI 3.0: A flexible and efficient framework for advanced atomistic simulations. J. Chem. Phys., 161(6): doi: 10.1063/5.0215860
062504, August 2024. ISSN 0021-9606, 1089-7690.
[85] Aidan P Thompson, H Metin Aktulga, Richard Berger, Dan Bolintineanu, W Michael Brown, Paul Crozier, Pieter In't 't Veld, Axel Kohlmeyer, Stan G Moore, Trung Dac Nguyen, et al. Lammps-a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales. Computer physics communications, 271:108171, 2022.
[86] Leopold Talirz, Snehal Kumbhar, Elsa Passaro, Aliaksandr Yakutovich, Valeria Granata, Fernando Gargiulo, Marco Borelli, Martin Uhrin, Sebastiaan P. Huber, Spyros Zoupanos, Carl S. Adorf, Casper Welzel Andersen, Ole Schütt, Carlo A. Pignedoli, Daniele Passerone, Joost Vondee,, Thomas C. Schulthess, Berend Smit, Giovanni Pizzi, and Nicola Marzari. Materials Cloud, a platform for open computational science. Sci Data, 7(1):299, December doi: 10.1038/s41597-000
2020. ISSN 2052-4463
[87] A Paszke. Pytorch: An imperative style, high segh-performanc deep learning library. arXiv preprint arXiv: 1rXiv:1912.0170:: 2019.
[88] Diederik Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv preprint
2014.
[89] Ask Hjorth Larsen, Jens Jørgen Mortensen, Jakob Blomqvist, Ivano E Castelli, Rune Chris- tensen, Marcin Dulak, Jesper Friis, Michael N Groves, Bjørk Hammer, Cory Hargus, Eric D Hermes, Paul C Jemnings, Peter Bjerre Jensen, James Kermode, John R Kitchin, Esben Leon- hard Kolsbjerg. Joseph Kubal, Kristen Kaasbjerg, Steen Lysgaard, Jón Bergmann Maronsson, Tristan Maxson, Thomas Olsen, Lars Pastewka, Andrew Peterson, Carsten Rostgaard, Jakob Schiøtz, Ole Schütt, Mikkel Strange, Kristian Thygesen, Tejs Vegge, Lasse Vilhelmsen, Michael Walter, Zhenhua Zeng, and Karsten W Jacobsen. The atomic simulation environment- a Python library for working with atoms. J. Phys.: Condens. Matter, (9727)2773002, July 2017. ISSN 0953-8984, 1361-648X. doi: 10. 111118//1311-65800 6488-660
[90] Nicholas Metropolis, Arianna W Rosenbluth, Marshall N Rosenbluth, Augusta H Teller, and Edward Teller. Equation of State Calculations by Fast Computing Machines. J. Chem. Phys.,
21 216)::088711092 1953.
[91] Benedict Leimkuhler and Charles Matthews. Robust and efficient configurational molecular sampling via Langevin dynamics. J. Chem. Phys., 138(17), 2013. ISSN 00219606. doi: N0.1063/1.4802990
16

---

[92] Giovanni Bussi, Tatyana Zykova-Timma Timan, and Michele Parrinello. Isothermal-isobaric molecular dynamics using stochastic velocity rescaling. J. Chem. Phys., 130(7):074101, February 2009. ISSN 0021-9606.
doi: 10. 10.1006313000 (003/88...
[93] Alex Kendall and Yarin Gal. What uncertainties do we need in Bayesian deep learning for computer vision? In NewrlPS, 2017.
[94] Michael F. Herbst, Antoine Levitt, and Eric Cancès. A Posteriori enror estimation for the non-self-consistent Kohn-Sham equations. Faruday Discuss. 224:227-246, 2020. ISSN 1359-6640, 1364-5498.
doi: 10. 1039/DDFDD0004S..
[95] Michael F. Herbst and Thomas Fransson. Quantifying the error of the core-valence separation approximation. J. Chem. Phys., 153(5):0541 14, August 2020. ISSN 0021-9606, 1089-7690. doi: 10. 00.1063/5.0013538.
[96] Michael Herbst and Antoine Levitt. Black-box inhomogeneous preconditioning for self- consistent field iterations in density functional theory. J. Phys.: Condens. Matter, 33(8): 085503, February 2021. ISSN 0953-8984, 1361-648X. doi: 1..1088/1311-648X280 6488//bbbbbb.
[97] Alexander Immer, Emanuele Palumbo, Alexander Marx, and Julia Vogt. Effective bayesian heteroscedastic regression with deep neural networks. Advances in Neural Information Processing Systems, 36:53996-54000 2023.
[98] Balaji Lakshminarayanan, Alexander Pritzel, and Charles Blundell. Simple and scalable predictive uncertainty estimation using deep ensembles. NIPS, 2017.
[99] Jasper Snoek, Oren Rippel, Kevin Swersky, Ryan Kiros, Nadathur Satish, Narayanan Sun- daram, Mostofa Patwary, Mr Prabhat, and Ryan Adams. Scalable Bayesian optimization using deep neural networks. In ICML, 2015.
[100] Agustinus Kristiadi, Matthias Hein, and Philipp Hennig. Being Bayesian, even just a bit, fixes overconfidence in ReLU networks. In ICML, 2020.
[101] Erik Daxberger, Agustinus Kristiadi, Alexander Immer, Runa Eschenhagen, Matthias Bauer, and Philipp Hennig. Laplace redux-effortless bayesian deep learning. Advances in neural information processing systems, 34:20089-20103 2021.
[102] Filippo Bigi, Sanggyu Chong, Michele Ceriotti, and Federico Grasselli. A prediction rigidity formalism for low-cost uncertainties in trained neural networks. Mach. Learn.: Sci. Technol., 5(4):045018, December 2024. ISSN 2632-2153, doi: 01088/2662-2153/aa880
[103] Victor Garcia Satorras, Emiel Hoogeboom. and Max Welling. B(n) equivariant graph neural networks. arXiv preprint arXiv:2102.09982 09844. 2021.
[104] Mahdi Hijazi, David M. Wilkins, and Michele Ceriotti. Fast-forward Langevin dynamics with momentum flips. J. Chem. Phys., 148(18):184109 May 2018. ISSN 00219606. doi: 10. 1063/1....00 .5029833.
[105] L. Martínez, R. Andrade, E. G. Birgin, and J. M. Martínez. Packmol: A package for building initial configurations for molecular dynamics simulations. Journal of Computational Chemistry, 30(13)2157-2164, 2009. doi: https://doi.org/100 http://d....rr.....00000/2.2222
[106] Zoran Šućur and Vojtech Spiwok. Sampling enhancement and free energy prediction by the flying gaussian method. Journal of Chemical Theory and Computation, 12(9):464444650, 2016. doi: 10.1022/aes..cte..6660 6bb.05555
[107] Michele Ceriotti, Giovanni Bussi, and Michele Parrinello. Langevin equation with colored noise for constant-temperature molecular dynamics simulations. Phys. Rev. Leut 102(2): 020601, January 2009. ISSN 00319007. doi: 10.1 10.100/hhysev..eett 10..1l.020.100020000000 102.020601.
[108] Laboratory of Computational Science and Modeling (COSMO), EPFL. UPET: Univer- sal interatomic potentials for advanced materials modeling. URL https: //ggt://///ttttt///// com/ lab- cosmo/upet. [Online; accessed Jan. 2026].
17

---

[109] Luis Barroso-Luque, Muhammed Shuaibi, Xiang Fu, Brandon M. Wood, Misko Dzamba, Meng Gao, Ammar Rizvi, Lawrence Zitncck, and Zachary W. Ulissi. Open materials 2024 (omat24) inorganic materials dataset and models, 2024. URI https: / /arxivvx.vv////////////x///v/av/ org/abs/ haabs/2410... 12771.
[110] Aaron D. Kaplan, Runze Liu, Ji Qi. Tsz Wai Ko, Bowen Deng, Janosh Riebesell, Gerbrand Ceder, Kristin A. Persson, and Shyue Ping Ong. A foundational potential energy surface dataset for materials, 2025, URL https: arxiv org/ abs/ab///2/ /200 2503.04070. 04070.
18

---

A Model details
A.1
Why predict gf 9, and p': small- and large-time limits
should be noted that the models for liquid argon in the early work by Zheng et al. [62] effectively predict qf 9, P,At/m; and pi P;- While the additional terms with respect FlashMI) (-p,At/m; and - P, respectively) reduce the variance of the targets for small time steps and/or very smooth potential energy surfaces (such Lennard-Jon for more complex systems and larger predicted time steps, which are the focus of the present work. This is consequence of more complex systems generally having smaller position and momentum correlation times. As a result, we do not shift the targets by these additional terms in our work.
argon in Ref. [62]), they instead increase it
A.2 Predicting mass-scaled positions and momenta
All models shown in this work are trained on, and therefore predict, mass-scaled displacements and momenta defined as Λq,- Aq; /m, and P; = P:/ ym,, respectively, where m; is the mass of atom i. This is aimed at making the scales of the training targets uniform across atoms of potentially very different mass, and it is of fundamental importance for models trained on the whole periodic table. This prevents, for example, the displacement of light atoms or the momenta of heavy atoms from dominating the loss, instead leading to good predictions for all atoms, independent of their mass. At prediction time, simple scaling using the masses recovers the conventional displacement and momenta. We found that a similar, but data-driven, approach can provide the same benefits. This consists of using different standardization factors for different chemical elements, so that training displacements and momenta are scaled to unit standard deviation during training. for each chemical element (i.e., two scaling factors used for every single chemical element in the dataset: one for displacements, one for momenta).
A.3 Center-of-mass enforcement
Using Eqs. 2 to evolve a system without external forces naturally leads to conservation of total momentum of the system, i.e.,
N
Σp'
p
=. =.
Since the total momentum constant, the center of mass of the system follows a uniform linear motion, i.e.,
N
N
Σ mag
Σ m.g. - At
Σ P
(5)
1=1
1=1
1=1
Both conditions enforced within the model avoid unphysical drift effects during molecular dynamics simulations (Fig. 1). We note that removing the center-of-mass motion entirely would not be correct in the general case, although many MD simulations performed with this additional constraint. Although we enforced these contraints within the model in this work, we recommend applying them inference time only in order not to break the locality assumption that underpins our approach.
A.4 Optional inference- filters
Within FlashMD, we have implemented "filter" (Fig. 1) that can be employed at inference-time to mitigate the artifacts of direct MD prediction. We refer to the dedicated sections for energy conservation enforcement (App. C) and thermodynamic ensemble control (App. D). Here we only provide a discussion of the random rotation filter.
Since equivariance is not exactly preserved and only learned via data augmentation in the case of unrestricted architectures such as PET [18]. simulations performed with the resulting FlashMD model would be prone to spurious effects. To mitigate this, we adopt the strategy proposed Langer et al. [20], where random rotations the simulated system are performed to average out any artifacts of non-equivariance along the MD trajectory. Implementation details of the random rotation filter is shown in Fig. We note that case of GNNs that preserve rotational symmetry, this filter is not needed and would actually have effect if applied.
19

---

(pi(t)
R
(T)
(ptor(r Ar)]
R
(pi At))
FlashMD
(when using non-eq. GNN)
(q:(t))
{qi t (T))
{qtt (t AT)3
+
Figure 4: Implementation of the random rotation filter. A random rotation matrix R is sampled and applied on all coordinates and momenta before the rotated inputs are provided to FlashMD. After model inference, R applied to rotate the system back to the original coordinate reference. Random rotation filter is only relevant for rotationally unconstrained GNNs.
B Training details
B.1 Loss function, optimization, normalization
All models are implemented in PyTorch [87] and, unless specified otherwise (see App. E), are trained using a loss function given by the sum of two mean square error terms, for the mass-scaled momenta and the positions respectively:
Nerain
N Σ
C
Σ
Po + (Aq) - A P
(6)
3N
where 8 is an index for structures in the training set and N, is the number of atoms in structure 8. In order ensure similar weight in the loss function between position and momentum terms, the mass-scaled positions and momenta are scaled by their standard deviation across the dataset before training.
Optimization carried out using the Adam [88] optimizer with an initial learning rate 3 10 10-4, Learning rate decay is applied at regular intervals of 100 and 50 epochs for the water and universal models, respectively. Training-tiin rotational augmentation for vectorial targets is carried out the same way in Pozdnyakov and Ceriotti [18].
B.2 Training- time energy conservation
We found that the degree of energy conservation on structures of the validation set correlates well with the quality of the models during molecular dynamics runs. As a result, during training, choose the best model as the one having the lowest product of three terms, evaluated across the validation set: (i) RMSE on the predictions of p'. (ii) RMSE the predictions of q' and (iii) RMSE on the energy of the resulting structure when compared to the energy of the target structure. While an energy term might also be included in the loss function, found that it slows down training significantly (due to evaluations of the energy model and its gradients), without improving the quality of the FlashMD models in any measurable way.
Indeed, as shown in App. models with similar accuracy on positions and momenta can predict MD states with highly varying degrees of energy conservation. In particular, energy misalignment is often observed if the error on the energies is ignored. Although we found this approach to improve the quality of the simulations afforded by our water models (both PET-MAD and q TITIPPP/)) we found its impact on universal models to be less dramatic.
We also found it useful to compare errors in total energies with familiar metrics such as "chemical accuracy" or the accuracy of the underlying energy model. For all models tested this work, such comparisons correlate extremely well with the quality of the models in the resulting MD simulations.
B.3 Reference MD trajectory generation
All reference MD trajectories were obtained from simulations performed with PET-MAD [39] (version 1.0), universal MIJIP capable of making reasonably accurate predictions of the potential energy surface across the entire periodic table of elements. All simulations were performed using the Atomic Simulation Environment [89] (ASE) software (version 3.24).
20

---

Water-specific models A water structure at experimental density (at 298 K and 1 atm) was equilibrated with PET-MAD (or q-TIP4P/E for the q-TIP4P/f-based water models discussed in Appendix J). Subsequently, two more structures were generated by increasing and decreasing the volume of the cell by 10%, scaling the atomic positions accordingly. For each of the three resulting structures, NVT equilibration runs were performed at all temperatures between 20 and 1000 K, in steps of 20 K, with a time step of 0.5 and a duration of 5 ps, using a Langevin thermostat with characteristic time of 10 fs. Subsequently, each equilibrated structure used to produce an NVE MD trajectory of ps with a time step of 0.25 fs. Structures for training were extracted from these trajectories every 100 fs, and augmented with their time-reversed version, for total of 5400 structures.
Universal models 10,000 structures from the MAD dataset, used in the training of PET-MAD [39], baseline MIIP, were randomly selected for reference MD trajectory generation (see Ref. [39] for further details on the MAD dataset). The initial geometry was first energetically optimized with the BFGS algorithm until the maximum force component threshold of 0.01 eV/À reached. The energy-optimized system was put through equilibration under the NVT ensemble for 10 ps with timesteps of 0.5 fs. characteristic time of 100 fs was used the Langevin thermostat. The final configuration from NVT equilibration was then taken for trajectory production under the NVE ensemble for 2.5 ps with finer timesteps of 0.25 fs.. Positions and momenta were recorded every timestep for FlashMD training. Simulations were repeated 10 times for each structure, with a randomly selected temperature between and 1500 K. Structures for training were chosen from these trajectories every 500 fs (5 samples per VVE trajectory to avoid correlate samples), and augmented with their reverreversed version, for total of million structures.
C Enforcing conservation of energy by momentum rescaling
Especially for NVE simulations, it important to avoid excessive energy drift. In practice, we found that enforcing total energy conservation after each step is beneficial (even for thermostatted runs, Sec. 4). This is achieved by rescaling the momenta in the following way:
E' E K"
p op', a= a-
where E' is the total energy after the step, E is the total energy before the step, and K' is the kinetic energy after the step.
should be noted, however, that enforeing energy conservation requires one or two energy evaluations per step (depending on the ensemble and integration scheme), potentially introducing significant overhead. Implementation details the energy conservation enforcement filter in FlashMD is visualized in Fig. 5.
(P)
Ip((t) 2m,
(T)
(((g((t)) E(x)
Energy model
ta P:(T AT))
(q:(1 + Ar))
E(r dr)
(p((r+An))
F(r FCO 1. /2//
(p. + Ar))
Figure 5: Implementation details of the energy conservation enforcement filter in FlashMD. Energy model can be any model of the interatomic potential (e.g. MIJP, classical force field, etc.) that can be used compute V((g,).
21

---

D Arbitrary thermodynamic ensembles with FlashMD
Molecular dynamics trajectories conserve, at least approximately, the classical total energy of the system, which makes it appropriate for sampling configurations under constant energy, volume and particle number (NVE) conditions. Several modifications to Hamiltonian dynamics have been proposed [25, 26] to generate configurations consistent with other thermodynamic conditions, such as constant temperature (N (N),), constant pressure (NpT). or constant chemical potential (nVT). These ensembles are often more relevant to compare with experimental conditions, that usually do not involve closed, isolated setups especially not on the small length scales that are used in simulations. Here we discuss two approaches that are routinely used to this end in MD, and how they can be can be combined with FlashMD to extend its constant-energy long-sttrn mnngsttriee integration to sample other ensembles.
First, given that NVE trajectories conserve only the total energy, but also the probability measures associated with most other ensembles, is possible to alternate segments of N trajectories with discrete Monte Carlo moves changing the particle velocities, the simulation cell size, or the nature of the atoms, using a Metropolis- Hastanps criterion [90] to accept or reject them in a way that is consistent with the desired ensemble. This approach can be applied straightforwardly and its reliability depends on the assumption that segments of FlashMD trajectories are symplectic and energy-pre preserving
FlashMD,
high degree of accuracy, which is why we give much emphasis these diagnostics in our study.
The second approach is slightly more subtle and requires some additional technical background, and we discuss and test in more detail. Integrators for Hamiltonian dynamics can be expressed in a Liouvillian formalism, in which the trajectory density P(q. p) is evolved in time according to operator that combines the time evolution of the different variables, e.g. for Hamiltonian dynamics
iL = OH am ØH an ity ttp
(8)
Öp. 0q, dq, op.
Finite-time propagation of the trajectory density can be formally achieved with an exponential operator ci +.y)) and the difficulty in developing an exact propagation algorithm for (q, p) can be understood as a consequence of the fact that L, and Lp do not commute. The error grows with the time step At, and can be reduced using symmetric Trotter factorizations such as At/2g, Af/..,At//2 This splitting corresponds, in the trajectory picture, to the VV integrator in Eq. (2). Continuous equations of motion that describe other thermodynamic ensembles can be derived with an extended Lagrangian formalism, in which additional structural parameters (e.g. the simulation cell volume V or shape) are associated with fictitious masses and momenta. Starting from the Lagrangian, one can derive a Liouvillian that describes their time evolution as an additional term, e.g. L.v.
Langevin-type thermostats [29] can also be described with an associated Liouvillian LG and are good example to explain the formalism. There are in fact multiple possible ways to factorize the overall Liouvillian: the so-called OBABO splitting reads 2 At / and corresponds to bracketing a velocity Verlet integrator (BAB) between two finite-time propagators for Uhlenbeck process (O, the Langevin equation for a free particle with inertia)
P; t >
p, + mm,kTT(1 - e-
1 10V
P: P: 2 åqq
At
4, P. At
(9)
mi
1 @V
P: < 2 0q:
At
P; 7e-1211 + + m,kkT(1 ---nat)
where & and & are vectors of uncorrelated, unit-variance Gaussian random numbers. Note that this splitting preserves the symmetry of the VV integrator. Many other splittings are possible, such as BAOAB (with the Uhlnebeek propagator sandwiched between two half-VV integrators), which has been found to be more accurate for position observables [91].
22

---

From this extensive overview of the Liouville operator formalism and its connection to the theory of integrators, one can see how to incorporate FlashMD into the integration schemes that are used for other thermodynamic ensembles. the full Liouvillian for a given integrator is É + I.g + Exps one can factor integration over Ar as
ei A/ (e) (l.. A A1//22
(10)
The central term is precisely the evolution that FlashMD aims approximate, and can therefore be readily replaced with one large step on (q.P).
There few considerations that should be made when designing one of these extended integrators for FlashMD. First, if there are multiple splittings available for the base integrator, one has to choose those that involve VV core. For instance OBABO can be used, but not BAOAB, and we cannot use the Mussi--yyovv-PParineelll splitting [92], but a more naive one that does not simultaneously update atomic positions and cell vectors. Second, if one wants to further split ل a At/2, it should be done in symmetric way, applying the factors in opposite order before and after the FlashMD step. Last, and most importantly, the integration of çiL'A+=1/2 should be accurate also with a large time step, and the factorization with the VV core be similarly accurate, or at least preserve the target ensemble. This implies choosing large effective masses for extended Lagrangian terms (e.g. the cell volume in a prressntt integrator). Long time scales for Langevin-type thermostats should also be chosen if one wants to preserve time-dependent properties of the original Langevin dynamics but this is lesser concern for sampling accuracy, because usually Langevin-type free-particle integrators preserve the velocity distribution for any time step. The workflow proposed herein for the integration of thermostats and barostats with FlashMD is shown in Fig. 6.
(p)(t).q(1))
(p + Ar).q((G Ar)]
(.(T)
Thermostat/ Barostat
p(((T
p1(T AT)
Thermostat/ Barostat
pT AT)
FlashMD
(T())
(T(
Aqi(T + AT)
Aq (I + At)
Figure 6: Integration of thermostats and barostats with FlashMD for thermodynamic ensemble control. Note that the integration is performed for the full stride in FlashMD, and the split operators both ends are applied for half strides.
E
Uncertainty quantification
ML models are inherently statistical, introducing different types of error at prediction time. In the context of FlashMD, uncertainty quantification (UQ) becomes an even bigger necessity as the potential error accumulation along a trajectory makes FlashMD simulations prone to many undesirable artifacts, leading to the incorrect sampling of the thermodynamic ensemble, which can manifest itself, for example, as unphysical bond forming/breaking behavior, uncontrolled expansion of the simulated system, etc. In this section, we identify the design choices that are required for robust UQ for MD prediction models and provide simple demonstrations.
In general, there largely exist two different sources of uncertainty for the model: aleatoric, irreducible uncertainty stemming from the "noise" in data, and epistemic, reducible uncertainty from the model's lack of knowledge [93]. In the training of MLIPs, it is widely assumed (although not necessarily true [94-96]) that the reference data is noise-free, and that is therefore appropriate to only account for epistemic uncertainty in UQ approaches for MLIPs. In the learning formulation for FlashMD that directly targets evoeleeolee positions and momenta, we note the potentially significant presence of aleatoric uncertainty due to the chaoticity of the underlying physical problem. This aleatoric contribution to the uncertainty is expected to be strongly heteroscedastic, since different chemical systems exhibit very different degrees of chaotic behavior in molecular dynamics (in other words, the Lyapunov exponent can vary significantly based on the system, see See. 2.2 ). Last but not least, the UQ scheme of choice should not result in significant prediction time overhead in the simulation, as
23

---

-= 16 fs
== 64 fs
101
10
10"0
18"1
19"7
19"1
10-
13-1
10.
13:1
Actual Error
Actual Error
12-1
13-1
12'4
(L))
OLLPR (p))
1214
OLLPR
OLLPR (p:)
OMVE
( (p))
MVE (q))
OMVE (p)) 10
10
20
10
10
10"
101
10
20
10 30
10
20
10'
131
10
10
Predicted Error
Predicted Error
b
training domain
Training domain
(ewrldj))
owr(())
al
10
(urn(G.)
OLux(A)
0.8 1.0
1.2
1.4
0.6
0.8
1.0
1.2
1.4
Cell scaling factor
Cell scaling factor
Figure 7: Uncertainty quantification diagnostics for the water-specific FlashMD model its test set. (a) shows the uncertainty plots (predicted variance on the axis vs squared residual on the axis) on mass-scaled positions and momenta (q; and P,) for two models trained the water dataset (with 16 fs and 64 fs strides, respectively). Units are eV for (f,) and A for (q;). ØLLPR is shown in green, and MVVE is shown in orange. The black lines corresponds to the parity line, as well as pairs of iso probbabbablly lines of the ideal distribution containing density equivalent to that contained within la, 2o, 30 Gaussian distribution. (b) shows the predicted uncertainty, in the same units, for diitr--distiitution predictions using the 16 fs model, function of the scaling of the cell and the atoms within it.
this would undermine the key advantages of FlashMD, For these reasons, we adopt UQ scheme that quantifies both types of uncertainties with a near-zero computational overhead. The method is sketched in the blue inset of Fig.
Taking inspiration from Immer et al. [97], we assume the overall uncertainty to arise from a sum of an epistemic and aleatoric term
2= +
(11)
which are predicted by different types of UQ estimator. For the aleatoric component, the prediction heads modified yield mean-variance estimators (MVEs), which the model predicts mean and variance of the target predictions. FlashMD in this mode is trained to the negative log-likelihood loss
1 (0 [Ino? L f-
(y
Frof)
(12)
2
o2
where of is the variance afforded by the mean- variance estimator, parametrized as described in Lakshminarayanan et al. [98]. In our case, the overall loss is obtained by summing one of such terms for mass- mass-scaled positions and one for mass scaled momenta. A Laplace approximation is usually considered good model for the epistemic uncertainty o and we implement as a last-layer approximation (LLPR) [99-102] on the mean part of the MVE, i.e., on y.
We now demonstrate the UQ capabilities of FlashMD, with a special focus on the need for aleatoric uncertainty within direct MD prediction models discussed in Sec. 3. To do slightly deviate

---

from Eq. (11) and consider separately the MVE and LLPR uncertainty predictions for liquid water and analyze their behavior (Fig. 7a). The ÅT 64 fs model, which very difficult to learn for the liquid water system (due to the fast correlation times of the physical system), displays the failure of epistemic uncertainty in isolation, exhibiting a narrow spread in values that does not provide meaningful insights. In contrast, the mean-variance estimator provides good uncertainty estimates. Perhaps more surprisingly, uncertainties of the 16 fs model are predicted reasonably by both uncertainty estimators. Since Åt = 16 fs is more learnable regime for water, where epistemic uncertainty is expected to dominate, this implies that the mean-variance estimator is also capable of capturing epistemic uncertainties to good accuracy, at least within the training distribution. As an out-of out-of-domain example, Fig. 7b shows the average LLPR and MVE uncertainties as a function of the compression (or expansion) factor of water cell. Even for very compressed or very stretched cells (up to a change of 50% in the cell length), the mean-variance estimator provides a qualitatively correct uncertainty profile, suggesting that MVE alone might be sufficient to quantify uncertainties in direct MD predictions. We leave a more thorough analysis of combining the two UQ metrics for future work.
F Ablation studies
10
10
1.25
RMSE (AVo]
10
RMSE (VeV)
10-1
RMSE [meV/aL]
10a
Symplecticity
1.00
10a
0.75
10:2
10a
10
10-1
109
0.25
water
10
10-3
10-2
0.00
universal
16
64 256
16
64
256
256
16
Time step [ts]
Time step [fs]
Time stop
Time step [fs]
Figure 8: Comparison of validation root mean square errors (RMSEs) in (a) mass-scaled positions, (b) mass-scaled momenta, (c) energy conservation, and (d) symplecticity as function of the stride for different FlashMD models. All errors are calculated on the respective validation sets, except for the symplecticity errors, which are evaluated using the left-hand side of Eq. 3 for 100 degrees of freedom of a liquid water structure. Dotted line in (c) is the RMSE of PET-MAD, an indirect metric of energy accuracy, and in (d) marks perfect symplecticity.
F.1 Effect of the stride length on training
As result of the chaoticity effects described in Sec. 2, it is desirable to examine the errors in the training as the predicted time stride increases. Fig. shows the errors on energies and symplecticity, as well as mass-scaled positions and momenta, for the training runs on the water and universal datasets for different time steps. As expected, training becomes extremely difficult after a certain number of steps. From these plots, it would appear that the increase in the machine learning error is not exponential with the time step, but rather polynomial. While the interpretation of this observation is not trivial, it is potentially promising improving the accuracy of the models using more data and/or learnable parameters, without the presence of hard limits to the accuracy. However, the rapid increase in the non-symplectic behavior of the predictions is alarming and worth future investigation.
it would make longer time scales accessible by simply
F.2 Enforcing energy conservation
We will now illustrate the effect of the procedure to enforce energy conservation presented in App. C on simulation targeting the NVE ensemble. Fig. shows that, while FlashMD exhibits a total energy drift that would quickly lead to unstable trajectories and large sampling errors, the proposed energy conservation enforcement technique allows for exact conservation of the total energy along the simulation, producing stable trajectory (even thought it is not guaranteed to sample the NVE ensemble because not exactly symplectic). We further show in Sec, that this technique improves the temperature control in NVT simulations.
25

---

Without E conservation filter
With conservation filter
2
Energy [eV]
potential
total
0.0
0.5
1.0
1.5
2.0
2.5
3.0
3.5
4.0
0.0
0.5
1.0
1.5
2.0
2.5
3.0
3.5
4.0
Time [ps]
Time [ps]
Figure 9: Comparison of the potential, kinetic, and total energies along short trajectories of water FlashMD model with 4 fs strides, with and without the energy conservation enforcement filter. The mean potential/kinettcchoaall energy is subtracted from all energy profiles.
-4.880
Checkpoint A
4,880
Checkpoint B
FlashMD energy [eV]
FlashMD energy [eV]
4.885
-4.885
-4.890
-4,890
-4.895
44. -4,900
4,900
-4,905
-4.900
-4.895
-4.880
-4.895
-4.890
-4.885
PET-MAD energy [eV]
PET-MAD energy [eV]
Figure 10: Parity plots of the predicted energy for two water-specific model checkpoints belonging to the same training The energies of the structure predicted by the model the left misaligned.
F.3 Energy criterion for model selection
App. B deseribes how the error on the energy of the predicted structures is also used during training to choose good models. Here, illustrate the utility of such approach, show models from different epochs in the same training run (a water-specific model trained to predict over a time stride of 16 fs). Despite the two models having similar position and momentum errors, the energies of the first model are misaligned (Fig. 10). In simulation, such an offset in model prediction would induce progressive and unphysical cooling of the system.
G Timings
In this appendix, we report the overall timings for the simulations that were performed in this work. When multiple FlashMD models with different time steps were used, we report the largest-stride model that affords qualitatively accurate results. The timings obtained in this way are compiled in Table 2. In comparing these results with the theoretical speed-ups given by Ar, one should keep several issues into considerations: (i) We did fine tune the time step of the base MD runs, nor attempt
grid of acceleration factors for the FlashMD models. This might skew results in either directions: for instance, water can be run stably with 0.5 fs MD, and LiPS could have probably been run stably with 32 fs FlashMD model. (ii) We did not systematically explore the Pareto frontier of FlashMD models in its architecture setup. However, in the current version, one evaluation of FlashMD is faster than MLIP energy evaluation due to hyperparameter differences. (iii) Some the extensions require additional energy evaluations; for instance energy scaling or NpT integrator require one or two energy evaluations each. There are obvious optimizations, such as only rescaling energy every few steps, which we did not consider order to keep the analysis clearer for this first demonstration of universal FlashMD model.
26

---

Table 2: Timings, strides, and acceleration factors of FlashMD compared to conventional MD on the systems investigated in this work. For completeness, we also include the speed-up factors we obtain without energy conservation enforcement (ECE), which would slightly degrade the quality of the simulations.
System (# atoms)
Liquid water (192)
AI(1 10) slab (560)
Solvated alanine dipeptide (622)
LigPS, (768)
MD timing [stride]
2.0 10 [0.25 fs]
3.2 104 [s]
10 [0.5 [s]
1.4 10 [2 fs] 4.1 [16 fs]
FlashMD timing [stride]
4.0 102 (16 [16
5.1- 102 [64 fs]
1.8 10" [16
Acceleration factor
50
48
12#
3.4*
Acceleration factor (no ECB)
195
186
12a
3.4"
*These simulations are run in the Np ensemble, leading to more energy model evaluations in order to compute stresses. Re-using these energy evaluation for energy conservation enforcement and printing the energy to output, opposed to recomputing them, would reduce the overhead, but it is not exploited in the present implementation. Furthermore, using non- conservative stresses [21] would speed up these simulations by factor between two and three.
H Liquid argon benchmark (MDNet)
The early work by Zheng al. [62] presents MDNet, simple architecture for fitting molecular dynamics trajectories in the NVE ensemble. Here, we compare FlashMD against MDNet, as well as Equivariant Graph Neural Networks [103] (EGNN, which was also explored as an alternative in Zheng et al. [62]). The lack of code and sufficient description makes it difficult to reproduce the workflow of the authors exactly; nonetheless, we attempt a similar set-up Zheng al. [62]:
All reference simulations are run using LAMMPS, using a system of 256 argon atoms and a Lennard-Jones potential. The time step for the reference simulations is 1 fs.
Ten runs are performed (eight for training, one for validation, one for testing), initially equilibrating in the NpT ensemble for 100 ps, and then performing production run for 10 in the NVE' ensemble using the velocity Verlet algorithm.
From each NVE trajectory, 25 equally spaced configurations are selected to be part of the training/validation/ees set.
A stride of 128 fs is used for training. Even though smaller strides are also explored in Zheng al. [62], we find that this physical system is not particularly challenging due to its very simple and smooth potential energy surface, and therefore we only test predictions the largest stride investigated in the original publication. Time-reversed targets are also added to the dataset, following Zheng et al. [62].
Using this set-up, we were able to reproduce the large-stride velocity Verlet results in Zheng et al. [62] exactly. The accuracies of velocity Verlet, EGNN, MDNet and FlashMD are shown in Table 3, where FlashMD is shown to outperform all methods. Note that our set-up involves less than one tenth of the training data used in Zheng et al. [62], and that FlashMD's accuracy is likely to be underestimated as result.
Table 3: Accuracies of different methods for molecular dynamics trajectory predictions, on a liquid argon system with a predictive stride of 128 fs. EGNN and MDNet results are from Zheng et al. [62]. Positions errors are in units of A, velocity errors are in units of Å/fs.
Method
Velocity Verlet
EGNN
MDNet
FlashMD
RMSE (g)
3.9 10
1.3 10 2
2.3 10 3
5.4 10
RMSE (v)
1.8 10-3
1.9
10-1
5.7 10-5
8.4 10-6
I SPC/E water benchmark (TrajCast)
27

---

TrajCast [63] published three datasets for the direct learning of molecular dynamics. Here, we compare our accuracies with the accuracies reported Thiemann al. [63], using the SPC/E water dataset they provide.
Table 4: Test-set accuracies of FlashMD and TrajCast [63] when trained SPC/E water dataset [63]. TrajCast results are reproduced from [63]. All errors are given in percentage MAE.
Architecture
TrajCast
FlashMD
Displacements 0.17
0.17
Momenta
0.37
0.22
The FlashMD results here make use of a slightly improved architecture compared to the one used to generate the results shown in this work. The same architecture was used to train the r'SCAN FlashMD models we currently recommend, and it corresponds the FlashMD implementation in metatrain [82]. The architecture used to produce the results in this work, for example, yields momentum MAE of 0.52%, indicating that FlashMD and TrajCast have similar accuracies and that relatively minor tweaks can tip the numbers in favor of one or the other. In general, we always found our models to show clear signs of underfitting, showing that larger models and/or longer training times are generally beneficial when training models for the direct prediction of molecular dynamics trajectories, especially when compared to machine-learned interatomic potentials. The design of more accurate models will be crucial challenge to achieve near-quantitatiee results using direct models for molecular dynamics in the future.
J Water simulations based on the q-TIP4P/f model
Due to the inaccurate description of water by PBEsol (the DFT functional used in the training of PET-MAD), also train models the T-TPP4P/f empirical water model [75] to investigate time-dependent properties in liquid water without raising the temperature, which in turns produces artifacts such frequent bond dissociations that significantly affect the dynamies.
The dynamical properties we focus on are the mean square displacement (MSD) of oxygen atoms, as well as the dipole-dipole correlation function, both function of time. In order to avoid the large temperature deviations shown in Sec. 4 for the SVR thermostat, we instead use a fast-forward Langevin [104] thermostat, which is modification of the Langevin thermostat aimed at reducing the effect of Langevin dynamics on dynamical properties, while being applied locally to each atom.
Fig. 11 shows the MSD and dipole-dipole correlation function for MD run with the q- T-TIP4P/f model, as well as FlashMD models of various strides. The results, shown in Fig. 11, establish the need for strong local thermostatting (time constant Tı < 100 fs) in order to obtain consistent statistical sampling between MD and FlashMD. Unfortunately, such thermostatting leads to an underestimation of the diffusion coefficient of water (which is proportional to the slope of the MSD curve) which evident comparing the reference MD results for the Ty = 100 fs with the gentler T = 1000 fs. The tradeoff between thermostat strength and aceuracy in enforcing correct sampling is similar to what was observed for explicit MD simulations based on non-conservatio This suggests that, as in that case, improving model accuracy and refining the thermostatting strategy might mitigate but not cure the underlying problems, and that one should also investigate methods to recover some of the conservation laws that obeyed by MD trajectories, in order achieve guaranteed quantitative accuracy in dynamics and sampling.
direct-force models [21].
K Simulation details
All MD simulations in Sec. 4 were performed with i-PI [84], employing the internal implementation of the thermostats and barostats in the case of reference MD, and using custom, FlashMDD-compatible integrators that adopt the integration schemes explained in Sec. D. In FlashMD simulations, inference was made with both the energy conservation enforcement filter and the random rotation filter, unless specified otherwise.
28

---

T, 100 fs
Tt. 11000 fs
1.0
0.8
0.6
15
16
2.6
f2.f%
MSD [A+]
666
KD
10
0.0
2.5
5.0
10.0
2.5
5.0
7.5
10.0
Time [ps]
Figure 11: Dipole-dipole correlation functions (top) and mean-square displacement (bottom) for simulation of water using fast-forward Langevin thermostats with time constants = 100 fs (left) and TL 1000 fs (right). The stronger thermostat leads to results from the different models that approach those of VV trajectories, but introduces its own spurious effects on dynamics.
Water All PET-MAD-basse water-specific model rums were performed for 100 ps using periodic box of 64 water molecules and starting from structure equilibrated with PET-MAD the NVT ensemble at 450 K. using a volume corresponding to the experimental density of liquid water at 300 K. NpT trajectories were started from the same structure, using a temperature of 450 and pressure of bar. The q-TIP4P/f-basse models were equilibrated and run in the NVT ensemble in same way, but 300 K and using the q- q-TIP4P/f for both equilibration and energy evaluations.
Solvated alanine dipeptide NpT simulations at 450 K and bar were performed with a single dipeptide solvated by 200 molecules water cubic cell of 20 20 20 A as done in Morrone al. [79]. The initial configuration of the system randomly initialized with packmol [105], and the starting conformations of the dipeptide were taken from Ref. [106]. Simulations were performed for PET-MAD MILP at 0.5 fs strides, universal FlashMD with 8 fs strides, and universal FlashMI) with 16 fs strides, for total duration of ns. Langevin thermostat was coupled to the system with 100 fs, and isotropic Bussi-Zykovar kova-Parrinello (BZP) barostat [92] used with fs. including Langevin thermostat coupled to the cell parameters with = 200 fs. For each MD engine, 10 simulations were parallelly performed with different starting configurations to optimize sampling. Despite the difference in time strides, equivalent number of snapshots were sampled across the simulation setups ensure fair comparison of the resulting free energy surfaces.
400
AI(110) surface AI(110) slab configurations were generated with ASE from a 5 x 8 superrell and 20 Á vacuum direction. Following Marzari al. [76], NVT simulations were performed for the system for the temperature range between 400 K and 900 K, for 0.5 ns. The system was coupled to SVR thermostat [30] with =10 fs. PET-MAD MLIP at strides, universal FlashMD at 16 strides, and universal FlashMD 64 fs strides were used for the simulations.
7 Li,PS Simulations details closely follow those of the original work by Gigli al. [80] and the starting configuration was also obtained from the reference. 3 ns NpT simulations at 0 bar were 25 intervals, using PET-MAD MLIP strides and universal FlashMD at 16 fs strides. SVR thermostat [30] was coupled to the system with т = 10 fs, and the BZP barostat [92] was used with 7 1000 fs, including Generalized Langevin Equation [107] (GLE) thermostat coupled to the cell parameters at the same value of r. From the resulting MDD trajectories, MSD of Li was first computed and used to calculate the Li conductivity equation.
performed for temperatures between 575 and 725
fs
with the
29

---

T
[K]
909
770
665
Go000
(0.5 fs)
(a 15)
116 fsi
10a
-1
R -1
10-1
b
10
MD
FlashMD (16 (S)
FlashMD (64 1s)
(1 (5)
10-3
MD
MSD [A4]
0.05
FlashMD 116 (s)
10
1.1
1
0.04
1000/T [K-
0.03
6
3
2
1
1
8
3 2
Layer
Figure 12: Results of case studies conducted for another set of universal FlashMD models trained with the trajectories from PET-OMATPES. (a) Ramachandran plots of the main backbone dihedrals for simulation of solvated alanine dipeptide at 300 K. (b) Mean square displacement (MSD) of the Al (110) surface atoms at 500 K, at different layers from the surface (B indicates the limiting value for the bulk). (c) Li conductivities of L varying T.
L Univerals FlashMD models trained with PET-OMATPES
Here, we demonstrate the performance of another set of universal FlashMD models that adopt different universal MLIP as the baseline. The new model of choice is PET-OMATPES [108], which is an MLIP that has first been trained on the OMat24 dataset [109], then fine-tuned on the MatPES- r2SCAN dataset [1] Using this model, we followed the same protocol detailed App. to generate the reference trajectories that can be used for training, then trained the new set of universal FlashMD models. These models have also been made available the FlashMD HuggingFace repository. The same set of case studies the original universal FlashMD models (see Fig. 3) have been performed. We note that the alanine dipeptide case study was newly conducted at 300 K. as the r2SCAN functional does not suffer from the inaccuracies in describing aqucous systems that PBEsol functional has. The new set results in Fig. 12 show that the PET -EAA-OM---P---EEE based universal FlashMD models successfully capture the expected trends in all three case studies.
M Chaos and Lyapunov exponents
Due to the chaotic nature of the Hamiltonian dynamics the vast majority of non-trivial microscopic systems, small enrors in the prediction of one step nearly always propagate exponentially (in the elapsed time) large errors with respect to the reference trajectory. Due to this phenomenon, calculating
in the reproduction of an individual trajectory is nearly meaningless, and what matters in practice is the correct description of collective dynamical modes. Therefore, instead of comparing the long-time accuracy of FlashMD against that of MD, basic (but still meaningful) experiment is that of comparing the Lyapunov exponents [44] of the respective dynamics.
In the following table, report the Lyapunov exponents of the dynamics of two different chemical systems: a liquid water system and a solid disordered GeTe system, both equilibrated at 300 K with standard MD with Langevin thermostat and then run targeting NVE conditions both with MD and FlashMD, using the newer r SCCAN models. We used 0.5 MD and 16 FlashMD for the water system, and 2 fs MD and 64 fs FlashMD for the GeTe system.
This result shows that FlashMD displays similar chaotic behavior that the baseline Hamiltonian dynamical system, although further evidence would be needed to establish these facts with certainty.
30

---

Table 5: Inverse of the largest Lyapunov exponents in the MD and FlashMD dynamics of water and Ge'Te
System
MD
FlashMD
Water
137 fs
114 fs
GeTe
294
238
N Computational resources
The use of computational resources in this work mainly stems from the generation of the universal dataset of NVE trajectories, which employed 20,000 GPU hours on an Nvidia GH200 cluster. Model training was performed on Nvidia H100 GPUs, for total of around 3,000 GPU hours. All other experiments, mostly molecular dynamics, were run either on H100 or L40S GPUs, and they do not contribute the overall total compute in a significant way. Overall, we estimate our total usage of computational resources as slightly under 25000 H100 GPU hours.
31

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: All claims in the abstract are justified by our analysis and findings. The time steps generally employed in MD simulations range from 0.5 to 2 fs, therefore our models (which generally tested in the 16 fs to 64 fs range) provide longer stride by a factor between one and two orders of magnitude. Our justification for the architecture, generalization arbitrary thermodynamic ensembles and analysis of failure modes can be found in the Theory section of this work. Our experiments and relative discussion, on which we base the last sentence of the abstract, are in the Results section.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results be expected to generalize other settings.
It is fine include aspirational goals as motivation as long is clear that these goals are not attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations of the work performed by the authors?
Answer:
[Yes]
Justification: Although we do not provide standalone Limitations section, we have been manifest about the potential shortcomings of our models, especially regarding their inter- action with the global SVR thermostat (Results section) and their occasional failures in reproducing time-dependent properties (in the Appendix, using q- TIP4P// water models).
Guidelines:
The answer NA means that the paper has limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification. asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low or images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and faimess.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best
32

---

judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community, Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [NA]
Justification: The paper does not include theoretical results.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced in the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but they appear in the supplemental material, the authors are encouraged provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless whether the code and data are provided or not)?
Answer: [Yss]
Justification: All the information needed to reproduce the results in the paper is available in the main text and in the Appendices. Based on these alone, a reader would be able to reproduce all of results (quantitatively, in most cases), reaching the same scientific conclusions as us.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is specific model and empirical evaluation, it may be necessary to either make possible for others to replicate the model with the same dataset, provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how replicate the results, access to a hosted model (e.g., in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make it clear how to reproduce that algorithm.
(b) the contribution is primarily a new model architecture, the paper should describe the architecture clearly and fully.
33

---

(c) If the contribution is a new model (e.g.. a large language model), then there should either way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome to describe the particular way they provide for reproducibility. the case of closed-source models, may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers have some path reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [Yes]
Justification: A repository is made available with all the code and instructions to reproduce all results the paper. This included as Supplementary Material.
Guidelines:
The answer NA means that paper does not include experiments requiring code. Please see the NeurIPS code and data submission guidelines (https://nips. cc/ public/ge cublic/guides//
comPolicy) for more details.
While encourage the release of code and data, understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https: / //nips. cctpcccpppcccccccccccccc weabblccggiidss.... hct////oooSSSssiiieeeeeee sesiisssssiioyyy for more details.
The authors should provide instructions data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: Our training details are fully explained in the corresponding Appendix.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core the paper to a level of detail that is necessary appreciate the results and make sense of them.
The full details can be provided either with the code, appendix, or as supplemental material.
7. Experiment statistical significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [Yes]
34

---

Justification: All results that are affected by significant statistical error are reported together with error bounds. These include the average water temperatures and pressures the Results section, as well as the mean square displacements and dipole autocorrelation functions discussed the Appendix.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper.
The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g.. Normally distributed errors).
should be clear whether the error bar is the standard deviation the standard error of the mean.
is OK to report I-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables figures symmetric error bars that would yield results that are out range (c.g. negative error rates).
If error bars reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification:
This information is reported in the corresponding Appendix.
Guidelines:
The answer NA means that the paper does include experiments.
The paper should indicate type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't 't make it into the paper).
9. Code of ethics
Question: Does the research conducted in the paper conform, in every respect, with the inneeeiiiideeeee
NeurlPS Code of Ethics https: http///eeurrp..eee
Answer:
[Yes]
Justification: All authors have reviewed, adhered to, and respected the Code of Ethics.
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., there a special consid- cration due to laws or regulations their jurisdiction).
10. Broader impacts
35

---

Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [Yes]
Justification: It is difficult to estimate the societal impact of this paper, as its main result is the acceleration of existing modeling techniques empowering practitioners access longer time scales in their computational studies. To the best of our knowledge, these modeling techniques are overwhelmingly used for constructive purposes (drug discovery, modeling of new materials) and so we would expect the outcomes to have an overall positive societal impact. We have expressed this view in the paper, by highlighting that our method would allow practitioners to access longer time scales in their experiments.
Guidelines:
The answer NA means that there is no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious unintended uses (e.g., disinformation. generating fake profiles, surveillance), fairness considerations (e.g.. deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility of MI.).
11. Safeguards
Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Answer: [NA]
Justification:
The paper poses no such risks.
Guidelines:
The answer NA means that the paper poses such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
36

---

Answer: [Yes]
Justification: All the employed software has been used properly within its terms, acknowl- edged and credited. We have not used any other assets.
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (e.g., CC-BY 4.0) should included for each asset.
For scraped data from a particular source (c.g., website), the copyright and terms of service of that source should be provided.
If assets released, the license, copyright information, and terms of in the package should be provided. For popular datasets, paperswithcode. thcode. com/ datasets
has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if has changed) should be provided.
this information is not available online, the authors are encouraged to reach out to the asset's ss creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes]
The paper introduces few general- models for molecular dynamics prediction. These are already accompanied by excellent documentation, which will be made public along with the models.
Justification:
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification:
The paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then much detail as possible should be included the main paper.
According to the NeurIPS Code of Ethics, workers involved in data collection, curation, other labor should be paid at least the minimum wage in the country of the data collector.
Institutional review board (IRB) approvals or equivalent for research with human subjects
37

---

Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: The paper does not involve crowdsoureing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with
human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent)
may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
16. Declaration of LLM usage
Question: Does the paper describe the usage of LLMs if it is important, original, or non-standard component of the core methods this research? Note that the LLM used only for writing, editing, or formatting purposes and does not impact the core methodology, scientific rigorousness, or originality of the research, declaration is not required.
Answer: [NA]
Justification:
We did not involve LLMs in any non-standard or original way.
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, non-standard components.
Please refer to our L.LM policy (https:// http:p//e/rrripp..cco/ooeeeeeeee for what should or should not be described.
38


Reviews:
[{"review_id": "KLwWPDkB9d", "reviewer": "Reviewer_7QyG", "summary": "This paper introduces FlashMD, a machine learning framework designed to accelerate Molecular Dynamics (MD) simulations. Instead of calculating forces at every minuscule time step, FlashMD uses a Graph Neural Network (GNN) to directly predict the positions and momenta of atoms over time strides that are up to two orders of magnitude longer than those used in traditional methods.", "questions": "(1) What is the total simulation time of the result in Figure 3(a) shows? \\\n(2) The experiments typically use hundreds of atoms - however, for practical purposes, it is common to use millions of atoms. What is the scalability of this framework?", "limitations": "yes", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 2, "significance": 3, "originality": 3, "strengths_and_weaknesses": "**Strengths:**\\\n(1) Traditional MD approaches are typically known for their requirement of fine-grained time steps, leading to slow simulations. Approaches toward increasing strides should be encouraged. \\\n(2) The connection between MD and GNN is natural and worth investigating. \\\n**Weaknesses:**\\\n(1) Imprecise statements. \\\n1a) Line 55 is vague and can be misleading. The Hamilton's equations require energy conservation and are improper to directly define MD. Although the authors emphasize \"In its simplest form\" to make it more precise, this can still be misleading for researchers outside the research domain. \n1b) VV approximately conserves energy only because it is applied to an energy-conserved system -- VV can be used without energy conservation, and it is misleading to only emphasize that VV is suitable for NVE due to the ability to conserve energy with sufficiently small time steps. \\\n1c) MD is not necessarily deterministic. e.g., some simulations take Brownian motion into consideration, and some include quantum effects. \\\n(2) Unclear central point. The paper should more focus on the key contribution, which I presume is the orders-of-magnitude increase in time step sizes. However, in the framework section, I only see scattered bullet points and fail to observe their connections. I would recommend focusing on a key point within the 9 pages or submit elsewhere with higher page limits. At current stage, the main text of this paper looks like a list of contents in appendix.", "was_revised": true, "final_justification": "The quality of the results is clarified by the authors, thus raising the score.", "rebuttal": "We thank the Reviewer for their time and effort. We are happy to read that the Reviewer has identified the push towards larger time strides in MD simulations using GNNs as a worthwhile initiative and hence a strength of our work. Regarding the weaknesses that the Reviewer has pointed out, we would like to first make the following clarifications to help resolve their concerns about our work:\n\n1a) The introduction we provide on the subject of molecular dynamics is in agreement with standard textbooks [1, 2] and graduate courses. When stating that our presentation is imprecise, we believe that the Reviewer refers to most practical applications of molecular dynamics, where thermostats are added to the time-integration algorithm, causing fluctuations in the energy which are consistent with the target thermodynamic ensemble. (We are not entirely sure that this is what the Reviewer was referring to and we invite them to inform us through a response if this was not the case.) We would like to emphasize that one of the advantages of our approach, compared to generative methods, is to disentangle Hamiltonian dynamics (which we learn exclusively from NVE data through a GNN) from thermostatting, which can be applied separately and in a very inexpensive fashion. This allows us to generalize a single model to make predictions for arbitrary thermodynamic ensembles and arbitrary thermodynamic state points, which is something that generative models targeting specific distributions (i.e., a specific ensemble at a single thermodynamic state point) cannot achieve.\n\n1b) Here, we believe that the Reviewer is again referring to thermostatted systems, or perhaps systems with external forces, and we would like them to correct us if this is not the case. We clearly state the assumptions that lead to our treatment (e.g., “in the absence of external perturbations”, line 57), and therefore we believe that it is incorrect to characterize our introduction as “imprecise” or “vague”. If the Reviewer is more familiar with other sets of assumptions, we would be happy to make a comparison or to mention them as a limitation of our method. For instance, it is true that velocity Verlet can be applied to non-Hamiltonian systems (or non-separable Hamiltonian systems), and that the symplecticity of the algorithm, and therefore its long-time energy conservation [3], does not apply in those settings.\n\n1c) We refer to 1a) for Langevin dynamics (which is equivalent to Brownian motion and can be seen as a form of thermostatting). Concerning nuclear quantum effects, path integral molecular dynamics can also be reduced to Hamiltonian dynamics and learned with FlashMD. We would be happy to discuss this point more in detail, both with the Reviewer and in the text of the manuscript.\n\n(2) We thank the Reviewer for pointing this out and we agree. Just for completeness, the central contributions are:\n(a) Speeding up molecular dynamics with large time steps; (b) Compared to generative approaches in the same domain, we learn from NVE data to generalize to arbitrary thermodynamic conditions with a single model; (c) Training a universal model (in the chemical sense, i.e., across the whole periodic table) for molecular dynamics predictions. We will definitely clarify and emphasize these points over more technical details in a revised version.\n\nHere are our responses to the Reviewer’s questions:\n\n(1) The total simulation time of the results in Figure 3(a) is 0.5 ns for each of the 10 runs, for all models. We thank the Reviewer for highlighting that this was missing. We will update the simulation details accordingly. We would like to emphasize that we chose these examples because they could be brought close to statistical convergence with a limited effort (approximately one GPU day per trajectory) even when performing the baseline MD simulations with an MLIP. As apparent from the timings in Table 2, the FlashMD models would allow for much longer simulation times with the same effort.\n\n(2) Our framework has the same scalability as current MLIPs. The computational effort scales linearly with system size, although it is limited to a few tens of thousands atoms on a single GPU due to memory reasons. We would be happy to add some scaling plots to the appendices if this Reviewer thinks it would be useful.  Current MLIPs can be applied to large-scale simulations by using multiple (10-100) GPUs in parallel [4], which has allowed machine-learning models to simulate up to around ten million atoms. We are currently working on the integration of our framework with LAMMPS, which involves a few technical hurdles, but no conceptual problem. Once this integration is complete, we expect to achieve a speed on the order of 10 ns per day, or slightly more, on such large systems. However, we would like to emphasize that, while these scales are sometimes needed to perform biochemical simulations, most research on materials using molecular dynamics does not need these scaling efforts and is performed on much smaller systems (between 500 and 5000 atoms).\n\n[1] Frenkel, Daan, and Berend Smit. Understanding molecular simulation: from algorithms to applications. Elsevier, 2023.\n\n[2] Allen, Michael P., and Dominic J. Tildesley. Computer simulation of liquids. Oxford university press, 2017.\n\n[3] Hairer, Ernst, et al. \"Geometric numerical integration.\" Oberwolfach Reports 3.1 (2006): 805-882.\n\n[4] Musaelian, Albert, et al. \"Learning local equivariant representations for large-scale atomistic dynamics.\" Nature Communications 14.1 (2023): 579."}, {"review_id": "P6T7xJW17A", "reviewer": "Reviewer_SKE5", "summary": "This manuscript investigates the problem of learning a “universal” model for directly predicting state changes (coordinates and momenta) as an alternative approach to molecular dynamics (MD) simulations. Unlike previous methods such as MLIPs that focus on energy/force prediction or direct MD propagators on specific systems, the proposed approach aims for broader applicability and faster acceleration. The authors discuss several critical factors for training such a propagator and propose solutions, including the choice of GNN architecture, data augmentation strategies, and energy conservation during both training and inference. In experiments, they demonstrate that the proposed FlashMD model can accurately recover MD ensemble statistics for both single-system cases (e.g., water) and across diverse chemical systems.", "questions": "1. Have the authors attempted to train a single model across multiple temporal scales? What are the potential challenges in learning such multi-timescale dynamics within one model?\n\n1. Can FlashMD be trained on MD datasets that do not follow the NVE ensemble or lack full observations (e.g., missing momenta, solvent state not saved)? How critical are the assumptions of determinism and strict Markovian dynamics in such cases?\n\n1. The authors evaluate energy conservation error for model selection but do not include it as part of the training loss. What are the main challenges in incorporating energy loss directly during training?\n\n1. All experiments are compared only against reference MD from PET-MAD. Why were no machine learning baselines, particularly other “direct MD propagators” discussed in Section 2.3, included for comparison?\n\n1. Scalability and generalization: Current evaluations seem to be limited to small systems (<1000 particles) and short simulations (ps-level). Have the authors tested the model on larger systems or longer timescales?", "limitations": "Yes.", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "## Strengths\n\n1. The manuscript is clearly written and well structured. The main objectives are effectively conveyed, and the authors provide a comprehensive review of the problem, in-depth discussions, and detailed results in the appendix. Overall, the clarity and  quality are high.\n2. The paper offers a solid review of the physical and mathematical foundations of molecular dynamics, demonstrating a strong grounding for the proposed method.\n3. The ablation studies are thorough, covering key factors such as training time steps and enforcement of energy conservation. Showing the importance of considering these key factors.\n4. Results show that a universal, direct MD propagator model can potentially capture the underlying physics across different systems.\n\n## Weaknesses\n\n1. A key limitation is the need to train separate models for different time strides, along with tedious hyperparameter tuning to find an appropriate stride that balances accuracy and speed. Especially for the “universal model,” optimal strides may vary across systems, this limits the model’s “universality” and practical applicability (significance). (Ref. Q1)\n\n1. Some modeling assumptions, though physically sound, may be overly restrictive. For instance, the insistence on determinism and strict Markovian behavior limits the model learn from NVE simulation data with complete observability. In practice, many MD datasets include missing data (e.g., coordinates only on non-solvent atoms) or follow alternative ensembles (e.g., NPT), which violates these assumptions. (Ref. Q2)\n\n1. The model is mainly evaluated against reference MD, with no comparison to other “direct MD propagator models”. This makes it difficult to assess the actual improvement of the proposed approach (Q4). Additionally, evaluations are conducted on small systems with short simulation times (Q5), limiting the scope of validation. (significance)", "was_revised": true, "final_justification": "The authors has addressed my main questions (multi-timescale, training on NVE, energy conservation, prediction strides). Remaining issue is minor and universal (scalability). The manuscript is pleasant to read and insightful. Recommend for acceptance.", "rebuttal": "We thank the Reviewer for their time and their feedback. On our side, it was pleasing to receive recognition for our efforts to formally ground our work in the fundamentals of molecular dynamics, which are often overlooked in similar works. Below are our responses to the Reviewer’s questions:\n\n1. What the Reviewer refers to would constitute a very natural extension of our work, and it is a very relevant question. We have indeed tried to train a single model across multiple temporal scales. We found that the learning exercise becomes very difficult due to the need to learn the larger steps, where potentially many systems in the training set are in a chaotic (and therefore unlearnable, with current techniques) regime. Unfortunately, this negatively impacts the quality of the predictions also on smaller time steps. One additional consideration that made us lean towards training a single model per time step is that of computational efficiency. Learning arbitrary time horizons with a single model might require larger (and therefore slower) models, in a research field where simulation time is often a limiting factor. In this case, we believe that training and delivering one model per time step is more advantageous, especially from the point of view of practitioners. We thank the Reviewer for this question and we would be happy to include a discussion along these lines in a revised version.\n\n2. We would like to begin by emphasizing that, from the point of view of the target ensembles, training on NVE data is not a limitation, but a deliberate choice that makes our model more transferable: we demonstrate the use of FlashMD in NVT and NPT sampling by combining the NVE “core” with completely inexpensive thermostat steps, and with a barostat that can be implemented with a couple MLIP evaluations. In other words, compared to generative methods, our approach disentangles Hamiltonian dynamics (which we learn exclusively from NVE data through a GNN) from, e.g., thermostatting, which can be applied separately and in a very inexpensive fashion. This allows us to generalize a single model to make predictions for arbitrary thermodynamic ensembles and arbitrary thermodynamic state points, which is something that generative models targeting specific distributions (i.e., a specific ensemble at a single thermodynamic state point) cannot achieve. That said, we see no major hurdle including the cell degrees of freedom to train on NPH trajectories to enable constant-enthalpy/constant-pressure runs without the additional MLIP evaluations that are included in our barostat. Removing the requirement of determinism and complete information from the training data should be possible (similar to how coarse-grained models are trained using the same functional form of a MLIP, but using only a subset of the degrees of freedom), but we have not tried it yet. It is certainly an interesting future direction, as it would require including an intrinsic stochastic component in the model architecture.\n\n3. We thank the Reviewer for this very good question. We have tried including an energy conservation term in the loss, but we found no significant improvement in the energy conservation of the resulting trajectories, which is instead dominated by highly non-trivial error accumulation phenomena over tens or thousands of steps. We found that long-time energy conservation can only be enforced reliably by having an exactly symplectic model, and this is in complete agreement with the literature on numerical integrators for Hamiltonian systems [1]. (We experimented with symplectic architectures, but we found them to be too slow for molecular dynamics applications. This is a very promising future direction and we will definitely include a discussion of this aspect in a future version with less strict page limits.) Given that encouraging energy conservation through the loss function has only very limited positive effects on the quality of the resulting dynamics, and that evaluating an MLIP at every training step makes the training procedure nearly twice as expensive, we decided to train without this feature. If the Reviewer thinks this is helpful, we can definitely add a few sentences to explain why we decided not to add an energy conservation term in the loss for our models.\n\n4. As the Reviewer points out, we identified two methods that are directly comparable to ours: MDNet and TrajCast. MDNet was presented in a very short, proof-of-principle, workshop paper [2] which only evaluated on what could be considered a \"toy system\": that of liquid argon with a Lennard-Jones potential. Since the MDNet code is not available (as far as we know), we had to use the instructions in the original paper to compare with the published results. This is what we did in Appendix H. We did not believe that this Lennard-Jones example was significant enough to appear in the main text, especially given the page limit. Regarding TrajCast [3], we did not compare to TrajCast explicitly because we had limited time to do so. TrajCast was released on arxiv six weeks before the submission deadline (incidentally, this means that it should be considered as “contemporaneous” work, according to the NeurIPS guidelines), and we came to know of it only three weeks after that. Nevertheless, we are happy to share some new results that we obtained during the rebuttal period with the Reviewer. We have trained FlashMD on the water trajectory dataset released by the authors of TrajCast (which is generated using the SPC force field as implemented in LAMMPS, and targets a 5 fs prediction step), and compared our results to the results that the TrajCast manuscript reports. FlashMD obtains out-of-the-box (i.e., without dataset-specific hyperparameter optimizations) a momentum MAE of 0.52%, while TrajCast achieves 0.37%. Although the results are comparable and both are accurate enough for stable direct MD propagation, it seems that TrajCast is more accurate on this particular dataset. Note that longer-time-step trajectory datasets were not made available by TrajCast authors, even though we were successfully able to push FlashMD water models for time strides up to 16 fs (i.e., more than three times as much as TrajCast). Still, from indirectly comparing the reported higher-time-step accuracies of TrajCast for 15 fs strides and FlashMD for 16 fs strides (trained on our own water trajectory dataset using the q-TIP4P/f force field), it is likely that our model performs better than TrajCast for longer strides. We would be happy to add these results and considerations in a revised version of our work. Finally, we would like to highlight that a direct MD prediction model is much more useful in practice if it is (1) pre-trained and applicable to the whole periodic table and (2) able to predict long strides to provide significant acceleration. These are two features that our models are designed for and manage to achieve in practice, while TrajCast might be, in our opinion, more limited in its scope.\n\n5. Our framework has roughly the same scalability as current universal MLIPs. The computational effort scales linearly with system size, although it is limited to a few tens of thousands atoms on a single GPU due to memory reasons. We would be happy to add some scaling plots to the appendices if this Reviewer thinks it would be useful.  Current universal MLIPs can be applied to large-scale simulations by using multiple (10-100) GPUs in parallel [4], which has allowed machine-learning models to simulate up to around ten million atoms. We are currently working on the integration of our framework with LAMMPS, which involves a few technical hurdles, but no conceptual problem. Once this integration is complete, we expect to achieve a speed on the order of 10 ns per day, or slightly more, on such large systems. However, we would like to emphasize that, while these scales are sometimes needed to perform biochemical simulations, most research on materials using molecular dynamics does not need these scaling efforts and is performed on much smaller systems (between 500 and 5000 atoms).\n\n[1] Hairer, Ernst, et al. \"Geometric numerical integration.\" Oberwolfach Reports 3.1 (2006): 805-882.\n\n[2] Zheng, Tianze, Weihao Gao, and Chong Wang. \"Learning large-time-step molecular dynamics with graph neural networks.\" arXiv preprint arXiv:2111.15176 (2021).\n\n[3] Thiemann, Fabian L., et al. \"Force-free molecular dynamics through autoregressive equivariant networks.\" arXiv preprint arXiv:2503.23794 (2025).\n\n[4] Musaelian, Albert, et al. \"Learning local equivariant representations for large-scale atomistic dynamics.\" Nature Communications 14.1 (2023): 579."}, {"review_id": "80Qw36RUOg", "reviewer": "Reviewer_5YjL", "summary": "The paper introduces a deep learning model that predicts up to ~100 times larger time steps than MD simulation with machine learning force fields. The model is based on a GNN, which uses positions and velocities of the current state to predict a new configuration ~10-100 MD steps ahead. This potentially allows being faster than MD. The model is trained on a wide variety of systems, and can, hence, be evaluated even on unseen systems, which is demonstrated in the experiments. The investigated evaluation systems are water, alanine dipeptide in water, a metal surface, and a solid-state electrolyte, where the model recovers a phase transition.", "questions": "- There is a lot of work on ML models that predict larger time steps than MD, which are not discussed in this work. Although the here presented approach is slightly different, they should be cited and discussed regardless. Examples include \"Timewarp: Transferable Acceleration of Molecular Dynamics by Learning Time-Coarsened Dynamics\" by Klein et al. and \"Implicit Transfer Operator Learning: Multiple Time-Resolution Surrogates for Molecular Dynamics\" by Schreiner et al. \n- The PET-MAD MLIP is used throughout the paper, I think it would help to have a short section discussing this model\n-  The FlashMD framework section is a little short, a lot of important things, like the training, are only mentioned in the appendix. I am not sure that this can be fixed given the page limit, though.  \n- For the metal surface experiment: From figure 3c it seems that the transition temperature is predicted differently, while the text states that it matches. Why the discrepancy? Maybe there are some error bars that are not shown?\n- I expect that the systems the model was evaluated on were not seen during training for the universal model. Is this correct? And how similar are the systems in the training set?", "limitations": "Yes", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 3, "strengths_and_weaknesses": "Strengths:\n- The paper is well written and clearly motivated. Especially, MD simulations and potential pitfalls are in depth discussed.\n- Although the idea of predicting larger time steps with ML models is not novel, the proposed architecture with the filtering is. Moreover, the model seems to be much more transferable than previous models. \n- They evaluate the same trained model on a variety of systems, highlighting the transferability of the universal model\n- The model can be applied to different ensembles, e.g.. NVE, NVT, NPT.\n\nWeaknesses:\n- For the alanine dipeptide experiments, the positive phi states are not shown.\n- Results are mostly qualitative\n- Wall-clock speed of the method is unclear. Although FlashMD works with larger time steps, it might still be slower than using an MLIP, e.g., PET-MAD. It most likely is slower than classical MD, but might be more accurate.\n- See also questions", "was_revised": true, "final_justification": "The authors addressed my questions and concerns in the rebuttal. Especially important was that the universal model was evaluated on systems unseen during training.", "rebuttal": "We thank the Reviewer for their feedback and their questions. As the Reviewer pointed out, most molecular dynamics or sampling predictors are only valid for a single thermodynamic state point, corresponding to the trajectories they were trained on. By disentangling the stochasticity of thermostatted dynamics from the underlying (and purely deterministic) Hamiltonian dynamics, a single instance of our models is transferable to arbitrary thermodynamic conditions. From the responses of the other Reviewers, it is clear that this aspect was underemphasized in our initial manuscript, and we will work to make it clearer.\n\nHere are our answers to the questions:\n- We thank the Reviewer for this point. Such pre-existing models that seemingly deal with larger time strides adopt a generative and/or time-recurrent approach, and deviate from the deterministic and Markovian nature of MD. As such, we put a lot more emphasis on works closer to our “direct MD prediction” model training setup in Section 2.3 where we highlight previous works. Nevertheless, we inform the reviewer that our manuscript already highlights ITO by Schreiner et al. [1], and we are happy to newly acknowledge Timewarp by Klein et al. [2] and other related works in the camera-ready version.\n\n- The PET-MAD universal MLIP [3] was used in this work to generate a dataset of molecular dynamics trajectories for our universal model (although any other universal MLIP could have been used). The Reviewer is correct that a description of PET-MAD would be beneficial. We will include an overview of the model as an Appendix.\n\n- As the Reviewer correctly pointed out, many details on the FlashMD architecture and its training had to be relegated to the Appendix given the page limit. In case of acceptance, we will include a more detailed description in the camera-ready version, making use of the additional allowed page.\n\n- In the text for Figure 3(c), we state that the conductivities are “reasonably matched” with systematic under/over-estimations, not the predicted transition temperature (Tc). Also, Tc calculation can be sensitive to the simulation setup and the calculation details (e.g. compare our PET-MAD results with that in Mazitov et al. [3]). Despite the slight discrepancy of 25 K between FlashMD and our own reference MD results, we confirm that the Tc of FlashMD simulations falls within the range of reference values for PET-MAD.\n\n- The training set for the universal FlashMD models does not deliberately include any of the examples we considered. As a matter of fact, the 10,000 structures from the MAD dataset [4] that were used for MD simulations to generate the training set does not include any configurations of Li3PS4, elemental Al, or solvated alanine dipeptide, and only contain 1 random cluster, 1 slab, and 3 bulk configurations with the same elemental stoichiometry as water. We further highlight that the exercises presented in Figure 3 are purely extrapolative exercises for FlashMD.\n\nWe would also like to touch on the weaknesses pointed out by the Reviewer, very briefly:\n- We only showed the negative-phi half of the Ramachandran plots as they are generally believed to be the most interesting part for solvated alanine dipeptide (e.g., see Fig. 4 of Morrone et al. [5], Fig. 3 of Koyama et al. [6]). During the rebuttal, however, we still confirmed that the simulations that we have used to generate the original plots show reasonable consistency also in the positive phi range. We would be happy to update our Figure 3a to also include the positive phi range in the camera-ready version of our manuscript.\n\n- Timings of the simulations are reported in Table 2 of Appendix G. FlashMD is consistently faster than MLIPs, although, as the Reviewer correctly points out, slower than MD with (less accurate and non-reactive) classical force fields.\n\n[1] Schreiner, Mathias, Ole Winther, and Simon Olsson. \"Implicit transfer operator learning: Multiple time-resolution models for molecular dynamics.\" Advances in Neural Information Processing Systems 36 (2023): 36449-36462.\n\n[2] Klein, Leon, et al. \"Timewarp: Transferable acceleration of molecular dynamics by learning time-coarsened dynamics.\" Advances in Neural Information Processing Systems 36 (2023): 52863-52883.\n\n[3] Mazitov, Arslan, et al. \"PET-MAD, a universal interatomic potential for advanced materials modeling.\" arXiv preprint arXiv:2503.14118 (2025).\n\n[4] Mazitov, Arslan, et al. \"Massive Atomic Diversity: a compact universal dataset for atomistic machine learning.\" arXiv preprint arXiv:2506.19674 (2025).\n\n[5] Morrone, Joseph A., et al. \"Efficient multiple time scale molecular dynamics: Using colored noise thermostats to stabilize resonances.\" The Journal of chemical physics 134.1 (2011).\n\n[6] Koyama, Yohei M., Tetsuya J. Kobayashi, and Hiroki R. Ueda. \"Perturbation analyses of intermolecular interactions.\" Physical Review E—Statistical, Nonlinear, and Soft Matter Physics 84.2 (2011): 026704."}, {"review_id": "LoNLYWbCjM", "reviewer": "Reviewer_ZniE", "summary": "This paper introduces FlashMD, a graph neural network-based framework for long-stride molecular dynamics (MD) prediction. Instead of integrating atomic trajectories with femtosecond-scale time steps, FlashMD directly predicts the positions and momenta of atoms at large time strides (up to 64 fs), bypassing both force evaluation and numerical integration. \n\nThe model is trained either on system-specific or universal MD trajectories generated using a machine-learned interatomic potential (PET-MAD), and it is designed to handle diverse thermodynamic ensembles (NVE, NVT, NpT). To stabilize predictions and maintain physical plausibility, the authors incorporate energy conservation at inference time, rotational and time-reversal augmentations during training, and uncertainty quantification to manage chaos and out-of-distribution behavior. \n\nFlashMD is evaluated across a wide set of systems—including water, biomolecules, metals, and solid-state electrolytes—and shows significant speed-ups over conventional MD with competitive accuracy.", "questions": "1.\tHow does FlashMD compare to TrajCast or other equivariant autoregressive GNNs on realistic systems (e.g., alanine dipeptide)?\n2.\tHow many steps can be chained in FlashMD before the trajectory diverges significantly from ground truth MD? Is there any cumulative error analysis?", "limitations": "yes", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 4, "significance": 4, "originality": 2, "strengths_and_weaknesses": "Strengths\n- FlashMD bypasses the standard time-integration paradigm by directly learning large time step updates in atomic trajectories—offering acceleration up to 34× in realistic settings.\n- Universal generalization: The model is trained on a chemically diverse dataset (MAD) and evaluated across multiple systems and ensembles without fine-tuning, supporting the claim of universality.\n- The model includes mass-scaled targets, center-of-mass conservation, optional rotational data augmentation, and energy conservation enforcement during inference.\n- Strong empirical scope: Results span small molecules, biomolecules, and materials systems across different thermodynamic settings, with ablation studies and timing benchmarks.\n\nWeaknesses\n- The model is not symplectic by design. Although energy is conserved via momentum rescaling, phase space volume conservation and time reversibility are not guaranteed, raising questions about long-term stability and correctness.\n- Key prior models such as TrajCast or Equivariant Graph Neural Operators are not evaluated head-to-head on realistic benchmarks. Comparisons to MDNet and EGNN are restricted to toy systems (liquid argon).\n- Limited evidence for stability over long time horizons: The paper does not systematically assess long-time trajectory divergence, especially for chaotic systems. There is no quantitative analysis of error accumulation across multiple FlashMD steps.", "was_revised": false, "final_justification": "", "rebuttal": "We thank the Reviewer for their comments, which are all very appropriate and relevant. Similar to the Reviewer, we also consider the lack of symplecticity by design to be the main limitation of this approach, and we hope that we will have more room to discuss it further in the camera-ready version of the manuscript, if it is accepted.\n\nHere are our answers to the questions:\n- We did not compare to TrajCast [1] explicitly because we had limited time to do so. TrajCast was released on arxiv six weeks before the submission deadline (incidentally, this means that it should be considered as “contemporaneous” work, according to the NeurIPS guidelines), and we came to know of it only three weeks after that. Nevertheless, we are happy to share some new results that we obtained during the rebuttal period with the Reviewer. We have trained FlashMD on the water trajectory dataset released by the authors of TrajCast (which is generated using the SPC force field as implemented in LAMMPS, and targets a 5 fs prediction step), and compared our results to the results that the TrajCast manuscript reports. FlashMD obtains out-of-the-box (i.e., without dataset-specific hyperparameter optimizations) a momentum MAE of 0.52%, while TrajCast achieves 0.37%. Although the results are comparable and both are accurate enough for stable direct MD propagation, it seems that TrajCast is more accurate on this particular dataset. Note that longer-time-step trajectory datasets were not made available by TrajCast authors, even though we were successfully able to push FlashMD water models for time strides up to 16 fs (i.e., more than three times as much as TrajCast). Still, from indirectly comparing the reported higher-time-step accuracies of TrajCast for 15 fs strides and FlashMD for 16 fs strides (trained on our own water trajectory dataset using the q-TIP4P/f force field), it is likely that our model performs better than TrajCast for longer strides. We would be happy to add these results and considerations in a revised version of our work. Finally, we would like to highlight that a direct MD prediction model is much more useful in practice if it is (1) pre-trained and applicable to the whole periodic table and (2) able to predict long strides to provide significant acceleration. These are two features that our models are designed for and manage to achieve in practice, while TrajCast might be, in our opinion, more limited in its scope.\n\n- We thank the Reviewer for this question and believe it is important to clarify the scope of the \"long-term stability\" issue that the Reviewer is referring to. In a statistical sense, our benchmarks demonstrate that FlashMD trajectories are very stable (they do not fail catastrophically) and in semi-quantitative agreement with the baseline short-time integrator (meaning that the ensemble they sample, and the average dynamical properties, are close to those of the reference MD). There are in all likelihood starting conditions or materials compositions that would lead to major instabilities (just as it would be the case for the parent MLIP), but we have not found any this far in our experiments, which involve a number of realistic and challenging systems and tens of millions of time steps. In the sense of “how long will it take until the actual trajectory will deviate completely from the specific MD trajectory” the answer is “a few steps, but it does not matter in practice”. MD is chaotic by nature and even with a VV integrator small perturbations of the initial configuration leads to exponential deviations in the trajectories (this can be modeled as $\\Delta q \\cdot \\textrm{exp}(-t/\\tau)$, where the time constant $\\tau$ is often referred to as the Lyapunov exponent). This means that even a minuscule error in the predicted trajectory would lead to exponentially divergence of the trajectories, which however does not imply that one cannot extract accurate time-dependent properties from the simulation. We ran some additional experiments determining the Lyapunov exponents of MD and FlashMD. Both are very similar (for example, around 400 fs on a GeTe system) meaning that FlashMD is about as chaotic as the baseline. We intend to add these experiments in a revision of our work. These results are very interesting and we thank the Reviewer for prompting our investigation.\n\n[1] Thiemann, Fabian L., et al. \"Force-free molecular dynamics through autoregressive equivariant networks.\" arXiv preprint arXiv:2503.23794 (2025)."}]

Author rebuttal:

Meta-review: This paper presents an interesting new approach to enable larger time-steps in molecular simulations. The paper is well written, and clearly presents what task it is solving. The strengths of the approach is that is transferable. While the approaches not compete with generative molecular dynamics emulators in terms of time-steps possible, e.g. implicit transfer operators and timewarp, it provides and interesting general purpose approach to the problem, which will likely have an immediate impact on the molecular simulation field. Simultaneously, the approach can likely inform and inspire other approaches in related areas. All the weaknesses, raised by reviewers were thoroughly addressed addressed during the rebuttal. We hope the authors, include a thorough discussion of other related approaches in the camera ready version of this manuscript. This discussion should, include methods, such as those mentioned above, and more clearly position their own work in relation to these other approaches, which are orthogonal, rather than competing.
Decision: Accept (spotlight)