Forces are not Enough: Benchmark and Critical Evaluation for Machine Learning Force Fields with Molecular Simulations
Authors: Xiang Fu
Zhenghao Wu
Wujie Wang
Tian Xie
Sinan Keten
Rafael Gomez-Bombarelli
Tommi Jaakkola
Venue: Accepted by TMLR ()

Abstract: Molecular dynamics (MD) simulation techniques are widely used for various natural science applications. Increasingly, machine learning (ML) force field (FF) models begin to replace ab-initio simulations by predicting forces directly from atomic structures. Despite significant progress in this area, such techniques are primarily benchmarked by their force/energy prediction errors, even though the practical use case would be to produce realistic MD trajectories. We aim to fill this gap by introducing a novel benchmark suite for ML MD simulation. We curate representative MD systems, including water, organic molecules, peptide, and materials, and design evaluation metrics corresponding to the scientific objectives of respective systems. We benchmark a collection of state-of-the-art (SOTA) ML FF models and illustrate, in particular, how the commonly benchmarked force accuracy is not well aligned with relevant simulation metrics. We demonstrate when and how selected SOTA methods fail, along with offering directions for further improvement. Specifically, we identify stability as a key metric for ML models to improve. Our benchmark suite comes with a comprehensive open-source codebase for training and simulation with ML FFs to facilitate future work.

Published in Transactions on Machine Learning Research (04/2023)
Forces are not Enough: Benchmark and Critical Evaluation for Machine Learning Force Fields with Molecular Simulations
Xiang Fu
ziongfu@mitt.cm edu
Massachusetts Institute of Technology
Zhenghao Wu
zhenghoo. wmL@morthaussttrr..ei
Northwesterr University
Wujie Wang
wuju@mit. edu
Massachusetts Institute of Technology
Tian Xie Microsoft Research
tianzie@@micrrooft.ccmm
Sinan Keten Northwesterr University
wketen&@northursttern.
Rafael Gomez-Bombarelli Massachusetts Institute Technology
rafagb@mit edu
Tommi Jaakkola
tommi@csail. mit edu
Massachusetts Institute of Technology
Reviewed on OpenReview: https: // openreview. form?? id= A8pg0ipukt
Abstract
Molecular dynamics (MD) simulation techniques are widely used for various natural science applications. Increasingly, machine learning (ML) force field (FF) models begin to replace ab-inilio simulations by predicting forces directly from atomic structures. Despite significant progress in this area, such techniques are primarily benchmarked by their prediction errors, even though the practical use would be to produce realistic MD trajectories. We aim to fill this gap by introducing a novel benchmark suite for learned MD simulation. We curate representative MD systems, including water, organic molecules, peptide, and materials, and design evaluation metrics corresponding to the scientific objectives of respective systems. We benchmark a collection slate-of-the-art (SOTA) ML FF models and illustrate, in particular, how the commonly benchmarked force accuracy not well aligned with relevant simulation metrics. We demonstrate when and how selected SOTA methods fail, along with offering directions for further improvement. Specifically, stability as key metric for MT models to improve. Our benchmark suite comes with comprehensive open-source codebase for training and simulation with ML FF's to facilitate future work.
identily
1
Introduction
Molecular Dynamics (MD) simulation is a widely used technique that provides atomistic insights into physical phenomena materials and biological systems (Alder Wainwright, 1939; Rahman, 1964; Schlick, 2010). MD simulations use force fields (FFs) characterize the underlying potential energy surface (PES) of the system and simulate long trajectories based Newton's second law (Frenkel & Smit, 2001). The PES itself

---

Published in Transactions on Machine Learning Research (04/2023)
(a)
(h)
(c)
I+1
Erway
Subbiyyy
Target
owest Farter Rersiio Glenfonc
MAB (4)
Modil
Forry MAE(LL)
Model
Radial Distribution
Diffusion
Dihedral Angles
State
Figure 1: (a) Results on simulating a water system with ML force fields. Models are sorted by force mean absolute error (MAE) in descending order. High stability and low radial distribution function (RDF) MAE are better. Performance in force error does not align with simulation-base may not reveal key factors in simulating MD with ML force ficlds. In this toy example, model 2 (green) has lower force error but likely leads to unstable simulations due to extreme forces from local pathological behavior. (c) Illustrations of MD observables.
metrics. (b) Force-only evaluation
is challenging to compute and would ideally be done through quantum chemistry which is computationally expensive. Traditionally, the altcrnative has been parameterized force fields built from empirically chosen functional forms (Halgren, 1996). Recently, machine learning (ML) force fields (Unke et al., 2021b) have shown promise to accelerate MD simulations by orders of magnitude while being quantum chemically accurate. The learned force field can then simulate MD or perform structural relaxation. However, despite MD simulations being a primary motivation for ML FFs, the evidence supporting the utility of ML FFs often based on their accuracy reconstituting forces across test. cases (Faber et al., 2017) without involving simulations. This paper highlights the importance of simulation-based evaluation, demonstrating that force accuracy alone does not suffice for effective simulation (Figure (a)).
In practice, the exact recovery of the trajectories given the initial conditions is not the ultimate goal of learning MD simulations. Instead, MD simulation quality is better evaluated through macroscopic obscrvables that characterize system properties (Leach & Leach, 2001; Tuckerman, 2010). These observables are designed to be predictive of material properties such as diffusivity in electrolyte materials (Webb et al., 2015), or reveal detailed physical mechanisms, such as the folding kinetics of proteins (Lane et al., 2011; Lindorff- Larsen et al., 2011). Although these observables are critical products of MD simulations, systematic evaluations have not been sufficiently studied in the existing literature. In particular, we show that even small errors in the force field parameters can lead to catastrophic failure in tong-time simulation ML force fields may exhibit pathological behavior such as extreme force/cnergy predictions on states not captured by the training data, causing the simulation to become unstable or explode (toy example illustrated in Figure 1 (b)).
A benchmark for MD simulations requires mindful design over the selection of systems, the simulation protocol, and the evaluation metrics: (1) A benchmark suite should cover diverse and representative systems to reflect the various challenges in different MD applications. (2) Simulations are computationally expensive (see Table 7 in Appendix B). An ideal benchmark must balance the cost of evaluation and the system's ss complexity to obtain meaningful metrics with reasonable time and computation.
(3) Selected systems should be well studied in the simulation domain, and chosen metrics should be based on physical observables that characterize the system's 's important degrees of freedom or geometric features to reflect practical utility. We carcfully curate four representative MD systems and design simulation protocols/ment /etrric considering all the abovementioned desiderata.
Are current state-of-the-art (SOTA) ML FF's capable of simulating a variety of MD systems? What might cause model to fail in simulations? This paper aims to answer these questions with novel benchmark study. The contributions of this paper include:
We introduce a novel benchmark suite for ML MD simulation with simulation protocols and quanti- tative metrics. We perform extensive experiments to benchmark collection of SOTA ML models.
2

---

Published in Transactions on Machine Learning Research (04/2023)
We provide complete codebase for training and simulating MD with ML FFFs to lower the barrier to entry and facilitate future work.
We show that many existing models are inadequate when evaluated on simulation based benchmarks, even when they show accurate force prodiction (as shown in Figure 1).
By performing and analyzing MD simulations, we summarize common failure modes and discuss the causes and potential solutions to motivate future research.
2 Preliminaries
Training. An ML FF aims to learn the potential energy surface E(x) E R as function of atomic coordinates æ E Nx3 (N is the number of atoms), by fitting atom-wise forces F(x) and energies from a training dataset: {mi, F;, F& J Naata, where 2; E RNx33 F E NR3 E E TR. For evaluation, the test force/ enrryyy prediction accuracy is often used as a proxy to quantify the quality of the learned PES. The force field learning protocol has been well established (Unke et al., 2021b).
MD simulation. Simulating molecular behaviors requires integrating Newtonian equation of motion dx//df" = m F(x)) with forces obtained by differentiating D(x): -00 To mimic desired thermodynamic conditions such as constant temperature/pree pressure, an appropriate thermostat and barostat are chosen to augment the equation of motion with extended variables. These conditions are system-dependent and task-dependent. The simulation produces a time series of positions (æz, E NRN333 Ji-o (and velocities), where is the temporal order index, and T the total simulation steps.
Observables. From the time series of positions (and velocities). observables Oxz;) can be computed characterize the state of the system at different gramularities. Under the ergodic hypothesis, the time averages of the simulation observables converge to distributional averages under the Gibbs measure (Frenkel & Smit, 2001): (0) = + lim E. O(z,) f dæ p(æ) O(z), where p(æ) x exp(- E(z)) with T as the bath lemperature and kn as the Boltzmann constant. Calculations of such observa require the system to reach equilibrium. Simulation observables connect simulations to experimental measurements and are predictive of macroscopic propertics of matter. Common observables (illustrated in Figure (c)) include radial distribution functions (RDFs), virial stress lensor, mean-squared displacement (MSD), dihedral angles, etc. We propose evaluation metrics based well-established observables in the respective types of systems (Section 5).
kg
3 Related Work
ML force fields learn the potential energy surface (PES) from the data by applying expressive regressors such as kernel methods (Chmiela et al., 2017) and neural networks on symmetry--reservive of atomic environments (Behler Parrinello, 2007: Khorshidi & Peterson, 2016; Smith et al., 2017; Artrith et al., 2017; Unke & Mcuwly, 2018; Zhang et al., 2018b;a; Kovács et al., 2021; Thölke & De Fabritiis, 2021; Takamoto et al., 2022). Recently, graph neural network architectures (Gilmer et al., 2017; Schütt et al., 2017; Gasteiger et al., 2020; Liu et al., 2021) have gained popularity as they provide a systematic strategy for building many-body correlation functions to capture highly complex PES. In particular, equivariant representations have been shown powerful in representing atomic environments (Satorras et al., 2021; Thomas et al., 2018; Qiao et al., 2021; Schütt et al., 2021; Batzner al., 2022; Gasteiger et al., 2021; Liao & Smidt, 2022; Chen & Ong, 2022; Gasteiger et al., 2022), leading to significant improvements in benchmarks such as MD17 and OC22/20. Some works presented simulation-based results (Unke al., 2021a; Park et al., 2021; Batzner et al., 2022; Musaelian et al., 2022) but do not compare different models with simulation-based metries.
representations
Existing benchmarks for ML force fields (Ramakrishnan et al., 2014; Chmnilaa et al., 2017) mostly focus on force/energy prediction, with small molecules being the most typical systems. The catalyst-focusee OC20 (Chamussot et al., 2021) and OC22 (Tran et al., 2022) benchmarks focus on structural relaxation with force computations, where force prediction is part of the evaluation metrics. The structural relaxation processes are around hundreds of steps, and the goal is to predict the final relaxed structure/energy. These
3

---

Published in Transactions on Machine Learning Research (04/2023)
(a)
(b)
(C)
(c)
Figure 2: Visualization of the benchmarked systems. (a) MD17 molecules: Aspirin, Ethanol, Naphthalene, and Salicylic acid. (b) 64 water molecules. (c) 512 water molecules. (d) Alanine dipeptide. (e) LiPS.
Table 1:
Dataset summary. PBC stands for periodic boundary conditions. uses Metadynamics with implicit solvation.
"Simulation of alanine dipeptide
Dataset
System Type
PBC
#Atoms
Simulation Length
Objectives
MD17 Water Water-large Alanine dipeptide
Small molecule
X
9-21
300 ps (600k steps)
Interatomic distances
liquid
V
192
500 ps (500k steps)
RDF, Diffusivity
liquid
V
1536
150 ps (150k steps)
RDF, Diffusivity
Peptide
X
22
5 ns (2.5M steps)"
Dihedral angle analysis
LiPS
solid-state materials
V 83
50 ps (200k steps)
RDF, Diffusivity
tasks do not characterize system propertics under a structural ensemble, which requires simulations that are millions of steps long. Several recent works (Rosenberger et al., 2021) have also studied the utility of certain ML FFs in MD simulations. In particular, Stocker et al. 2022 uses GemNet (Gasteiger et al., 2021) to simulate small molecules in the QM7-x (Hoja et al., 2021) dataset, with focus on prediction robustness. Zhai et al. 2022 applies the DeepMD (Zhang et al., 2018a) architecture to simulate water and demonstrates its shortcoming in generalizat ion across different phases. However, existing works focus on a single system and model, while evaluation protocols and quantitative metrics for model comparison remain subject to debate. The experimental results in past rescarch are often presented qualitatively or as figures, which can be hard to analyze or used as scalar metric for comparing performance of different models. As previous works are mostly motivated by specific different data generation processes, making it hard to disentangle the impact from the ML model and data. Since the experiments are designed for a single model, the simulation protocols in previous works are often too expensive for large-scale evaluation. Systematic benchmarks for simulation-based metries are lacking in the existing literature, which obscures the challenges applying ML FF for MD applications.
cases, they often lack coverage over diverse MD systems and have very
4 Datasets
Popular benchmark datasets, such as MD17, focus the force prediction task for gas-phase small molecules. However, successes in these tasks are not sufficient evidence for (1) capturing complex interatomic interactions that are critical for condensed phase systems; and (2) recovery of critical simulation observables that cannot be directly indicated by force prediction accuracy. This work focuses on atomic-level MD simulations that manifest complex intermolecular interactions at multiple scales. We choose systems that (1) have been frequently used force field development (Henderson, 1974); (2) cover diverse MD applications such as materials and biology; and (3) can be simulated within reasonable time and compute. The simulation conditions such as temperature and time step sizes (Tuckerman, 2010) are configured to replicate realistic settings in MD literature. Beyond force predictions, we conduct simulations and benchmark observables that reflect the actual simulation quality, along with stability and computational efficiency. The selected systems are summarized in Table 1.
MD17 (Chmicla al., 2017) dataset contains AIMD calculations for eight small organic molecules and is widely used as a force prediction benchmark for ML FFs. We adopt four molecules from MD17 and benchmark the simulation performance. In addition to force error, we evaluate the stability and the distribution of
4

---

Published in Transactions on Machine Learning Research (04/2023)
interatomic distances h(r). For each molecule, we randomly sample 9,500 configurations for training and 500 for validation from the MD17 database. We randomly sample 10,000 configurations from the rest of the data for force error evaluation. We perform five simulations of 300 ps for each model/molecute by initializing from 5 randomly sampled testing configurations, with a time step of 0.5 fs, at 500 K temperature, under a Nosé Hoover thermostat.
Water is arguably the most important molecular fluid in biological and chemical processes. Due to its complex thermodynamic and phase behavior, it poses great challenges for molecular simulations. In addition to force error, we evaluate simulation stability and recovery of both equilibrium statistics and dynamical statistics, namely the clement-conditioned RDFs and liquid diffusion cocfficient. Our datasct consists of 100,000 structures collected every 10 Is from a 1 ns trajectory sampled at equilibrium and a temperature of 300 K. We benchmarked all models with various training validation dataset sizes (1k/10k/900) randomly sampled structures) and used the remaining 10,000 structures for testing. We performed 5 simulations of 500 ps by initializing from 5 randomly sampled testing configurations, with a time step of 1 fs, at 300 K lemperature, under Hoover thermostat. Additionally, we evaluate model generalization to a larger system of 512 water molecules with 5 simulations of 150 ps.
Alanine dipeptide features multiple metastable states, making it a classic benchmark for the development of MD sampling methods and force fields (Head-Gordon et al., 1989; Kaminski et al., 2001). Its geometric flexibility well represented by the central dihedral (torsional) angles & and sb. Our reference data are obtained from simulations with explicit water molecules, with detailed protocols described in Appendix A. For faster simulation, we learn an implicitly solvated FF following a protocol similar to Chen et al. (2021). Our task is more challenging in that it aims to learn the implicitly solvated atomistic FF rather than the implicit solvation correction in Chen et al. (2021). To facilitate accelerated sampling, we apply metadynamics with & and & as the collective variables. We evaluate force prediction, simulation stability, and free energy surface (FES) reconstruction F(6, v). Our dataset consists of 50,000 structures dumped every 2 ps from a 100 ns trajectory at a temperature of 300 K. We used 38,000 randomly sampled structures for training, 2,000 for validation, and the rest as a lest set. We performed 6 simulations of ns by initializing from six local minima on the FES (Figure 5) with a time step of 2 fs at 300 K, and under a Langevin thermostat to mimic random noise from solvation effects. More information on our simulation protocols can be found in Appendix A. All simulations use hydrogen bond length constraints to enable the 2 fs time step.
LiPS is a crystalline superionic lithium conductor relevant to battery development and a representative system for MD simulation usage in studying kinetic properties in materials. We adopt this dataset from Batzner et al. 2022, and benchmark all models on their force error, stability, RDF recovery, and Li-ion diffusivity coefficient. The dataset has 25,000 structures in total, from which we use 19,000 randomly sampled structures for training, 1,000 structures for validation, and the rest for computing force error. We conduct 5 simulations of 50 ps by initializing from 5 randomly sampled testing configurations, with time step of 0.25 fs, at 520 K temperature, under Hoover thermostat.
5 Evaluation Metrics
We deisgn benchmark metrics based on physical observables that describe the most fundamental structural and dynamical propertics of MD system. These observables are relatively easy to understand and compute, and their accurate recovery is often deemed a prerequisite for computing more sophisticated observables. We note there other observables such as thermal conductivity, viscosity, etc. of practical interest. However, these observables are often used for specific applications and require extra domain knowledge or non-standard MD set ups to compute. Our benchmark aims to simplify the MD procedure while being comprehensive, and focus on analyzing the performance of the ML models.
Distribution of interatomic distances is a low-dimensional description of the 3D structure and has been studied in previous work (Zhang et al., 2018a). For a given configuration æ, the distribution of interatomic distances h(()) is computed with:
5

---

Published in Transactions on Machine Learning Research (04/2023)
N Σ Σ j#=i
M()) = N(N - 1)
8(r
22 22yl)
(1)
where r is the distance from a reference particle; is the total number of particles; i,j indicates the pairs of atoms that contribute to the distance statistics; 8 is the Dirac Delta function to extract value distributions. To calculate the ensemble average, h(r) is calculated and averaged over frames from equilibrated trajectories. The MAE is then calculated by integrating r: MAEE()) = -, Kh(r)) - (h(r)) where (-) is the averaging operator, (b(r)) is the reference equilibrium h(r), and (h(r)) is the model preeel-predited equilibrium h(r).
RDF. As one of the most informative simulation observables, the radial distribution function (RDF) describes properties of the system and is experimentally measurable (Yarnell et al., 1973; Filipponi, 1994). It has been widely used in force field development (Henderson, 1974). By definition, the RDF describes how donsity varies as function of distance from particle. For given configuration æ, the RDF can be computed with the following formula:
the /
N Σ +=
RDF()) 4 Np
Σ
8(r - ||= - -xyll)
where r is the distance from a reference particle; N is the total number of particles; í,j indicates the pairs of atoms that contribute to the distance statistics; p is the density of the system; is the Dirac Delta function to extract value distributions. To calculate the ensemble average, RDF(r) is calculated and averaged over frames from equilibrated trajectories. The final RDF MAE is then calculated by integrating r: MAERDF Srza (RDF(r) - (RDF(r)) dr, where (-) is the averaging operator, (RDF(r)) is the reference equilibrium RDF, and (R)F(r)) is the model-predicted RDF.
Diffusivity coefficient is relevant to many practical applications such as battery design (Xie et al., 2022). The diffusivity coefficient D quantifies the time-correlation of the translational displacement, and can be computed from the mean square displacement:
N (t) (0)2
lim +00
(3)
Gt N
where (t) is the coordinate of particle i at time t, and N' is the number of particles being tracked in the system. For the water system, we monitor the liquid diffusivity coefficient and track all 64 oxygen atoms. For the LiPS system, we monitor Li-ion Diffusivity and track all 27 Li-ions. Accurate recovery D requires sufficient long trajectories as it converges with long simulation time. In this paper, we only compute diffusivity for stable trajectories of at least 100 ps for water and 40 ps for LiPS. As we simulate multiple runs for each model, we average the diffusivity coefficient extracted from each valid trajectory to obtain the final prediction.
Free energy surface. Given the probability distributions over configurations p(z) and chosen geometric coordinate& & transformed from x and the marginalized density p(E), the free energy (Tuckerman, 2010) can be calculated from F(E) = T p(E). In the specific case of alanine dipeptide, there are two main conformational degrees of freedom: dihedral angle & of C-N C -C C and dihedral angle vi of N-C, C. C-N. Therefore, the FES w.r.t ó and &o is the most physically informative. We propose our quantitative metric MAEF(e), F(e) based on the absolute error in reconstructing the FES along the ø and & coordinates. We integrate the absolute difference between the reference free energy F and the model predicted F from [-x,r)) .): MAEF(() = 1 FF()) - F()) dó. Error for F(e) is defined similarly.
Quantifying simulation stability. ML FF's can produce unstable dynamics, as the learned force field may not extrapolate robustly to the undersampled configuration space and predict extreme forces. The trajectory can enter nonphysical states that would never happen in a realistic simulation (c.g., bond brcaking should never happen for non-reactive dynamics a low temperature). Predictions can then become increasingly pathological as the input states are far from being captured in the training data. Such unstable trajectories
6

---

Published in Transactions on Machine Learning Research (04/2023)
Table Models benchmarked this work. The paanlation/rotatiinn symmetries are respected by the feature representation at every layer. Number of paramcters on the MD17 datasct are reported.
Model
Symmetry Principle of Geometric Features Energy Conservation
#Parameters
DeepPot-SE (Zhang al., 2018b)
E(3)-invarianl -inviiiari
1.04M
SchNet (Schütt al., 2017)
0.12M
DimeNet (Gasteiger et al., 2020)
E[3) (i3)-ivvriant
2.1M
PaiNN (Schütt ct al., 2021)
SE(3) -equiv-auvar
0.59M
SphereNet (Lin et al., 2021)
E(3)-invariant
<
1.89M
ForceNet (Hu al., 2021)
Pranslation-invaiiant
X
11.37M
GemNet-T (Gasteiger al., 2021)
E(3)-invariant
V
1.89M
GemNet-dT (Gasteiger et al., 2021)
SE(J) -qquivraaaa
X
2.31M
NequlP (Batzner et al., 2022)
1.05M
are not meaningful for observable calculations. For a fair comparison of different models, we need to quantify stability, and compute ensemble statistics only over the stable part of the simulated trajectories.
Although energy drift has been used in previous work in classical MD to monitor stability (Tuckerman et al., 1992), it is not ideal for our benchmark. As data generation using quantum chemistry is expensive, we cannot assume access to the ground truth energy function at evaluation time. Monitoring the energy drift predicted by the ML model is inconsistent across models and unreliable. Since all MD systems studiod in this paper are in equilibrium, we instead keep track of stability by closely monitoring equilibrium statistics. Such equilibrium statistics have range of sensible values that a realistic simulation will never go outside the range. We can characterize the physical range by computing the equilibrium statistics from reference data and assert a simulation becomes "unstable" when the deviation from equilibrium exceeds the realistic range by too far.
Stability criterion. For systems with periodic boundary conditions, we monitor the RDF and say a simulation becomes "unstable" at time T when:
((DDF()) (RDF ( A
(4)
where (+) is the averaging operator, T is a short time window, and A is the stability threshold. In this paper we use == 1 ps, A = 3.0 for water, and == 1 ps, A = 1.0 for LiPS. For water, we assert unstable if any of the three cloment-conditioned RDFs: (D.((,,), RDF (H,H), RDF RDF(H.,) (H.O) exceods the threshold.
For flexible molecules, we keep track of stability through the bond lengths and say a simulation becomes "unstable" time T when:
[T] .T((T)|l - brr)) A (K,))
(5)
where B is the set of all bonds, i.j are the two endpoints of the bond, and big is the equilibrium bond length. For both MD17 molecules and alanine dipeptide, we use A = 0.5 Á.
Our chosen stability thresholds are set to relaxed values to dotect catastrophic failure that the model cannot recover from. We include more details on threshold selection in Appendix B. All metrics related to observable prodiction are computed only over the stable part of the trajectories decouple accuracy and stability.
6 Experiments
Benchmarked models. In our experiments, we aim to cover diverse models with different design such as architecture and symmetry principles. We adopt the Open Catalyst Project implementation of SchNet (Schütt et al., 2017), DimoNet (Gasteiger et al., 2020), ForceNet (Hu et al., 2021), PaiNN (Schütt al., 2021), CemNet-T/dT (Gasteiger et al., 2021), and the official implementation of DeepPot-SE. (Zhang al., 2018b), SphereNet (Liu et al., 2021), and NequlP (Batzner et al., 2022). A summary of all benchmarked models is in
7

---

Published in Transactions on Machine Learning Research (04/2023)
Table 3: Results on MD17. Darker green color indicates better performance. For all results, force MAE is reported in the unit of [meV/A], and stability is reported in the unit of [ps]. The distribution of interatomic distances h(r) MAE is unitless. FPS stands for frames per second. For all metrics (4) indicates the lower the better, and (f) indicates the higher the better. Standard deviation from 5 simulations is in subscript for
applicable metrics.
Model
DesPot-SE
SchNett
DimeNet
PaiNN
SphereNet
ForceNet
GemNet-T
GemNet-dTT
NequlP
OPIS
Molecule
Aspirin
Foree (1)
21.0
35.6
10.0
9.2
3.4
22.1
3.3
5.1
2.3
Stabiiity (1)
9(1
64gm
(10)
141
182
T2:50)
192
300
hir) (3)
0.6510.47)
0.36 on
0.04
0.04
0.03g.00
0.56gc 15)
0.04,
0.04g0.0.
0.02:a.c::
0.21
FPS (t)
88.0
108.0
20.8
86.8
17.5
137.3
28.2
56.8
8.4
Ethunol
Foree
Stabiiit
(00
247 (10s)
26 noy
(6)1)0))
33
300
(00
300gg)
300,ey
h(r)
0.0910.00)
0.21
0.15
0.15 ne)
0.13
0.86 05)
D. 1.100000)
0.09
0.088pcc
0.36
PPS
101.0
112.6
21.4
87.3
30.5
141.1
27.1
54.3
Naphthalene
Force
13.4
22,5
5.7
3.8
1.5
9.9
1.5
1.9
1.1
Stability
2006
18(2)
85
300lgo)
Gtz)
3000
S:21
25,30)
300
h(r)
g0.00)
0.09
0.100
0.13
0,14
(0.00)
B. 1.
0.12p.00
1..1
0.30
FPS
109.3 3
110.9
19.1
92.8
18.3
130.2
27.7
53.5
8.2
Salicylic Acid
Force
14.9
2.6
Stability
300gg)
300g;
281 281(
36/20)
1
26
94,a09)
300
MCr)
0.00y0.00)
0.033
(.06
0.033 20)
0.06go.02)
0.35 03)
(.005
(.07g0.03)
0.03gocco
0.25
FPS
94.6
19.4
90.5
143.2
28.5
52.4
Salicyti
(b)
Ethano
Naghthelsne
(d)
Aspiin
Dialiiy
SANDIY (T)
gentee bae cooeeee
NO L))
Noshel Den seattea
MAT
pecrutten Dest Bonatiee seeiit Screectee
Parts MAE (4)
(1)
Fora
Poptitite da Heroll
CTT
Forms MAG (4)
Force MAKE (4)
(4)
Figure 3: Head-to-head comparison of force MAE vs. Stability and Mr) MAE on MD17 molecules. Models are on the x-axis and are sorted according to force error in descending order. High stability and low h(r)
MAE mean better performance. Error bars indicate 95% confidence intervals.
Table 2. These models have been popular in previous benchmark studies for force/energy prediction. They use different model architecture (feed-forward neural networks message-passing neural networks), different representations for atomistic interactions (e.g., use angular/torsione information or not) and respect different levels of euclidean symmetry. We follow all original hyperparameters introduced in the respective papers and only make minimal adjustments when the training is unstable. More details on the hyperparameters can be found in Appendix B.
Key obsorvations. We make two key observations as evidenced in the experimental results:
1. Despite being widely used, force prediction is not sufficient for evaluating ML. FFs. It generally does not align with simulation stability and performance in estimating ensemble properties.
2. While often neglected, stability is crucial prerequisite for practical usage and major bottleneck for ML FFs. Lower force errors do not necessarily give rise to more stable simulations, suggesting stability is fundamental consideration for comparison and model design.
We next go through the experimental results of all four datasets in detail to demonstrate the key observations while making other observations.
MD17. As shown in Table 3, more recent models that lic on the right side of the table generally achieve a lower force error, but may lack stability. Figure 3 selects and rearranges results from Table 3 demonstrate the non-aligned trends of force prediction performance simulation performance, which supports key
8

---

Published in Transactions on Machine Learning Research (04/2023)
Table 4: Results on Water-1Ok. RDF MAE is unit-less. Diffusivity is computed by averaging 5 runs from random initial configurations and its MAE is reported in the unit of [10-9 m² /s]. The reference diffusivity coelficient is 2.3 x 10 m²/s.
DeepPot-SR
SchNet
DimeNet
PaiNN
SphereNet
ForceNet
T GemNet-dT GemNet-T
NequIP
Force (4)
5.8
9.5
1.4
5.1
16.1
10.9
0.7
1.3
1.5
Stability (1)
p447
(00
350(15a)
500
9r0;
328
500
500go)
RDF(0ooo (1)
(.07
0.58ço.02)
(.07
0.12
0.90
(.68)).06)
0.09ga.03)
0.07:o.0a
0.06
RDF (L)
0.06go.02)
0.26ro.01)
0.05
(..7
0.40
0.47
0.06
0.0 0.05ço.01)
0.05
RDF (D.O) (.)
0.191o.05)
0.58
(.23
0.38
1.15
1.29
0.27 12)
0.30
0.27
Diffusivity (L)
0.04
1.78
0.23
0.87
2.24
0.36
0.71
0.18
FPS (t)
91.0
78.9
17.9
64.2
3.1
67.6
11.3
33.7
3.9
500
Diffusivity MAE [10 1m
2.0
0.8
400
RDF MAE
Stability [ps]
1.5
300
200
0.4
1.0
100
0.2
0.5
0.0
10
20
30
10
30
10
15
Force MAE [mev/A]
Force MAE [meV/A]
Force MAE [meV/A]
Model
SchNet DimeNet
PaiNN SphereNet
ForceNet GemNet-T
GemNet-dT NequiP
DeepPot-SE
Figure 4: Comparison of force MAE vs. stability (Left), force MAE vs. RDF MAE (Middle), and force MAE vs. Diffusivity MAE (Right) on the water benchmark. Each model trained with three dataset sizes. The color of point indicates the model identity, while the point size indicates the training dataset size (small: Ik, medium: 10k, large: 90k). Metrics infeasible to extract from certain doael/dataset size (e.g., Diffusivity for unstable models) are not included.
observation 1. SphereNet and GemNet-T/dT can attain a very low force error for all four molecules, but often collapse before the simulation finishes. This observation constitutes key observation 2. We note that although the stable portion of simulated trajectories produced by SphereNet and GemNet-T/diT can recover the h(r) relatively accurately, stability will become bigger issue when the statistics of interest require long simulations, demonstrated in other experiments. On the other hand, despite having relatively high force error, DeepPot-SE. performs very well on simulation-based metrics on all molecules except for Aspirin (Figure 3). With the highest molecular weight, Aspirin is indeed the hardest task in MD17 in the sense that all models attain high force prediction errors on it. PaiNN also attains competitive simulation performance, while its force error is not among the best. As reference, we also include the h(r) error using a popular classical force field: the optimized potentials for liquid simulations (OPLS, Jorgensen & Tirado-Rives 2005) in Table More details on the classical force field results in Appendix A and Figure 7).
We further observe that good stability alone does not imply accurate recovery of trajectory statistics. Although ForceNet remains stable for Ethanol and Naphthalene, the extracted hrr) deviates a lot from the reference (Table 3), indicating that ForceNet does not learn the underlying PES correctly, possibly due to its lack of energy conservation and rotational equivariance. Overall, NequlP is the rest-performing model on MD17. It achieves the best performance both force prediction and simulation-based metrics for all molecules while requiring the highest computational cost. More detailed results on MD17 inchuding study on stability's relation with training epochs and individual h(r) are included in Appendix B.
9

---

Published in Transactions on Machine Learning Research (04/2023)
Table 5: Results on alanine dipeptide. #Finished is the number of simulations stable for 5 ns. MAF of F(6)
and F(2) are reported in the unit of [kJ/mol].
DeepPot-SE
SchNet
DiteNet
PaiNN
SphereNet
ForceNet
GemNet-TT
GemNet-dTT
NequlP
Force (4)
272.1
217.0
239.0
266.2
256.3
284.7
233.5
219.7
213.6
#Finished (f)
0/6
0/6
0/6
0/6
0/6
0/6
6/6
0/6
5/6
Stability (1)
((o)
0g0)
0(0)
0(0)
804(788)
O1o)
5000
O(o)
(189
F(Ø) (1)
104g2)
104(14)
F(2) (4)
110
113
FPS (1)
42.4
12.1
42.2
15.0
36.5
8.3
Table 6: Results on LiPS. Li-ion Diffusivity coefficient computed by averaging runs from random initial configurations. The reference Li-ion diffusivity coefficient 1.35 10 m²/s
DeepPot-SE
SchNet
DimeNet
PaiNN
SphereNet
ForceNet
CemNet-T
ComNet-dTT
NequlP
Force (4)
40.5
28.8
3.2
11.7
8.3
12.8
1.3
1.4
3.7
Stability (t)
4:3)
50/o;
45(10)
50g))
50(m)
14(3)
50(0)
50
50
RDF (4)
(.27
0.05ço.01)
0.05
0.05
0.05
0.42
0.05
0.05
0.04...01)
Diffusivity (L)
0.37
0.57
0.67
0.38
0.35
0.34
FPS (t)
66.1
35.2
14.8
75.7
18.1
72.1
16.9
43.5
8.2
Water. Under a condensed phase system, key observation 1 and 2 are still evident according to Table 4: GemNet-T/dT and DimeNet are the top-3 models in terms of force prediction, but DimeNet lacks stability. The water diffusivity cocfficient requires long (100 ps in our experiments) trajectories to estimate and thus cannol. be extracted for unstable models. DeepPot-SE does not achieve the best force prediction performance but demonstrates decent stability and highly accurate recovery of simulation statistics. Interestingly, SchNet and SphereNet have high force errors but are highly stable. However, the properties are not accurately recovered.
Figure 4 further compares model performance with different training dataset sizes. Key observations and 2 are clearly shown: Models located on the left of each scatter plot have very low force error but may have poor stability. Error in simulation statistics demonstrates a positive correlation with force error. We also observe that more training data can often improve force prediction performance, well simulation stability and accuracy. We observe clear improvements in simulation stability for DcepPot-SE, SchNet, PaiNN, and GemNel-dT when more training data is used. NequlP is the best-performing model, achieving very low force error, excellent stability, and accurate recovery of ensemble statistics, even under the lowest data budget of 1,000 +
structurcs. However, when the training dataset is sufficiently large (90k), many other models have equally good results as NequIP while being more efficient dalaset size also influences the model of choice for a certain dataset. The detailed results for water-1k/901 and the impact of dataset. size are included in Appendix B. Furthermore, we also conduct ablation studies on model sizes and training epochs in Appendix B.
Alanine dipeptide poses unique challenges in sampling different metastable states separated by high free energy barriers and modeling the implicit solvent. Table 5 shows all models have high force errors due to the random forces introduced by the lack of explicit account water molecules. Although the force errors are in the same order of magnitude, all models except GemNet-T and NequTP fail to simulate stably. The FES roconstruction task requires stable simulation for the entire 5 ns. GemNet-T and NequlP are the only two models that manage to finish at least one simulation but produce inaccurate statistics. All other models are not stable enough to produce meaningful results. As the key difference between GemNet-T and GemNet-dTT is energy conservation, this result suggests a strong stabilization effect of enforcing energy conservation through the prediction of scalar energy. We further analyze the results of this task in Section 7.
LIPS. Compared to flexible molecules and liquid water, this solid material system features slower kinetics. From Table 6 we observe that most models are capable of finishing 50-ns simulations stably. In this datasct, the performance on diffusivity estimation and force prediction align well. We observe that both GemNet-T and Te show excellent. force prediction, stability, and recovery of observables, while Gem-et-
10

---

Published in Transactions on Machine Learning Research (04/2023)
(a)
Reference
Delmmo
NequIP
F Dulmmo]
GemNel-T
F Bellng
190
60
40
30
(oograa)
00
30
Fiulmmll
10
-90
NaquiP
GomNet
180
90
90
180
180
90
80
180
[degree]
(degree]
[degree]
[degree]
[degree]
(b)
Referance
[cttmo]
GamNet-T
Referepioe
90
GomNict
[ocgroc
90
20
20
10
10
100
-30
0
100
-100
-00
0
90
160
-180
-00
50
10U
[degree]
[degree]
$ [degreee
[degree]
[degrea]
(d)
GemNet-T, Water
160
Training Data
-346
LiPS, ensemble
1o
Predicted Erotal [eV]
GemNeetT GemNettdT
Stability [pS]
400
-348
[degree]
300
-350
100
352
-18 80
-00
160
10"
0
20
30
40
50
0.5 Time Step [fs]
0.25
[dogree]
Simu ation Time ([s)
Figure 5: (a, b) For each row, the first three panels show Ramachandran plots of the alanine dipeptide FES reconstructed from 5-ns reference vs. 5-ns NequlP /GemNet-T simulation, all using MetaDynamics. The last two panels show F(Ø) and (()) of alanine dipeptide extracted from reference simulation vs. from NequIP simulation. The two rows result from different initialization, annotated as yellow star. (c) (ø, v) distribution of the alanine dipeptide training dataset. The six initialization points marked with stars. NequlP fails to total
remain stable when the simulation starts from the point marked with black color. (d) Model-preditee energy a function of simulation time when simulating the LiPS system using the NVE ensemble. (c) On water-10k, stability does not improve when the time step is reduced for GemNet-T.
is 2.6 times faster. The better efficiency comes from the direct prediction of atomic forces instead of taking the derivative F= = =dE/öa, which also makes GemNet-dT not energy-conserving - a potential issuc we further discuss Section
Implications on model architecture. More recent models utilizing SEE(3) -())-vvaraauu representations and operations such GemNet-dTT and NequlP more expressive and capture interatomic interactions more accurately. This is reflected by their very low force error and accurate recovery of ensemble statistics when not bottlenecked by stability. Moreover, excellent accuracy and stability can be simultaneously achieved. The stability may come from parity-cquivariance and the use of tensor products (Thomas et al., 2018) in manipulating higher-order geometric features. We believe further investigations into the extrapolation behavior induced by different equivariant geometric representations and operations (Batatia et al., 2022) is a fruitful direction in designing more powerful ML FF'.
7 Failure Modes: Causes and Future Directions
A case study on alanine dipeptide simulation. Figure (a) demonstrates how NequlP and GemNet-T fail reconstruct the FES: the ML potentials do not manage to sample much the transition regions and the configuration space with & E [0, 180%] when initialized from configuration with & -655. Figure 5 (b) shows
11

---

Published in Transactions on Machine Learning Research (04/2023)
(a)
NequlP, Alanine dipeptide
Max (ev/A]
10"
10'
10'
10
3622
t= 3702 fs
t= 3704
I= 3706 fs
t= 3716 fs
3550
3600
3650
Simulation Time [fs]
(b)
GemNet-T, T, Water
Max Force (eV/A]
10
10"
220
225
230
235
240
t 220.1 ps
224.2 fs
t 227 ps
t= 228.0 ps
235.0 ps
Simulation Time [ps]
Figure 6: Examples of simulation collapse when applying (a) NequlP to alanine dipeptide and (b) GemNet-T to water. The y-axis shows the maximum force observed on any atom in the system at a certain time step. An orange cross indicates visualized time steps. Notable nonphysical regions are circled. The collapse usually happens within a very short period of time after the initial local errors.
NequlP fails to simulate stably when initialized from the low-density meta-stable state. The reconstructed FES significantly deviates from the reference in both cases. This failure can be partially explained by Figure 5 (c), the training data distribution produced by the reference potential. The relatively high-energy (low-density) regions and transition regions are exactly those that are not sampled enough by the learned force fields. Even though our MD Irajectory is well-equilibrated, the difference in populations ions of different meta-stable states creates data imbalance, making it more challenging for the model to learn well for higher-energy configurations where density relatively low. The lack of energy labels for force matching in the implicit solvent setting exacerbated this problem and poses challenges learning the relative free energy difference between the meta-stable states. This issue implics that generalization across different regions in the conformational space is important challenge for ML FFs. The rest of the FES results included in Figure 13.
Energy conservation. Models that directly predict forces may not conserve energy. Figure 5 (d) demon- strates the evolution of model-predicte total energy for selected models on LiPS, in microcanonical (NVE) ensemble. The energy of an isolated system in the NVE ensemble is, in principle, conserved. We observe that GemNet-T conserves energy, whereas GemNee-dT fails to conserve the predicted total energy. The existence of non-conservative forces breaks the time-reversal symmetry and, therefore, may not be able to reach equilibrium for observable calculation. Energy conscrvation as inductive bias can also affect simulation stability. For alanine dipeptide, observe GemNet-T can achieve stable simulations while GemNet-dT is unstable. On the other hand, GemNet-dTT performs well on the LiPS dataset when coupled with a thermostat. Previous works (Kolluru et al., 2022) also found that energy conservation is not required for SOTA performance on OC20. The usability of non-conservative FFs in simulations requires further careful investigations.
Simulation instability major bottleneck for highly accurate models fail on several simulation tasks. We further experiment with smaller simulation time steps for GemNet-T water (Figure (e)), but stability does not necessarily improve. We observed more training data improves model stability Stocker et al. (2022) also demonstrates that the stability of GemNet improves with larger training sets on QM7-x, which includes high-energy off-equilibrium geometries obtained from normal mode sampling. Normal mode sampling was also used to generate large datasets in earlier ML force fields such as ANNI- (Smith et al., 2017), and MD simulation results are demonstrated. These results indicate that dataset's coverage over the conformational space has crucial implications on the behavior of the final model, regardless of the underlying neural network model. We also observe that the simulation can collapse within short time window after long period of
12

---

Published in Transactions on Machine Learning Research (04/2023)
stable simulation, visualized in Figure 6. In both cases, we observe that the nonphysical configurations first emerge at local regions, which cascade to the entire system very quickly as extremely large forces are being predicted and subsequently integrated into the dynamics. At the end of the visualization, the bonds in the alanine dipeptide system are broken. Therefore, the local-dcseriptor-bass based NequlP model prodicts very small forces. For the water system, the particles are packed in a finite periodic box. The uonphysical configurations exhibit incorrect coordination structures and extremely large forces. To prevent ML. FFs from sampling nonphysical regions, which is common precursor to failed simulation (Figure 6), one can deliberately include distorted and off-equilibrimm geometries in the training data to improve model robustness (Stooce (Stocker et al., 2022). Alternatively, one can resort to active learning (Wang et al., 2020b; Vandermause et al., 2020; Schwalbe-Kode et al., 2021) to acquire new data points based on model uncertainty. Past works also found adding noise to data paired with a denoising objective during training helpful in improving out-of-distributiin test performance on OC20 (Godwin et al., 2021), and in stabilizing learned simulations (Sanchez-Gonzalez et al., 2020). Another relevant line of work in coarsc-graince MD simulation has studied regularization with an empirical "prior energy" (Wang al., 2019b) and post-prediction refinement (Fu et al., 2022) to battle simulation instability.
8 Conclusion and Outlook
We have introduced a diverse suite of MD simulation tasks and conducted a thorough comparison of SOTA ML FF's to reveal novel insights into ML for MD simulation. As shown in our experiments, benchmarking only force error is not sufficient, and simulation-based based metrics should be used to reflect the practical utility of model. We demonstrate case studies on the failure of existing training schemes/mone /oodels to better understand their limitations. The performance of a model can be highly case-dependent. For more challenging MD systems, more expressive atomistic representations may be required. For example, recent work has explored non-local descriptors (Kabylda et al., 2022) aiming at capturing long-range interactions in large molccules. Strictly local equivariant representations (Musaclian et al., 2022) are studied for large systems where computational scalability critical. New datasets (Eastman et al., 2022) and benchmarks have been playing an important role in inspiring future work. Learning coarse-grainen MD (Wang et al., 2019b; Wang & Bombarlom,arllli, 2019; Fu al., 2022) is another avenue to accelerate MD at larger length/ mugth/time time scales.
The possibility of ML in advancing MD simulation is not limited to ML force fields. Enhanced sampling methods enable fast sampling of rare events and have been augmented with MT. techniques (Schneider et al., 2017; Sultan et al., 2018; Holdijk et al., 2022). Differentiable simulations (Schoenholz & Cubuk, 2020; Wang et al., 2020a; Doerr et al., 2021; Ingraham et al., 2018; Greener & Jones, 2021) offer a principled way of learning the force field by directly training the simulation process to reproduce experimental observables (Wang et al., 2020a; 2022; Thaler & Zavadlav, 2021). We hope our datasets and benchmarks will encourage future developments in all related aspects to push the frontier in ML for MD simulations.
Author Contributions
X.F. led the project. X.F., Z.W., W.W., and T.X. designed the benchmark and evaluation protocols; Z.W. performed MD simulations to generate datasets; X.F implemented the models and performed the experiments; X.F.Z.W., W.W.WW. T.X., and T.J. analyzed the results; X.F., Z.W., W.W., T.X., S.K., R.G.B., and T.J. wrote the manuscript.
Acknowledgments
We thank Adrian Roitberg and the Roitberg Research Group for their valuable suggestions rectifying our alanine dipeptide evaluation protocol. We thank Adeesh Kolluru, Simon Batzner, and Alby Musaclian for helping with reproducing previous results. X.F. is supported by the MIT-GIST collaboration. W.W. is supported by Toyota Research Institute. Z.W. is supported by the Center for Hierarchical Materials Design at Northwestern University and Army Rescarch Office. The authors also acknowledge the MIT SuperCloud and the Lincoln Laboratory Supercomputing Center for providing HPC resources. We also acknowledge Torch and PoTorc Geometric.
13

---

Published in Transactions on Machine Learning Research (04/2023)
References
Mark James Abraham, Teemu Murtola, Roland Schulz, Szilárd Páll, Jeremy C. Smith, Berk Hess, and Erik Lindahl. GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers. SoftwareX, 1-2:19-25 September 2015. ISSN 23527110. doi:
URL
elsevier http://www....eeeaaa/a/a////////
Berni Alder and Thomas Everctt Wainwright. Studies in molecular dynamics. general method. The Journal of Chemical Physics, 31(2):459-466 1959.
Nongmuch Artrith, Alexander Urban, and Gerbrand Ceder. Efficient and accurate machine-learning in- terpolation of atomic energies in compositions with many species. Physical Review B, 96(1 (6)011214,,,, 2017.
Tyes Batatia, Simon Batzner, Dávid Péter Kovács, Albert Musaclian, Gregor NC Simm, Ralf Drautz, Christoph Ortner, Boris Kozinsky, and Gábor Csányi. The desigu space of e (3)-equivariant atom-centered interatomic potentials. arXiv preprint arXiu:220500644 066443, 2022.
Simon Batzner, Albert Musaelian, Lixin Sun, Mario Geiger, Jonathan P Mailoa, Mordechai Kornbluth, Nicola Molinari, Tess E Smidt, and Boris Kozinsky. E (3)-oquivariant graph neural networks for data-efficient and accurate interatomic potentials. Nature communications,, 13(1):1-11, 2022.
Jörg Behler and Michele Parrinello. Generalized neural-network representation high-dimensional potential- energy surfaces. Physical review letters, 98(14):144400) 2007.
Lowik Chamussot, Abhishek Das, Siddharth Goyal, Thibaut Tavril, Muhammed Shuaibi, Morgane Riviere, Kevin Tran, Javier Heras-Domingo, Caleb Ho, Weihua Hu, et al. Open catalyst 2020 (oc20) datasct and community challenges. ACS Calabysis, 11(10):6059 6072, 2021.
Chi Chen and Shyue Ping Ong- A universal graph deep learning interatomic potential for the periodic table. Nature Computational Science, 2(11):718-77) 728,7728885 2022.
Yaoyi Chen, Andreas Krämer, Nicholas E Charron, Brooke Husic, Cecilia Clementi, and Frank Noé. Machine learning implicit solvation for molecular dynamics. The Journal of Chemical Physics, 155(8): 084101, 2021.
Stefan Chmicla, Alexandre Tkatchenko, Huziel E Sauceda, Igor Poltavsky, Kristof T Schütt, and Claus-Robert Müller. Machine learning of accurate energy-conservine molecular force fields. Science advances, 3(5): c1603015, 2017.
Danicl A. Colón-Ramos, Patrick La Riviere, Hari Shroff, and Rudolf Oldenbourg. Transforming the development and dissemination of cutting-edge microscopy and computation. Nal Methods, 16(8): 667-669, August 2019. ISSN 1548-7091, 1548-7105. doi: 010////441122--1-0777550 URL http: wwwure. com
014-4----
Wendy D Cornell, Piotr Cieplak, Christopher Bayly, lan R Gould, Kenneth M Merz, David M Ferguson, David C Spellmeyer, Thomas Fox, James W Caldwell, and Peter A Kollman. second generation force field for the simulation of proteins, nucleic acids, and organic molecules j. am. chem. soc. 1995, 117, 5179- 5197. Journal of the American Chemical Society, 118(9):2309 2309, 1996.
Stefan Doerr, Macicj Majewski, Adrià Pércz, Andreas Kramer, Cecilia Clementi, Frank Noc, Toni Giorgino, and Gianni De Fabritiis. Torchnd: A deep learning framework for molecular sinulations. Journal of chemical theory and computation, 27(5):2355-2 2263, 2021.
Peter Eastman, Jason Swails, John D Chodera, Robert T McGibbon, Yutong Zhao, Kyle Beauchamp, Loo-Ping Wang, Andrew Simmonett, Matthew P Harrigan, Chaya D Stern, et al. Openmm Rapid development of high performance algorithms for molecular dynamics. PLoS computational biology, 13(7):
e1005659, 2017.
14

---

Published in Transactions on Machine Learning Research (04/2023)
Peter Eastman, Pavan Kumar Behara, David L. Dotson, Raimondas Galvelis, John E. Herr, Josh T. Horton, Yuezhi Mao, John D. Chodera, Benjamin P. Pritchard, Yuanqing Wang, Gianni De Fabritiis, and Thomas E. Markland. Spice, dataset of drug-like molecules and peptides for training machine learning potentials, 2022.
Felix A Faber, Luke Hutchison, Bing Huang, Justin Gilmer, Samuel S Schoenholz, George Dahl, Oriol Vinyals, Steven Kearnes, Patrick F Riley, and o Anatole Von Lilienfeld. Prediction errors of molecular machine learning models lower than hybrid dft error. Journal of chemical theory and computation, 13(11): 5255 5264, 2017.
Michael Feig. Kinetics from implicit solvent simulations of biomolecules as a function of viscosity. Journal of chemical theory and compulation, 3(5):1734 115:11734-1748 2007.
Adriano Filipponi. The radial distribution function probed by x-ray absorption spectroscopy. Journal of Physies: Condensed Matter, 6(41):8415, 1994.
Daan Frenkel and Berend Smit. Understanding molecular simulation: from algorithms to applications, volume 1. Elsevier, 2001.
Xiang Fu, Tian Xie, Nathan J Rebello, Bradley D Olsen, and Tommi Jaakkola. Simulate lime-ineegrated coarso-grained molecular dynamics with geometric machine learning. arxiv preprint arXiv:2204. 10348, 2022.
Johannes Gasteiger, Janek GroB, and Stephan Günnemann. Directional message passing for molecular graphs. In International Conference Learning Representations, 2020.
Johannes Gasteiger, Florian Becker, and Stephan Günnemann. Gemmet: Universal directional graph neural networks for molecules. Advances in Newrol Information Processing Systems, 34:6790 1466990-602,, 2021.
Johannes Gasteiger, Muhammed Shuaibi, Anuroop Sriram, Stephan Günnemann, Zachary Ward Ulissi, C. Lawrence Zitnick, and Abhishek Das. Gemnet-OC: Developing graph neural networks for large and diverse molecular simulation datasets. Transactions on Machine Learning Research, 2022. URL
https://openreviev..u http:://oppaueeeeeeeeee/////tt/d
Justin Gilmer, Samuel S Schoenholz, Patrick F Riley, Oriol Vinyals, and George E. Dahl. Neural message passing for quantum chemistry. In International conference on machine learning, pp. 1263-1272. PMLR, 2017.
Jonathan Godwin, Michael Schaarschmidl, Alexander I, Gaunt, Alvaro Gonnaez, Yulia Rubanova, Petar Veličković, James Kirkpatrick, and Peter Battaglia. Simple gun regularisation for 3d molecular property prediction and beyond. In International conference on learning representations, 2021.
Joe G Greener and David T Jones. Differentiable molecular simulation can learn all the parameters in a coarse sarsse-rriied force field for proteins. PloS one, 16()))e0256990, 2021.
Thomas A Halgren. Merck molecular force field. basis, form, scope, parameterization, and performance of mmff94. Journal of computational chemistry, 17(5-6):490 519, 1996.
Teresa Head-Gordon Martin Head-Gordon, Michael J Frisch, Charles Brooks 111, and John Pople. A theoretical study of alanine dipeptide and analogs. International Journal of Quantum Chemistry, 36(S16):
311-322, 1989.
RL Henderson. A uniqueness theorem for fluid pair correlation functions. Physics Letters A, 49(3):197-115) 198, 1974.
Johannes Hoja, Leonardo Medrano Sandonas, Brian Ernst, Alvaro Vazquez-Mayagoitta, Robert A DiSta- sio Jr, and Alexandre Tkatchenko. Qm7-x, a comprehensive dataset of maanntum-mecanncal properties spanning the chemical space of small organic molecules. Scientific dala, 8(1):1-11, 2021.
15

---

Published in Transactions on Machine Learning Research (04/2023)
Lars Holdijk, Yuanqi Du, Ferry Hooft, Priyank Jaini, Bernd Ensing, and Max Welling. Path integral stochastic optimal control for sampling transition paths. arXiv preprint or arx::22077.02440 2022.
Weihua Hu, Muhammed Shuaibi, Abhishek Das, Siddharth Goyal, Anuroop Sriram, Jure Leskovec, Devi Parikh, and C Lawrence Zitnick. Forcenet: A graph neural network for large-scale quantum calculations. arXiv preprint arXix:1000.011 0166,, 2021.
John Ingraham, Adam Riesselman, Chris Sander, and Debora Marks. Learning protein structure with differentiable simulator. In International Conference on Learning Representations, 2018.
William Jorgensen and Julian Tirado-Rives. Potential energy functions for atomic-level simulations of water and organic and biomolecular systems. Procecdings of the National Academy of Sciences, 102(19): 6663-6670, 2005.
Adil Kabylda, Valentin
Stefan Chmicla, Igor Poltavsky, and Alexandre Tkatchenko. Towards linearly scaling and chemically accurale global machine learning force fields for large molecules. arXiv preprint Kkkr::2009009988 2022.
George A Kaminski, Richard A Friesner, Julian Tirado- Rives, and William L Jorgensen. Evaluation and reparametrization of the opls-aa force field for proteins via comparison with accurate quantum chemical calculations on peptides. The Journal of Physical Chemiisry B, 105(28):6474 6487, 2001.
Alireza Khorshidi and Andrew A Peterson. Amp: A modular approach to machine learning in atomistic simulations. Computer Physics Communications, 207:310 324, 2016.
Adeesh Kolluru, Muhammed Shuaibi, Aini Palizhati, Nima Shoghi, Abhishek Das, Brandon Wood, C Lawrence Zitnick, John R Kitchin, and Zachary W Ulissi. Open challenges in developing generalizable large scale machine learning models for catalyst discovery. arXiv preprint Xrxi:22006. 02005, 2022.
Dávid Péter Kovács, Cas van der Oord, Jiri Kucera, Alice EA Allen, Daniel J Cole, Christoph Ortner, and Gábor Csányi. Linear atomic cluster expansion force fields for organic molecules: beyond rmse. Journal of chemical theory and computation, 17(12):7696-7711, 2021.
Alessandro Laio and Michele Parrinello. Escaping free-energy minima. Proceedings of the National Academy of Sciences, (92)))122 12566, 2002.
Thomas J Lane, Gregory R Bowman, Kyle Beauchamp, Vincent A Voelz, and Vijay Pande. Markov state model reveals folding and functional dynamics in ultra-long md trajectories. Journal of the Americon Chemical Society, 133(45):18413- 18419, 2011.
Ask Hjorth Larsen, Jens Jorgen Mortensen, Jakob Blomqvist, Ivano E Castelli, Rune Christensen, 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 Leonhard Kolsbjerg, Joseph Kubal, Kristen Kaasbjerg, Steen Lysgaard, Jón Borgmann Maronsson, Tristan Maxson, Thomas Olsen, Lars Pastewka, Andrew Peterson, Carsten Rostgaard, Jakob Schiøtz, Ole Schütt, Mikkel Strange, Kristian S Thygesen, Tejs Vegge, Lasse Vilhelmsen, Michael Walter, Zhenhua Zeng, and Karsten W Jacobsen. The atomic simulation environment python library for working with atoms. Journal of Physics: Condensed Matter, 27777),,,,00,, 2017. URL http:/ /ttt/ksssskksss/// wtw.ooo/o......oooo org/0953- (1553-9984/9/1-221 273300.
Andrew R Leach and Andrew R Leach. Molecular modelling: principles and applications. Pearson education, 2001.
Jonas Lederer, Michael Gastegger, Kristof T Schütt, Michael Kampffmeyer, Klaus-Robert Müller, and Oliver T Unke. Automatic identification of chemical moieties. arXiv preprinl arXin:2203. 0et02203..200.. 16205, 2022.
Yi-Lum Liao and Tess Smidt. Equiformer: Equivariant graph attention transformer for 3d atomistic graphs arXiv preprint 2xXt::2206.11900 2022.
Kresten Lindorff-Larsen Stefano Piana, Ron O Dror, and David E Shaw. How fast-folding proteins fold. Science, (3055))55) 520, 2011.
16

---

Published in Transactions on Machine Learning Research (04/2023)
Yi Liu, Limei Wang, Meng Liu, Yuchao Lin, Xuan Zhang, Bora Oztekin, and Shuiwang Ji. Spherical message passing for 3d molecular graphs. In International Conference on Learning Represenlations, 2021.
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. arXiv preprint ar orXi::2204.05249, 2022.
Cheol Woo Park, Mordechai Kornbluth, Jonathan Vandermause, Chris Wolverton, Boris Kozinsky, and Jonathan P Mailoa. Accurate and scalable graph neural network force field and molccular dynamics with direct force architecture. ugj Compulational Malerials, 7(1):1-9, 9, 2021.
Jay W Ponder and David A Case. Force fields for protein simulations. Advances in protein chemistry, 66:
27 85, 2003a.
Jay W. Ponder and David A. Case. Force Fields for Protein Simulations. In Advances in Protein Chemistry, volume 66, pp. 27-85. Elsevier, 2003b. ISBN 978-0-12-0342866-2 doi:
URL
elaevier wlt.rvewr.ccmeeeeeeeeeeee///lles
Zhuoran Qiao, Anders S Christensen, Matthew Welborn, Frederick R Manby, Anima Anandkumar, and Thomas F Miller III. Unite: Unitary n-body tensor equivariant network with applications to quantum chemistry. arXie preprint arX 14655, 2021.
Aneesur Rahman. Correlations the motion of atoms in liquid argon. Physical review, 136(2A):A405, 1964.
Raghumathan Ramakrishnan, Pavlo O Dral, Matthias Rupp, and O Anatole Von Lilienfeld. Quantum chemistry structures and properties of 134 kilo molecules. Scientific data, 1(1):1-7, 2014.
David Rosenberger, Justin S Smith, and Angel E Garcia. Modeling of peptides with classical and novel machine learning force fields: A comparison. The Journal of Physical Chemistry B, 125( 3612, 2021.
Alvaro Sanchez-Gonzalee Jonathan Godwin, Tobias Pfaff, Rex Ying, Jure Leskovec, and Peter Battaglia. Learning to simulate complex physics with graph networks. In International Conference Machine Learning, pp. 8459-8468. PMLR, 2020.
Victor Garcia Satorras, Emiel Hoogeboom, and Max Welling. E (n) equivariant graph neural networks. In International conference on machine learning, pp. 9323 9332. PMLR, 2021.
Tamar Schlick. Molecular modeling and simulation: an interdisciplinary guide, volume 2. Springer, 2010.
Elia Schneider, Luke Dai, Robert Q Topper, Christof Drechsel-Grau, and Mark E Tuckerman. Stochastic neural network approach for learning high-dimensional free energy surfaces. Physical review letters, 119 (15):150601, 2017.
Samuel Schocnholz and Ekin Dogus Cubuk. Jax md: framework for diffcrentiable physics. Advances in Neural Information Processing Systems, 33:11428 11441, 2020.
Kristof Schütt, Picter-Jan Kindermans, Huziel Enoc Saucoda Felix, Stefan Chmicla, Alexandre Tkatchenko, and Mlaus-Robert Müller. Schnet: A continuous filter convolutional neural network for modeling quantum interactions. Aduances neural information processing systems, 30, 2017.
Kristof Schütt, Oliver Unke, and Michael Gastegger. Equivariant message passing for the prediction of tensorial properties and molecular spectra. In International Conference on Machine Learning, pp. 9377-9388. PMLR, 2021.
Daniel Schwalbe-Koda, Aik Rui Tan, and Rafael Cómer-Bombarelli. Differentiable sampling of molecular geometries with uncertainty-baset adversarial attacks. Nature communications, 12(1):1-12, 2021.
Hythem Sidky, Wei Chen, and Andrew Ferguson. Machine learning for collective variable discovery and enhanced sampling in biomolecular simulation. Molecular Physics, 17777))))737742, 2020.
17

---

Published in Transactions on Machine Learning Research (04/2023)
Justin S Smith, Olexandr Isayev, and Adrian E Roitberg. Ani-1: extensible neural network potential with dft accuracy at force field computational cost. Chemical science, 8(4):3192-322) 222-220,,-- 2017.
Sina Stocker, Johannes Gasteiger, Florian Becker, Stephan Günnemann, and Johannes T Margraf. How robust are modern graph neural network potentials in long and hot molecular dynamics simulations? Machine Learning: Science and Technology, 3(4):045010, 2022.
Mohammad M Sultan, Hannah K Wayment-Steele, and Vijay S Pande. Transferable neural networks for enhanced sampling of protein dynamics. Journal of chemical theory and computation, 14(4):1887 1894, 2018.
So Takamoto, Chikashi Shinagawa, Daisuke Motoki, Kosuke Nakago, Wenwen Li, Iori Kurata, Taku Watanabe, Yoshihiro Yayama, Hiroki Iriguchi, Yusuke Asano, et al. Towards universal neural network potential for material discovery applicable to arbitrary combination of 45 elements. Nature Communications, 13(1): 1 11, 2022.
Stephan Thaler and Julija Zavadlav. Learning neural network potentials from experimental data via differentiable trajectory reweighting. Nature Communications, 12(1):1-10, 2021.
Philipp Thölke and Gianni De Fabritiis. Equivariant transformers for neural network based molecular potentials. In International Conference Learning Representations, 2021.
Nathaniel Thomas, Tess Smidt, Steven Kearnes, Lusann Yang, Li Li, Kai Kohlhoff, and Patrick Riley. Tensor field networks: Rotation-and translation-eqquivarianl neural networks for 3d point clouds. arXiv preprinl arX2v:1802.088 2018.
Richard Tran, Janice Lan, Muhammed Shuaibi, Siddharth Goyal, Brandon M Wood, Abhishek Das, Javier Heras-Domingo, Adeesh Kolluru, Ammar Rizvi, Nima Shoghi, et al. The open catalyst 2022 (oc22) dataset. and challenges for oxide electrocatalysis. arXiv preprint arXiv:2206. 02906.0.8177 7, 2022.
Gareth A. Tribello, Massimiliano Bonomi, Davide Branduardi, Carlo Camilloni, and Giovanni Bussi. PLUMED 2: New feathers for old bird. Computer Physics Communications, 185(2):604-613, February 2014. ISSN 00104655. doi: 10.1016/j.cpc.2013.. pis/80010466513031966.
URL https://1inkingghh elsevier com/retrieve/ weweeee......eeeerrrrreeeeeee
Mark Tuckerman. Statistical mechanics: theory and molecular simulation. Oxford university press, 2010.
MBBJM Tuckerman, Bruce Berne, and Glenn Martyna. Reversible multiple time scale molecular dynamics. The Journal of chemical physics, 97(3):1990-2000 1992.
Oliver Unke and Markus Meuwly. reactive, scalable, and transferable model for molecular energics from a neural network approach based on local information. The Journal of chemical physics, 148(24)2241708, 2018.
Oliver T Unke, Stefan Chmiela, Michael Gastegger, Kristof T Schütt, Huziel E Sauceda, and Klaus-Robert Müller. Spookynet: Learning force fields with electronic degroes freedom and nonlocal effects. Nature communications, 12(1):1-14, 2021a.
Oliver T Unke, Stefan Chmicla, Huziel E Sauceda, Michael Gastegger, Igor Poltavsky, Kristof T Sch utt, Alexandre Tkatchenko, and Klaus-Robert Müller. Machine learning force fields. Chemical Reviews, 121 (16):10142-10186, 2021b.
Jonathan Vandermause, Steven B Torrisi, Simon Batzner, Yu Xie, Lixin Sun, Alexie M Kolpak, and Boris Kozinsky. On-the-fly active learning of interpretable bayesian force fields for atomistic rare events. nyj Computational Materials, 6(1):1 11, 2020.
Ercheng Wang, Huiyong Sun, Junmei Wang, Zhe Wang, Hui Liu, John Z. H. Zhang, and Tingiun Hou. End-Point Binding Free Energy Calculation with MM/PBSA and MM/GBSA: Stratogies and Applications in Drug Design. Chem. Rev., 119(16):9478 9508, August 2019a. ISSN 0009-2665, 1520-6890. doi: 10.1021 /0.1021/acs.chemmres www.scchemee..99000000 URL htttps://pule acs www.ooog/oo//o..oooo 1021//co0 cwemrev. 9b00055.
18

---

Published in Transactions on Machine Learning Research (04/2023)
Jiang Wang, Simon Olsson, Christoph Wehmeyer, Adrià Pérez, Nicholas E Charron, Gianni De Fabritiis, Frank Noé, and Cecilia Clementi. Machine learning of coarse-grained molecular dynamics force fields. ACS central science, 5(5):755 767, 2019b.
Wujie Wang and Rafael Kömez-Bombarele Coarse-graining auto-encoders for molecular dynamics. npj Computational Materials, 5(1):1-9, 2019.
Wujie Wang, Simon Axelrod, and Rafael KomezzBoomaarlll Differentiable molecular simulations for control and learning. arXiv preprint Xiv:20003 00868, 2020a.
Wujie Wang, Tzuhsiung Yang, William Н Harris, and Rafael Gómezz-ombarelli. Active learning and neural network potentials accelerate molecular screening of ether-based solvate ionic liquids. Chemical Communications, 56(63):8920 8923, 2020b.
Wujie Wang, Zhenghao Wu, and Rafael Maanex--onbarllli. Learning pair potentials using differentiable simulations. arXiv preprint arXiv:2209.0700 2022.
Michael A Webb, Yukyung Jung, Danielle M Pesko, Brett M Savoie, Umi Yamamoto, Geoffrey W Coates, Nitash P Balsara, Zhen-Gang Wang, and Thomas F Miller III. Systematic computational and experimental investigation of lichium-ion transport mechanisms in polyester-based polymer electrolytes. ACS central science, (4):198-2005 2015.
Yujie Wu, Harald L. Tepper, and Gregory A. Voth. Flexible simple point-charge water model with improved liquid-state properties. The Journal of Chemical Physics, 124(2):024503, January 2006. ISSN 0021-9606, 1089-7690. doi: 10. 1..063/..2136777 URL http://aip.sciieee http://iitaaiiiooorooo///o////// 1167777777777
Tian Xie, Arthur France-Lanord, Yanming Wang, Jeffrey Lopez, Michael A Stolberg, Megan Hill, Gra- ham Michael Leverick, Rafael Gomez- Bombbarlla,rlii Jeremiah Johnson, Yang Shao-Horn, et al. Accelerating amorphous polymer electrolyte screening by learning to reduce errors in molecular dynamics simulated properties. Nature communications, 13(1):1-10 10, 2022.
Li Yarnell, MJ Katz, Ro Go Wenzel, and SH Koenig. Structure factor and radial distribution function for liquid argon k. Physical Review A, ((6):2130, 1973.
Shuwen Yue, Maria Carolina Muniz, Marcos F. Calegari Andrade, Linfeng Zhang, Roberto Car, and Athanassios Z. Panagiotopoulos. When do short-range atomistic machino-learning models fall short? J. doi:
Chem. Phys., 154(3):034111, Jamuary 2021. ISSN 0021-9606, 1089-7690.
URL
http:://aap..citaatt.. http://ataaaaaooooor//o///......
org//oo//..oooooooo 1063/5.6..0. 0031215.
Yaoguang Zhai, Alessandro Caruso, Sigbjörn L Bore, and Francesco Paesani. A "short blanket" dilemma for a state-of-thc-art ncural network potential for water: Reproducing propertics or learning the underlying physics? 2022.
Linfeng Zhang, Jiequn Han, Han Wang, Roberto Car, and EJPRL Weinan. Deep potential molecular dynamics: a scalable model with the accuracy of quantum mechanics. Physical review letters, 120(14): 143001, 2018a.
Linfeng Zhang, Jiequn Han, Han Wang, Wissam Saidi, Roberto Car, et al. End-to-end symmetry preserving inter-atomic potential energy model for finite and extended systems. Advances in Neural Information Processing Systems, 31, 2018b.
19

---

Published in Transactions on Machine Learning Research (04/2023)
Aspirin, Classical force field
Ethanel, Classical force field
Naphthai one, Clsssical force field
Salicyic Acid, Classical force fiold
1.2
Classical Reference
Cassica Reference
2.5
Classical Reference
Classical Reference
0.8
2,0
0.8
5. 1.5
1.0
0.4
h/0
1000
0.5
0.5
M
0.0
0.0
0.0
0
100
200
300
400
600
0
100
200
300
400
600
0
100
200
300
400
500
0
100
200
300
400
500
[[4]
r(A]
[[A]
[(A]
Figure 7: Comparison between a classical force field and the reference simulation on the h(r) of MD17 molecules.
A
Dataset details
The MD17 ' (Chmiela et al., 2017) and the LiPS dataset? (Batzner al., 2022) are adapted from previous works and are publicly available. The MD17 datasct is generated from path-integral molecular dynamics simulations that incorporate quantum mechanics into the classic molecular dynamics simulations using Feynman path integrals. The LiPS datasets are generated by ab-initio molecular dynamics simulations with a generalized gradient PBE functional and projector augmented wave pseudopotentials. We refer interested readers to the respective papers for more details on the data generation process. The water dataset and alanine dipeptide dataset are generated by ourselves, and the generation process explained in the next. paragraphs.
To demonstrate how ML force fields can improve simulation accuracy over classical force fields, also simulated the four MD17 molecules with popular classical force field: the optimized potentials for liquid simulations (OPLS, Jorgensen & Tirado-Rives 2005) with OpenMM (Eastman et al., 2017). In OPLS, interatomic interactions are described by a combination of simple functional forms, such as harmonic and Lennard-Jones potentials, for bond stretch, bond angle, torsional angle, and pair-wise non-bonded interactions. We sample the molecule conformations in vacuum under an NVT ensemble with a timestep of femtosccond and lemperature of 500 K, controlled with a Langevin thermostat. The resulting h(r) from the classical force field is compared to the reference in Figure 7. The errors it attains are reported in Table 3, which are much higher than most ML force fields, whose h(r) curves are shown in Figure 14.
Water. Our water dataset is generated from molecular dynamics simulations of a simple classical water model, namely, the flexible version of the Extended Simple Point Charge water model (SPC/E-fw) (Wu et al., 2006) at temperature T 300 K and pressure P = atm. For this model, the interaction parameters (e.g., O-H bond stretch and H-O-H bond angles), are parameterized to match extensive experimental properties such the self-diffusion and dielectric constants at bulk phase. This classical model has been well-studied in previous work (Wu et al., 2006; Yue et al., 2021) and has shown reasonable predictions of the physical properties of liquid water. It provides a computationally inexpensive way to generate a large amount of training data. The experience and knowledge gained from the benchmark based on the simple model can be readily extended to systems with higher accuracy, such as the ab-initio models.
Alanine dipeptide. Our datasel is generated from the MD simulation of an alanine dipeptide molecule solvated in explicit water (1164 water molecules) performed in GROMACS (Abraham al., 2015) using the AMBER-03 (Ponder & Case, 2003b) force-field. In the AMBER-03 force field, the potential energy parameters such as van der Waals and electrostatics with minor optimization on the bonded parameters to reproduce the experimental vibrational frequencies and structures (Cornell et al., 1996; Ponder & Case, 2003a). The NPT ensemble is applied in simulations, with hydrogen bond length constraints using LINcar Constraint Solver (LINCS) and a time step of fs. The temperature and pressure of the system are controlled at T = 300 K and P = 1 bar using a stochastic velocity rescaling thermostat with damping frequency f, = 0.1 ps and Parrinello-Rane Rahman barostat with coupling frequency 4p 2.0 ps, respectively. The Particle Mesh Ewald approach is used to compute long-range
mostly derived from quantum mechanical methods
http:::///www..nn...m 2https://archive. wtt.:e/ttrrraaeeeoooooeoeeeeeraa crt//recood///222..5
20

---

Published in Transactions on Machine Learning Research (04/2023)
(a)
(b)
40
40
10
30
30
F[KJimol]
20
10
[[J]moll]
20
20
10
10
10
MB
10 ms
0
-180
-90
160
-180
-90
180
-180
-90
90
180
-180
-90
0
90
180
[degrees]
4 [degrees]
[degrees]
[degrees]
Figure 8: F()) and F(D) have converged for the reference force field (a) and NequlP (b) at time 5 ns under Metadynamics.
electrostatics with periodic boundary conditions applied to the x, y, and z directions. The conformational modes can be characterized by six free energy local minima, which have been used in previous work (Lederer et al., 2022). We initialize six simulations for each model from each of the six free energy local minima.
Implicit solvation. The explicit solvent of 1164 water molecules not the subject of study but adds a significant computational burden. In this task, we attempt to learn an implicit solvent model (ISM) of the alanine dipeptide, in which the explicit solvent environment is incorporated in the learned FF. The ISM is commonly used in drug design (Wang et al., 2019a) because it can speed up the computation by dramatically decreasing the number of particles required for simulation. In general, the mean-field estimation in ISM ignores the effect of solvent, thermal fluctuations, and solvent friction (Feig, 2007). Thus, molccular kinctics is not directly comparable to the explicit solvation simulation. However, the equilibrium configurations can be explicitly compared, as conducted in Chen al. (2021).
Metadynamics simulation. Simulating energy barrier jump usually requires long sampling of the trajectory in MD simulations. The conformational change of alanine dipeptide in water involves such a process, making it difficult to extract the complete free energy surface, i.e., the Ramachandran plot, in normal MD. In order to examine the learned ML FFs within a reasonable time limit, metadynamics (Laio & Parrinello, 2002) is employed to explore the learned FES of the solvated alanine dipeptide. Metadynamics is a widely used technique in atomistic simulations to accelerate the sampling of rare events and estimate the FES of a certain set of degrees of freedom. It is based on iteratively "filling" the potential energy using sum of Gaussians deposited on set of suitable collective variables (CVs) along the trajectory. At evaluation time, perform metadynamics with dihedral angles é and vi CVss', starting from the configurations located at one of the six energy minimums in the free energy surface indicated in Figure 5 (c). The Gaussians with height h = 1.2 and sigma a = 0.35 are deposited every 1 ps centered on d and Φ. As shown in Figure 8, the estimated FES of both @ and eh do not significantly change after ns. In addition, the height of the bias gaussian potential smoothly converges to
the time limit of ns. Therefore, simulation time of 5 ns is sufficient for the convergence of the metadynamics. This metadynamics simulation of alanine dipeptide with AMBER force fields is carried out using GROMACS (Abraham et al., 2015) integrated with the PLUMED library (Tribello et al,, 2014; Colón-Ramos et al., 2019) of version 2.8.
B Experimental details
Selection of stability threshold. Our stability thresholds are chosen to be relaxed, so a simulation is only flagged as "unstable" when the system has already went into highly non-realistic configurations. In Figure 9 (a) and (c), we show the natural fluctuation existing in the reference data of RDF for water and bond lengths for Aspirin, demonstrating that realistic simulations would never exceed the stability thresholds. In our experiments, unstable simulations cannot recover from catastrophic failure become stable again. Example trajectories using SchNet on Aspirin and Gem-et on water are shown in Figure 9 (b) and (d).
In practice, the selection of suitable collective variables can be a case-hy-case challenge (Sidky al., 2020).
21

---

Published in Transactions on Machine Learning Research (04/2023)
(1)
Plactuation of Asp ein bond langhs
(b)
Aspitn, Schible
(c)
Plachatiion of RDF
(d)
Water 1k. GartNNtt
Band Langth Diviafoon
10
SchNet
Reference
RDF Doviafta
Referonce
GarNee-T
Danasyy
Threshold
10
011
022
0.9
0.4
0.6
20
40
60
80
100
0
0.5
1.0
1.5
2.0
2.5
3.0
60
100
150
200
250
Devistion
Simutation Time (pa)
Devitton
Samu inion Time [po]
Figure (a) The distribution bond length deviation for the Aspirin reference dataset. Black dashed line is our chosen stability threshold. (b) An example simulated trajectory of Aspirin using SchNet that becomes unstable. Y-axis is log-scaled. (c) The distribution of RDF deviation for the water reference dataset. Black dashed line is our chosen stability threshold. (dl An example simulated trajectory of Water- 1k using GemNet-T that becomes unstable.
Translation
Rotation
Reflection
Figure 10: 3D symmetries that ML force fields should respect include translations, rotations, and reflections.
Symmetry principles. the learned foree fields should respect the symmetry principles: energy is permu- tational invariant and E(3)-invariands and forces are permutational invariant and E(3)-equivariant. E(3) comprises translations, rotations, and reflections 3D (Figure 10). A function f is G-equivariant if foT = Lof for some operator I,, where T an operator for transformation in G. f is invariant f oT f. For example, the energy prediction of a E(3))invariann model will be invariant with represet to 3D rotation, translation, and reflection of the input structure. Our benchmark includes various models with different levels of symmetry principles (Table 2) and thus have different expressive power.
Experimental procedures. Bascline models are trained on the datascts described in Section 4 according to the experimental settings described in Appendix B. At evaluation time, we simulate MD trajectories using the learned models, with thermostats and simulation length described in Section 4. The simulated trajectories are recorded as time series of atom positions, along with other information including atom types, temperature, total energy, potential energy, and kinetic energy. All observables described in this section can be computed from recorded trajectories. We use the stability criterion described above to find the time step when the systems become "unstable", and only use the trajectory before that time step for the computation of obscrvables. Among the observables, the distribution of interatomic distances and RDF are computed for each frame and averaged over the entire trajectory. Diffusivity coefficients are computed by averaging over the diffusivity coefficient computed from all applicable time windows where the time window length is predefined to be long enough. For example, for a trajectory of T steps and a time window of steps, we average over ................. T K + 1, We use 100 ps as the
the diffusivity computed from the time windows: KJ, time window size for water and 35 ps the time window size for LiPS. We also remove the first. 5 ps of the simulated trajectories of LiPS for equilibriun. As do multiple simulations per model and dataset, we compute the metrics for each trajectory and report the mean and standard deviation. When reporting the efficiency of different models, all frames per second (FPS) metrics are measured with an NVIDIA Tesla V100-PCIe GPU. We present FPS as a reference for models' computational efficiency but also note that code speed can be affected by many factors and likely has room for improvement. Further details on observable computation can be found in our code submission: observable. ipynb. The Open Catalyst Project codebase
22

---

Published in Transactions on Machine Learning Research (04/2023)
Table 7: Default training-related hyperparameters for cach dataset. "We adopt the original batch size from respective papers when available for MD17. DeepPot-SE: 4; SchNet: 100; DimeNet: 32; GemNet-T/dT: 1;
NequlP: 5. We use a batch size of 8 for ForceNet.
Dataset
Training dataset size
Batch size
Max epoch
LR patience
Longest simulation time
MD17
9,500
1-100"
2,000
5 epochs
20 hours
Water-1k -kk
950
10,000
50 epochs
28 hours
Water-10k
9,500
1
2,000
5 epochs
28 hours
Water-90k
85,500
1
400
3 epochs
28 hours
Alanine dipeptide LiPS
38,000
5
2,000
5 epochs
75 hours
19,000
1
2,000
5 epochs
7 hours
Table 8: Results on Water-1k. Force MAE is reported in the unit of IneV/Aj, Stability is reported in the unit of [ps]; Diffusivity MAE is reported in the unit of [10-9 m² /s]; RDF MAE and FPS are unitless.
DeepPot-SE
SchNet
DimeNet
PaiNN
SphereNet
ForceNet
GemNet-dT T
NequIP
Force (4)
6.7
13.1
27.5
13.6
29.2
Stability (f)
108(117)
(18))
21
111:63)
500gg)
[19
(91
388
500go)
(o. (L)
0.17
(.52.0.03)
0.18
(.08
(.67
0.95
0.25 22)
0.50
0.07
(D.( (L)
0. 0.13
(.23
0.1 0.14
(.07
0.29
(.60
0. 1.19
0.36
(.07
(n. (C)
(.28
(..7
0.24
(.19(o.06)
(.87
1.53
0. (.42(0.20)
1.01:o.10)
0.26
Diffusivity (L)
0.24
1.70
0.50
2.09
2.29
0.06
1.11
0.37
FPS (t) 61.8
09.2
16.4
54.5
3.0
68.1 1
15.4
34.5
3.9
and the official codebases of DeepPot-SE", SphereNet and NequIP" are all publicly available. We build our MD simulation framework based on the Atomic Simulation Environment (ASE) library (Larsen et al., 2017).
Hyperparameters. We adopt the original model hyperparameters in the respective papers and find they can produce good force prediction results that match the trend and numbers for MD17 reported in previous work. As we introduce new datasets, we set training hyperparameters such as the batch size and summarize them in Table 7. For water and LiPS, each structure already contains reasonable number of atoms and interactions. Following previous work, we use an initial learning rate of 0.001 for all experiments except for NequlP, which uses 0.005 as the initial learning rate in the original paper. For models that minimize a mixture of force loss and energy loss, we set the force loss coefficient ÄF to be 1000 and the energy loss coefficient Xø to be 1, if not specified in the original paper. higher force loss cocfficient common previous work (Zhang et al., 2018a; Batzner et al., 2022) as simulations do not directly rely on the energy.
batch size of like in previous work (Batzner et al., 2022)
Notably, NequlP' proposed several sets of hyperparameters for different datasets, including MD17, water+ice dataset from Zhang et al. 2018a, LiPS, etc. We follow the MD17 hyperparameters of NequlP for our MD17 and alanine dipeptide datasets; the water+iie hyperparameters of NequlP for our water dataset; and the LiPS hyperparameters of NequIP for the same LiPS dataset. For DeepPot-SE, we adopted hyperparameters introduced in Zhang et al. 2018a. The only architectural adjustment we made is because we observed training instability for ForceNet on water using the original hyperparameters. We resolve this issue by reducing the network width from 512 to 128 for ForceNet in our water experiments.
To facilitate benchmarking with a reasonable computational budget, stop the training of a model if either of the following conditions is met: (1) maximum training time of 7 days is reached on NVIDTA Tesla V100-PCle GPU; (2) a maximum number of epochs specified in Table 7 is reached; (3) The learning rate drops below 10-6 6 with a ReduceLROnPPateau scheduler with factor 0.8 and learning rate (LR) patience specified in Table 7. We also report the longest time for an ML model to finish our benchmark simulation in Table 7. All numbers are results of NequIP. The high computational cost for evaluating MD simulations has
'https://gitaee httg:hhhwb.ccm//eemmmda////////e httpp/////ggaubb..cm////llllaaa
hhttpe://gtthub..comiii-rr httap:///gitaab.com/astt
23

---

Published in Transactions on Machine Learning Research (04/2023)
Table 9: Results on Water-90k. Force MAE is reported in the unit of IncV/A) Stability is reported in the unit of [ps]; Diffusivity MAE reported in the unit of [10-9 m² /s]; RDF MAE and FPS are unitless.
DeepPot-SE
SchNet
DimeNet
PaiNN
SphereNet
ForceNet
G
GemNet-dT
NequIP
Force (4)
5.9
14.8
0.7
1.1
Stability (f)
(00
(00
1877
500
g00
20/2)
500(o)
500
500ga)
(o.0) (L)
0.07
0.64
0.08
0.
(.90)
(.77(0.04)
0.04,0.00)
0.04:o.01)
0.06go.01)
(m, (U)
0.05
0.28 0.28
0.06
0.08
0.41
0.30
0. 0.03g0.00)
0. 0.04ço.00)
0.04:o.01
RDF (L)
0.29(o.08)
0. 0.06)o.02)
(.
0.
1.26
1.53
0. 0.19/o.o2)
0. 0.19
0.25
Diffusivity (L)
0.35
1.89
1.06
2.
2.24
0.39
0.94
0.18
FPS (1)
62.1
103.0
16.3
71.9
3.1
43.8
15.3
32.7
3.0
Table 10:
Water- 1k results on NequlP with various model sizes and radius cutoffs.
Force
Stability
RDF (0.0)
RDF (H.H)
RDF (H,O)
Diffusivity
FPS
Width-64, r-4
3.5
(00
(.07()
(.05
(0.00)
0.38
8.2
Width=32, r=6
1.5
500(0)
0.06
(.05(..01)
0.26(0.06)
0.25
5.2
Width-64, r=5
1.6
(00
(..0))
(.01))
0.27
0.31
4.9
Width=64, r=6
1.4
(00
(.07()
0.07().02)
(.26
0.37
3.9
Width=128, r-
1.5
(00
0.07()
0.05(0.01)
0.29
0.37
2.5
been major consideration in designing our benchmark datasets and metrics. Training of DeepPot-SE. is efficient and we follow the training selup specified in Zhang et al. 2018a.
Complete water results. We present results water- 1k in Table 8 and results on water-90k in Table 9. Results on water-10k is presented in Table in the main text. All models generally achieve lower force error, higher stability, and higher simulation accuracy when trained with more data. In particular, DeepPot-SE shows clear improvement with more training data and becomes good as NequlP water-90k while being highly efficient. SchNet demonst strong stability, but the estimation of ensemble statistics does not improve with more data. This may be due to the limited accuracy of SchNet coming from the limited expressiveness of the invariant atomic representation.
Force MAE
35
Stability
300
Force MAE [meVIA]
30
250
Stability [PS]
25-
200
20
150
15-
100
10-
50
100
200
300
400
500
Training Epoch
Large water system of 512 molecules. To study model performance in generalizing to larger system and model scalability, we evaluate models trained on the water-10k dataset on a dataset with 512 water molecules simulated for ns, using the same reference force field. Given the high cost simulating a large system, we simulate trajectories of 150 ps long for each model. The results are shown in Table 11. We observe that all models suffor slightly higher force crrors compared the evaluation of the 64-molccule water system. In terms of stability, NequlP, SchNet, and SphereNet always remain stable for the entire 150 ps. However, SchNet and SphereNet do not produce correct ensemble properties. GemNet-T and GemNet-dT are not stable throughout the entire simulation but can produce decent ensemble properties. Noticoably, the stability of DeepPot-SE drops siguificantly. We hepothesize that the lack of message passing limits its capability to capture long-range interactions and thus limits the performance in a larger system.
Figure 11: The force error and stability of SchNet. for simulating Salicylic acid, as training progress.
Influence of model size. Table 10 shows an ablation study over the model size and radius cutoff of NequlP over water-1k: We observe that all models are highly stable and attain equally good performance in simulation-basel metrics. Although small radius cutoff of leads to worse performance in force prediction,
24

---

Published in Transactions on Machine Learning Research (04/2023)
Table 11: Results on the large water system with 512 molecules, with models trained on the water-10kt dataset (64-molecule water system). "SphereNet on Water-large requires more memory than Tesla V100 supports. We run its simulations on faster NVIDIA A100 cards so the FPS is not entirely comparable to other models.
DeepPot-SE
SchNet
DimeNet
PaiNN
SphereNet
ForceNet
GemNet-TT GemNet-dTT
Nequll
Force (4)
10.6
12.1
5.1
9.7
18.4
13.2
5.6
4.2
7.7
Stability (f)
19(2))
150(o)
19,24)
24(19)
(50
8to)
84,48)
124
150
oo. (1)
0 0.23g%.0)
(.62(0.01)
(.21
0.14
(.92
0.69
0. 0.12
0.10
0.10
am. (L)
0.24%%.00
0.28/0.00)
0.15
0.10go.o2
0.41(0.01)
0.48
0.08
0.07:o.02)
0.07:o.00)
RDF (H.(H. (-)
0.67
0.55
(.
0.
0.96
(.32
0. 0.14(0.03)
0. 0.13
0.12
Diffusivity (L)
2.63
2.99
0.81
0.60
0.89
FPS (1)
80.7
23.1
3.5
17.4
0.8
11.9
2.2
5.3
0.7
Table 12: Results on Water-10k Wateer-0k-time-split. Force MAE is reported in the unit of [meV Stability is reported in the unit of [ps]; Diffusivity MAE is reported in the unit of [10-9 9 m²/s] RDF MAE and FPS are unitless.
DeepPot-SE
SchNet
DimeNet
PaiNN
SphereNet
ForceNet
GemNet-T
GemNet-dTT
NequlP
Force (4)
5.8
10.1
1,2
6.6
18.4
12.3
0.7
1.3
1.7
Stability (f)
(152)
414
117
154:22)
(00
2g0;
(00
500gg)
500gg)
(o. (L)
0.08go.03)
(o73
Q.14
0.07
0.94
0.75
0.06go.01)
0.13
0.05
(R, (L)
0.06
(o.3))
(.
0.05
(o.1
0.55
0.04g0.00)
0.18
0.04:0.00)
(m.o (-)
0.18
0.76
(a.ss
0.24,0.05)
1.22(0.03)
1.28
0.24g0.0)
0. 0.62
0.21,0.00
Diffusivity (L)
0.15
1.97
0.60
0.57
2.25
0.34
0.77
0.25
FPS (f)
62.0
100.5
17.1
57.3
2.9
62.2
15.1
32.5
3.5
it is more computationally efficient and preserves the trajectory statistics. These results show that there exists a trade-off between accuracy and elficiency when choosing the hyperparameters of an ML. force field, and force error may not be the preferred criterion for model selection.
Stability's relation with dataset size. We extract the force and stability results for each model from Figure to make Figure 12 to better illustrate the relation between stability and dataset size for each model. We observe that more data almost always reduce force error and improves stability. In particular, NequIP is highly stable across all dataset sizes. DeepPot-SE, SchNet, DimeNet, PaiNN, GemNet, and GemNet-T have significant improvements in stability with more data. Section 7 contains more detailed discussions on the causes of instability and potential solutions to improve stability.
Stability's relation with training epochs. We study the evolution of simulation stability in the training process of an MT. force field. We take the SchNet model on the MD17 molecule salicylic acid and save the checkpoint at 100, 200, 300, 400, and 500 epochs. We conduct simulations of 300 ps using each checkpoint. Figure 11 shows the force error and stability of the model at different stages of training. We observe that the force error decreases as training progresses, and the stability improves to be stable across the entire 300 simulation and training epoch 300. This result reveals that thorough training is important to both the accuracy and stability of ML force fields.
Water-10k with a time split. We investigate a time split of the water dataset by using the first 10,000 structures for training and the last 10,000 structures for testing. The results are reported Table 12. We find most models perform slightly worse compared to the random split results in Table 4, while the trend of performance ranking stays the same.
Distribution of intoratomic distances for MD17. Figure 14 shows the h(r) curves for all models and molecules benchmarked. We randomly selected one simulation out of the five simulations we conducted for each model and molecule. We observe that due to lack of stability, DeepPot-SE produces noisy h(r) on Aspirin. ForceNet does not manage to learn the correct interatomic interactions and produces incorrect h(r) curves. Most models are able produce h(r) that match well with the reference, with SchNet being less accurate Aspirin and Ethanol.
25

---

Published in Transactions on Machine Learning Research (04/2023)
Force MAE (meV/A)
Water, DeepPot-SE
Water, SchNet
Water, DimeNet
Foice
Sibllty
10
Sualiyy
Sibilyy
Stability (ps)
4-
200
200
200
100
Water, PaiNN
Water, SphereNet
Water, ForcaNet
Force MAE [melWAA
20
Ferce
Sdally
Etbiity
200
Stability (ps)
200
200
100
Force MAE [meW[A]
Water, GamNet-T
30
Water, GemNet-dT
1.5
Water, NequiP
Fasze
Btal
400
Stality
Builky
Stability [pe]
800
10
0.5
202
Training Data Size
Training Data Size
Training Data Size
Figure 12: The force error and stability of all models on the water datasct, with varying training datasct size: 1k, 10k, and 90k.
NequiP
ComNeet
Reference
NequiP
Fullmo
GomNetT
Fummo]
[00grod)
10
-00
[eggees
190
6R
100
(degree)
100
NesulP
Duammol)
GemNus-
Brogree)
20
20
10
400
20
130
22
-30
13
(02
-0
30
182
.103
40
30
130
[degres]
(degree)
[degrea]
(degree)
Figure 13: Alanine dipeptide FES results from different initialization.
RDFss for water. Selected RDF curves for water- fater-1k/10k, 10k//0k are in Figure 15, Figure 16, and Figure 17. Most noisy curves are due to insufficient sampling time, which results in small number of frames to be averaged in obtaining the RDF curves. We observe that SchNet, SphereNet, and ForceNet produce inaccurate curves that are not very noisy, showing that their failure not entirely due to a lack of stability but because of inaccurate modeling of interactions caused by limited expressiveness and lower sample efficiency. Further, we that the reference curves have zero values below a certain threshold, any pair of atoms cannot get too close to each other. However, GemNet-T water-10k exhibits abnormally high values for very small distances, indicating the simulations have gone into nonphysical configurations and collapsed.
RDFs for LiPS. As shown in Figure 18, DeepPot-SE does not manage to stay stable on LiPS. ForceNet learns inaccurate interactions and produces inaccurate RDFs. All other models can produce highly accurate RDF and can reproduce Li-ion diffusivity relatively accurately, as demonstrated in Table
26

---

Published in Transactions on Machine Learning Research (04/2023)
DeepPot-SE, Aspirin
DeepPot-SE, Ethanel
DeepPot-SE, Naphthalene
DeepPot-SE. Salicylle Adi
Refurnne Predittion
1.5
Refarnne Prediction
0.8
Refarnnce Pnadittion
Refrence
h(r) [4
0.4
0.5
0.2
0.2
0.0
0.0
0.0
0.0
10
10
16
SchNet, Aspirin
SchNet, Ethanol
SchNet, Naphthalene
SchNet, Salicylic Acid
Raterence
Relerence Predision
0.6
Reteranc Predidion
0.8
Balerrcc
0.6
Precediio
hig [4-
0.4
0.2
0.5
0.2
0.2
DimeNet, Aspirin
DimeNet, Ethanol
DimeNst. Naphthalone
DimeNet, Salicylic Acid
0.8
Reforence Prodidtion
Reference Prodittion
Referunoe Prodidioon
Reforanc Prodiction
1.5
0.6
0.4
0.4
0.
0.2
0.2
0.8
0.0
0.0
0.0
10
PaINN, Aspirin
PaiNN, Ethanol
PaiNN, Naphthalene
PaiNN, Salicylic Acid
0.8
0.8
1.5
0.0
0.6
1N0 (A
6-
0.4
0.4
0.5
0.6
0.0
10
S Netere Aspirin
SphereNet, Ethanol
SphereNet Naphthalene
SphereNet, Sallcylic Acid
0.8
Haference Prediction
Reference
Haferance Prediction
Halorarce
Prediction
0.6
Pradidion
0.6
0.4
C.4
0.5
0.2
0.0
ForceNet, Aspirin
ForceNel, Ethanol
ForcaNat, Naphthalone
ForceNet, Salicylic Acid
1.5
Refernnce
Reforance
Refernnce
1.00
Refurenc
h(r) [A -1
0.75
1.0
0.00
10
10
GemNeetTT Aspirin
GemNet-T, Ethanol
GemNet-T. Naphthalene
GemNettT, Salicylic Acid
Rofornce
Rofurnoe
Rofurunce
0.8
Rofurnce
0.4
Prodition
Prodiction
0.6
Pnoditioo
on
1.0
044
0.4
02
0.2
0.2
GemNet-dT, Aspirin
GomNettdT, Ethanol
GemNet-dT, Naphthalane
GomNet-dT, Salicylic Acid
Relerenc Predicioo
1.5
Relerenc Prediciion
0.8
Ralennca Precidion
Relemece
0.6
Predidion
Ney (4
1.0
0.4
02
0.5
0.2
0.2
0.0
0.0
0.0
0.0
10
NequiP., Aspirin
NequiP. Ethanol
NequiP, Naphthalene
NequiP. Salicylic Acid
0.8
Referance Prodiction
Reforance Prodictio
Referance Prodition
Reforunce
0.6
Procicioo
0.4
0.8
0.0
10
T [A]
(A]
[A]
[A]
Figure 14: h(()) of md17.
27

---

Published in Transactions on Machine Learning Research (04/2023)
DeepPot-SE, 0-0
DeepPot-SE, H-H
DeepPot-SE H-O
Reference Prediction
Reference Prediction
Reference Prediction
SchNet, O-O
SchNet, H-H
SchNet, H-O
Reference Prediction
Reference Prediction
Reference Prediction
DimeNet. O-O
DimeNet, H-H
DimeNet. H-O
Reference Prediction
Reference Prediction
Reference Prediction
RDDD
PaiNN, O-O
PaiNN, H-H
PaiNN, H-O
Reference Prediction
Reference Prediction
Reference Prediction
M 00.004
SphereNet, O-O
SphereNet, H-H
SphereNet, H-O
[ MDgH
Reference Prediction
Reference Prediction
Reference Prediction
ForceNet,
ForcaNet. H-H
ForceNet,
Tr A 1elaaa
Reference
Reference Prediction
Reference Prediction
Prediction
GemNet-T, 0-O
GemNet-T,
GemNet-T, H-D
Reference
Reference Prediction
Reference Prediction
Prediction
GemNet-dT. 0-0
GemNeetdT,, H-H
GemNet-dT. H-O
Reference Prediction
Reference Prediction
Reference Prediction
NequiP, O-O
NequiP, H-H
NequiP, H-O
Reference Prediction
Reference Prediction
Reference Prediction
A
RDF
[A]
Figure 15: RDFs of Water-1k. 1k. GemNet-dT does not remain stable for more than 1 ps and therefore not feasible for RDF computation.
28

---

Published in Transactions on Machine Learning Research (04/2023)
DespPot-SE, 0-0
DeepPot-SE, H-H
DeepPot-SS SE, H-O
Reference Prediction
Reference Prediction
43
Reference Prediction
RDFI0 (8
SchNet. O-O
SchNet, H-H
SchNet, H-O
Moriy (A -1
Reference Prediction
Reference
Reference Prediction
Prediction
DimeNat,
DimeNet,
DimeNet, H-O
Reference Prediction
RDF (
Prediction
Pred ction
PaiNN, O-O
PaiNN, H-H
PaiNN, H-O
Reference Prediction
Reference Prediction
Reference Prediction
SpheroNet, O-O
SphereNot, H-H
SpheroNet, H-O
Reference Prediction
Reference Prediction
Reference Predi
ForceNat, 0-O
ForceNat,
ForcaNet, H-O
-1
Reference Prediction
Reference Prediction
Reference Prediction
14
RDF
20
GemNat-T,
GemNet-T,
GerNeet-,, H-O
Reference Prediction
Reference Prediction
Reference ction
GemNet-dT, 0-O
GemNet-dT, H-H
GemNet-dT, H-O
b
Reference
Reference
Reference Prediction
RDF A4
Prediction
Prediction
NequiP 0-O
NequIP, H-H
NequiP. H-O
RDR() [A-3) []
Reference
Reference Prediction
Reference Predictionn
Prediction
[A]
[A]
[A]
Figure 16: RDFs of Water-10k.
29

---

Published in Transactions on Machine Learning Research (04/2023)
DespPot-SE, 0-0
DeepPot-SE, H-H
DeepPot-Se SE, H-O
Reference Prediction
Reference
Reference
PDFW0 (4
Prediction
Prediction
SchNet. O-O
SchNet, H-H
SchNet, H-O
Mora (A -1
Reference Prediction
Reference
Reference Prediction
Prediction
20
DimeNat,
DimeNet,
DimeNet, H-O
Reference Prediction
RDF (
Prediction
Pred ction
20
PaiNN, O-O
PaiNN, H-H
PaiNN, H-O
Reference Prediction
Reference
Reference
PoRRe 14
Prediction
Prediction
20
SpheraNet, 0-0
SphereNot, H-H
SpheroNet, H-O
Reference Prediction
Reference Prediction
Reference Predi
Porvy 44
ForceNat, 0-0
ForceNat, H-H
ForcaNet, H-O
-1
Reference Prediction
Reference Prediction
Reference Prediction
IDFMM 14
GemNat-T,
GemNet-T,
GerNeet-,, H-O
Reference Prediction
Reference Prediction
Reference ctionctionn
10r1 14
20
GemNet-dT, 0-O
GemNet-dT, H-H
GemNet-dT, H-O
Reference
Reference Prediction
Reference
RDF A4
Prediction
Pred ction
20
NequIP. 0-O
NequIP, H-H
NequiP. H-O
RDR() [A-3) --
Reference
Reference Prediction
Reference Predicttonn
Prediction
[A]
[A]
[A]
Figure 17: RDFs of Water-90k.
30

---

Published in Transactions on Machine Learning Research (04/2023)
DeepPot-SE
SchNet
DimeNet
2.5
Refarence
Reference
2.5
Rafarence
Prediction
Prediction
Prediction
2.0
2.0
2.0
RDF( 14
1.5
1.0
1.0
1.0
0.5
0.0
0.0
PaiNN
SphereNet
ForceNet
2.5
Rafarence Prediction
Reference
Refaronce
RDF(1) [A-1]
2.0
2.0
Prediction
Prediction
1.5
1.5
1.0
0.5
0.5
2
2
GemNet-T T
GemNet-dT
NequlP
2.5
Reference Prediction
2.5
Reference
RDF(1) [A- -A-1
2.0
Prediction
Prediction
1.5
1.5
1.5
1.0
1.0
0.6
0.5
0.0
0.0
0.0
r [A]
r [A]
[A]
Figure 18: RDFs of LiPS.
31


Reviews:
[{"review_id": "EJ2cGY68Mk", "reviewer": "Reviewer_21Df", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "One major strenght of this paper is that it highlights a critical issue in the field of MD simulations, which is that even small errors in the force field parameters can lead to instabilities in certain integration schemes.  These instabilities can cause the simulation to become invalid, rendeirng all subsequent results meaningless. \n\nHowever, the paper's definition of \"stability\" in the appendix may be confusing to non-experts, leading them to believe that this problem is not a significant concern.  A more detailed explanation of the underlying casues of instabilities in symplectic integrators, or a connection to sampling frequency theorems, or a stronger connection to the early literature of molecular dynamics simulation that addressed such questions could have helped to a degree in strengthening the paper's main argument.  Additionally, a more rigorous demonstration that once an \"instability\" occurs, the MD simulation will basically collapse, would have reduced any possible concerns about the importance of these observations.  The current criteria for assessing stability may be misleading to non-experts: the RDF criterion appears to be related to the kinetics of sampling, while the bond-breaking criterion, while sufficient for detecting undesirable bond-breaking behavior, would also flag a valid tautomerization reaction as unstable (e.g. the hopping of the proton in aspirin).  The discussion of the accuracy of the RDF very much obscures the main objective of the paper, which is the characterization of the ability to perform MD simulations with these models, which is of fundamental importance and potentially orthogonal to the perceived accuracy.  It is not very clear if any remaining discrepancy between force errors and RDF errors are not simply due to the lack of the ability to perform stable simulations.\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "\nThis paper examines several machine learning force fields, with an emphasis on modern graph-neural network force fields that are commonly used in the materials science community. Specifically, the paper evaluates the stability and accuracy of radial distribution functions coming from molecular dynamics (MD) simulations that use these force fields in a number of different benchmark systems.  The stability of MD simulations is a critical evaluation criterion for accurate simulations over extended timeframes.  With the exception of studies in an isolated system, previous comprehensive benchmarks have only evaluated the mean force errors of such force fields. To address this gap in previous studies, the authors have developed a comprehensive benchmark that includes multiple models from the Open Catalyst benchmark, as well as Nequip and demonstrate that the force error metric alone is insufficient for predicting the ability to perform stable long-running MD simulations with such machine-learning force fields.  \n", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Perhaps the authors could clarify from their data if all of the unstable events lead to a complete breakdown of the subsequent simulations, or if some of these instabilities they've observe might be self-correcting fluctuations (which I very much doubt).\n\nAs a side note I'd point out that in addition to the self-correcting force fields, people have already used variants of the Behler Parrinello, and ANI to run long MD simulations of much larger systems, so arguably solutions to all these problems already exist in earlier MD-simulation focused force fields.  I don't expect the authors to add any model under the planet, however, it may help if the authors could sharpen the definition of their criteria for selecting the starting materials for their benchmarks (even if these criteria currently only involve practical considerations, such as the uniformity of the opencatalyst benchmark models.)\n", "broader_impact_concerns": "No ethical implications.  \n", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "TFiMz8rUoR", "reviewer": "Reviewer_vsLs", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n\nI appreciate how thorough the authors are in comparison across all methods. As a researcher, it is useful for me to read the paper and understand which algorithm I would like to use for my task of interest.\n\nI also really liked that there is a three bullet point summary outlining recommendations and outcomes of the study that was easy to follow.\n\nI think the stability metric is really the differentiating factor across different algorithms for comparison, so I think that is a useful contribution.\n\nWeaknesses:\n\nI do not think the stability metric is consistent across datasets, and it is difficult to compare use in each. Please see the “Requested Changes” section for more information.\n\nThe 3D bar charts are difficult to interpret, and somewhat redundant. While in the tables, up and down arrows indicate “good” or “bad” performance, these are absent on the bar charts. Generally, since the graphs are projected in 3D, comparison within and across figures is difficult.\n\nThe water data set is randomly split into 1k/10k/90k random structures. However, is random sampling not “difficult” enough for testing? Is there a way to first cluster structures/configurations, and hold out those that are distant from those in the training set?\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The authors compare state-of-art molecular simulation algorithms across a number of different benchmark datasets and evaluations metrics. \n\nThe different datasets vary from small to large molecules in relevant tasks. The metrics are those that are used to train models and for typical comparison, as well as new metrics that the authors have created to help compare methods.\n\nThe authors are careful about performance of each, and speculate why models will be performing in certain ways.\n\nThe code will also be available on github for others to run and compare against.\n\n", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "While I think the “Stability” metric is useful, I feel that it is inconsistently used across datasets. For some datasets, it is, to my understanding, ensuring that the current state of the simulation does not differ from some set of states in the past. For other datasets, (i.e. macromolecules), it uses bond angles to ensure realistic simulation conditions. If the metrics mathematically parameterize different things, they should be different metrics.\n\nI do not think the title is justified by the work of this paper: “Forces are not enough”. Did the paper specifically address this question? Moreover, weren’t there ML models that did use forces that were successful on many desired criteria? While catchy, I do not think the paper thoroughly and concretely outlines this statement.\n\nI found Figure 1 to be confusing and out of place. Only after reading the whole paper and understanding the metrics did I realize what this figure was trying to tell me, because I needed to understand if larger or smaller metric values were “good”. Ideally, the first figure should condense the authors’ approach and relevant data (condensing Figures 2 and 3). It also takes quite a bit of reading to understand what “water-10k” means out of context.\n\nWhile I realize the authors are comparing a number of different methods, 3D bar charts are notoriously difficult to interpret. I’m not sure they add much more versus colored heatmap tables.\n\nNice-to-have changes:\n\nFigures:\n\nIt would be great to have a figure that described the different classes of models: SE(3)-equivariant, E(3)-invariant, and translation-invariant. As a reader, I have to read each of these papers to understand how their capacity differs.\n\nIf you are tight on real estate for the paper, I think figures 2 and 3 could be condensed into a single figure.\n", "broader_impact_concerns": "Not Applicable", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "RofND3axh2", "reviewer": "Reviewer_nktP", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Overall I think this is a good paper. The key contributions are impactful, relevant, and generally well-executed. To me the main weaknesses are related to presentation. Many things in the paper could be explained or argued better, but something that was particularly lacking was the motivation and justification of the metrics. Without this, I think the key claims of the paper are not fully supported since they depend on the metrics. Details below.\n\n**Strengths**\n\n- this kind of benchmark is the _good kind of obvious_: it did not occur to me that ML models should be evaluated using MD, but when I read the abstract it felt immediately obvious to me that this is better than MAE on forces. Given that most papers _do not_ use MD for  evaluation, such a contribution is very welcome.\n- the main conclusion that force MAE $\\neq$ MD performance is also _obvious (in a good way)_: when exploring the state spaces edge cases matter, and average performance does not highlight edge cases. However, to my knowledge this intuitive result was not empirically demonstrated before on a wide range of models. The results in this paper not only show that this can happen in practice, but that the mismatch is _systemic_ across many models. I think Figure 4 demonstrates this very well\n- benchmark is well-executed: admittedly I think I lack expertise to be a great judge of this, but the choice of datasets, evaluation conditions, and metrics seems well-thought out and thorough, even if the rationale for some choices is unclear (see weaknesses below)\n\n**Weaknesses**\n\n- Metrics not discussed sufficiently: $h(r)$ is mentioned only in an indirect way when discussing the MD17 dataset, and the stability metric is only described qualitatively. I thought this was odd since the main contribution of this paper is the benchmark, which is defined by the datasets, MD setup, and metrics. The conclusions also depend on the metrics, so for a reader to believe the claims it is critical to show or argue that these metrics are the right ones. Not only are the choice of metrics not discussed much (e.g. which alternatives did you choose from), _the definitions are not even in the main paper_!\n- Motivation and discussion of stability is particularly lacking.\n  - Not only is it not defined in the main text (see point above), the authors only describe it in a vague way with phrases like \"may not extrapolate robustly to the undersampled configuration space\" and \"trajectories can enter nonphysical states that are not meaningful for observable calculations\". I did not know what \"nonphysical\" and \"meaningful\" specifically meant in this context. I think it needs a much more precise description.\n  - In the appendix the authors more or less define \"stability\" as \"large deviations in RDF metric\", using a seemingly arbitrary threshold to define \"large\", which differs for each datasets. How were these thresholds chosen? I felt disappointed that this was put into the appendix like a minor design choice when it is actually a critical part of the benchmark which is essential to support the authors' claims.\n  - Outside of section 4, the authors refer to stability as if it were a measure of when the simulation \"fails\" (indeed the name \"stability\" seems to imply this), when _that is not how it was defined_. I presume that when the model outputs \"unstable\" trajectories running the simulation for longer does not help, but as far as I can tell this was never actually stated or shown by the authors. Without this, I felt it was inappropriate to discuss stability in this way, and overall I was not convinced of its utility as a metric.\n  - Because of these things, I decided to consider all claims stemming from stability to be \"unsupported\", which I marked in the yes/no question below. In reality I believe these claims are probably true, but I think the issues above need to be fixed so that the claims are rigorously justified.\n- Other MD design choices are not discussed: e.g. why were specific choices for temperatures and sampling times made? This is a key contribution of your paper so it should at least be stated why certain design choices were made.\n- The MD performance metrics to some degree involve comparing the performance of the ML FF trajectories with those from another method. However, it is not stated very clearly what these methods are, and it is not argued why they are suitable for use as reference data. For example, if the reference data was another ML FF method then arguably it may be unreliable. I'm guessing that the choice of reference data was good, I just think it should be explained and argued to convince a skeptical reader.\n- Discussion of related work could be improved. It seemed like the most relevant works which use force fields for MD are dismissed saying \"they didn't compare multiple models\". While this is indeed true, the stated contributions of this paper are not only _running_ the benchmark on multiple models, but _designing_ the experiments themselves, so I think such a dismissal is inappropriate. The experiments of many of these papers actually seem quite similar to those of this benchmark. I picked Zhai et al 2022 (at random) and saw that they also do MD on water and report a variety of observables, comparing the ML FF to a classical one and to experimental values. The key difference here seems to be which observables are used and other MD parameters. I think the authors either need to explain more clearly how their experimental setup is different/better than those or other works, _or_ state more clearly that the experimental procedure is fairly standard/in line with prior work and therefore discuss this work less dismissively. Remember that this is TMLR: if aspects of the work are standard or not very novel, the authors do not need to pretend otherwise to get accepted.\n- Scope of models considered could be improved: the authors seem to focus exclusively on the \"latest and greatest\" neural-network based models in their evaluation. While I think including these models is good, \"machine learning\" encompasses much more than just neural networks, for example body-order decompositions with learned coefficients. It would be nice to see at least one more classical model included, even something as simple as a fitting a Lennard-Jones potential model. This would give the reader a sense of whether neural networks are actually improving over other methods in any meaningful way.\n- Rhetoric could be improved. In general, I feel like the authors _undersold_ the importance of their benchmark. I think there are a lot of good arguments for the use of an MD-based benchmark which the authors did not fully exploit. I've given some suggestions in the \"requested changes\" section.\n- Better background information on MD: a reader unfamiliar with the key ideas of MD would not be able to get much out of this paper. MD is not a very complicated idea but people with an ML/math background might not be very familiar with it. TMLR papers should not assume a physics background for the readers, so I think it would help if the authors expanded this.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The authors create a benchmark which assesses models that outputs force fields (FFs) by the accuracy of molecular dynamics (MD) simulations using those FFs. The benchmark consists of a curated set of datasets, MD setups, and evaluation metrics. They evaluate a number of existing models using their benchmark. Their main conclusion is that the models whose force predictions have a lowest mean absolute error (MAE) are _not_ in general the models with the best MD performance, suggesting that FF MAE is not a good proxy for usefulness in MD.\n\nI would list the contributions and new knowledge of this paper as the following:\n\n1. The datasets, MD settings, and evaluations metrics which define their benchmark\n2. The evaluation of many popular ML models on their benchmark\n3. Analysis of the results from point 2 to demonstrate empirically that lower FF MAE does not consistently result in better MD (at least using their metrics)\n\nConversely, one thing which is definitely **not** a contribution of this paper is the _idea_ of evaluating FF models using MD. Most papers using FFs are motivated by MD (see e.g. the review by Unke et al), and evaluating models based on their ability to do this is an obvious conclusion. I say this only because the authors do not seem to make an explicit claim about whether this idea is novel (it is not in their contribution list, but I also couldn't find a sentence stating that it is not a novel idea).", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "I hope the authors can make the following changes which I think are critical to ensure that the paper's main claims are well-supported:\n\n- Move definitions of metrics into the main text since they are important.\n- Add some sort of justification of why these metrics make sense. For example, why the RDF and not other observables? Why is stability important?\n- Reconcile the mismatch between the description of stability and its mathematical definition, in addition to arguing why it is an important metric. This could take many forms depending on what the authors actually observed during \"unstable\" MD simulations. For example, if the authors observed that once the trajectory became \"unstable\" it would never recover, it would be good to show this empirically, and also show that this can be detected by changes in the RDF. The authors should also think about finding a more descriptive name than \"stability\" if the trajectories do not in fact diverge or fail, since I do not think that the inclusion of non-physical states in a trajectory makes it \"unstable\"\n- Add at least a minimal discussion of the way the MD simulations are set up. For example, were times and time steps chosen to give each experiment a similar computational cost? Were the temperatures chosen to be in line with the training data? Currently this is completely unclear. I think it is important to show the reader that you made reasonable design choices when constructing the benchmark.\n\nI also suggest that the authors consider the following changes which I believe would strengthen the work, even though I would be willing to accept the paper without them. Given that the paper is currently _below the 12 page limit_, I think there is plenty of room for extra content.\n\n1. Argue the importance of your benchmark more directly in your introduction. The introduction seemed to more or less say \"people do MD, some people want to learn force fields for MD, so we make a benchmark for this.\" I think it should instead say \"MD simulations are _actually the main reason why people work on ML force fields_, so performance on MD is not just a nice benefit, it is a critical requirement!\" I think this would convince a reader that this is an important way to measure performance, and also strengthen the impact of your main conclusion: optimizing FF MAE may actually be _misleading_ the field to research models which aren't actually useful for anything! (although if you make these changes be sure to acknowledge that MD is not the only thing that people use FFs for, another one is relaxing structures to the nearest energy minimum).\n2. Provide a better introduction to MD, which at the bare minimum should contain the key MD equation $d^2x/dt^2=-m^{-1} dE/dx$ (i.e. Newton's equation). I suggest moving section 3 to right after the introduction, or maybe even merging it with the introduction. I would also state that MD is a well-established field with tons of prior work (maybe citing a few textbooks), and add some references to your discussion of observables which currently has none. You could also state some basic challenges of MD such as accuracy/time tradeoff for modelling E(x), accuracy/speed tradeoff for the time step, challenges scaling MD to systems with many atoms, etc. This would also be a good time to state, at a high level, where ML models make different design decisions than classical methods. For example, many classical methods try to approximate the energy in a sparse way (e.g. only including atoms which are bonded together, or in a local region). My understanding is that many ML methods make slightly different approximations, which has implications for both scaling and accuracy, while other ML methods are much more similar to classical methods but just replace fixed functions with neural networks. Overall, I think these changes would help convince readers unfamiliar with MD that your benchmark makes sense and that research in this area is important.\n3. Give the reader an intuitive understanding of why a low MAE when producing forces will not necessarily imply good performance in MD. I think this is an obvious possibility to anybody with physics knowledge but may not be obvious to those without. My suggestion for doing this: choose a 2-body system in 1D Lennard-Jones potential which can be visualized in a nice 1D plot. Show two approximations of the potential: one which has errors of moderate magnitude everywhere, and one which matches very closely almost everywhere except in one location where it has an erroneous, but very deep second minimum. The MAE of the first one will be larger, but the MD of the second one will be less realistic because the trajectories could get stuck in the second erroneous minimum.\n4. Expand discussion of most relevant related work as discussed above.\n5. Add some classical or semi-classical methods as baselines. These should be quick to implement and run if they are sufficiently simple.\n6. Small thing: I suggest changing the colouring of the tables to not be all shades of green, but instead use a color map with various shades so that it is easier to differentiate performance visually. Personally I found it difficult to differentiate the various shades of green. Be careful to choose a color map which is friendly to color blind people though.", "broader_impact_concerns": "No concerns about ethics/broader impact.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper studies machine learning force fields commonly used in materials science. While previous comprehensive benchmarks have only evaluated the mean force errors of such force fields, the current paper went beyond this evaluation and analyzed the resulting stability of molecular dynamics simulations. 

The authors successfully addressed initial concerns during the rebuttal period. These changes included a stronger motivation of the stability metrics and their connections to the earlier literature, a restructuring of the introduction and related work sections, additional experiments on Aspirin molecules, as well as more details on dataset generation and the choice of benchmark models. Extensive changes were carried out during the review period. It is worth pointing out that the authors and reviewers worked collaboratively and constructively on further improving the manuscript. 
Decision: Accept (poster)