The dark side of the forces: assessing non-conservative force models for atomistic machine learning
Authors: Filippo Bigi
Marcel F. Langer
Michele Ceriotti
Venue: ICML 2025 oral (2025)
Keywords: geometric machine learning
energy conservation
atomistic modelling
molecular dynamics
statistical mechanics
Primary area: applications->chemistry_physics_and_earth_sciences
TL;DR: An assessment of non-energy-conserving geometric machine learning models for atomic-scale systems

Abstract: The use of machine learning to estimate the energy of a group of atoms, and the forces that drive them to more stable configurations, have revolutionized the fields of computational chemistry and materials discovery.
In this domain, rigorous enforcement of symmetry and conservation laws has traditionally been considered essential. For this reason, interatomic forces are usually computed as the derivatives of the potential energy, ensuring energy conservation. 
Several recent works have questioned this physically constrained approach, suggesting that directly predicting the forces yields a better trade-off between accuracy and computational efficiency -- and that energy conservation can be learned during training.
This work investigates the applicability of such non-conservative models in microscopic simulations. 
We identify and demonstrate several fundamental issues, from ill-defined convergence of geometry optimization to instability in various types of molecular dynamics.
Contrary to the case of rotational symmetry, energy conservation is hard to learn, monitor, and correct for.
The best approach to exploit the acceleration afforded by direct force prediction might be to use it in tandem with a conservative model, reducing -- rather than eliminating -- the additional cost of backpropagation, but avoiding the pathological behavior associated with non-conservative forces.

The dark side of the forces: assessing non-conservative force models for atomistic machine learning
Filippo Bigi" Marcel F. Langer'
Michele Ceriotti
Abstract
1. Introduction
Interatomic potentials model the microscopic interactions between atoms, and determine directly by modulating the response to thermal excitations the stability and reac- tivity of molecules and materials. Over many decades, inter- atomic potentials have been used in Monte Carlo (Metropo- lis et al., 1953) and molecular dynamics (MD) (Andersen, 1980) simulations, geometry optimization, and other tech- niques, allowing mechanistic study of the atomic-scale behavior and properties of molecules, materials, and biolog- ical systems (Allen & Tildesley, 2017; Tuckerman, 2008). While traditional interatomic potentials are based on simple physically inspired functional forms, the last decade has seen machine-learned interatomic potentials (MLIPs) ob- tain remarkable accuracies at a high level of computational efficiency, most often by learning reference energies and forces from much slower first-principle quantum mechanical calculations (Behler, 2021; Unke et al., 2021).
The use of machine learning to estimate the en- ergy of a group of atoms, and the forces that drive them to more stable configurations, has revolu- tionized the fields of computational chemistry and materials discovery. In this domain, rigorous en- forcement of symmetry and conservation laws has traditionally been considered essential. For this reason, interatomic forces are usually computed the derivatives of the potential energy, ensur- ing energy conservation. Several recent works have questioned this physically constrained ap- proach, suggesting that directly predicting the forces yields a better trade-off between accuracy and computational efficiency, and that energy con- servation can be learned during training. This work investigates the applicability of such non- conservative models in microscopic simulations. We identify and demonstrate several fundamen- tal issues, from ill-defined convergence of geom- ctry optimization to instability in various types of molecular dynamics. Given the difficulty in monitoring and correcting the lack of energy con- servation, direct forces should be used with great care. We show that the best approach to exploit the acceleration they afford is to use them in con- junction with conservative forces. A model can be pre-trained efficiently on direct forces, then fine-tuned using backpropagation. At evaluation time, both force types can be used together to avoid unphysical effects while still benefitting al- most entirely from the computational efficiency of direct forces.
At first, such machine-learned interatomic potentials were trained on quantum mechanical data specific to the chem- ical system of interest. Over the last years, more diverse datasets have emerged, consisting of up to billions of tar- gets, and spanning much of the periodic table (Chanussot et al., 2021; Tran al., 2023; Barroso-Luque Luque al., 2024; Schmidt et al., 2024). This has not only led to increase in the complexity of the models, going from linear or kernel regression to deep graph neural networks, but it has also pushed many recent models to abandon some of the under- lying physical symmetries of interatomic potentials in favor of simpler and potentially more scalable and efficient architectures (Pozdnyakov & Ceriotti, 2023; Neumann et al., 2024; Qu & Krishnapriyan, 2024). In such models, symme- tries are learned from data, aided by data augmentation.
Some recent models for interatomic potentials also disregard the property of energy conservation (Gasteiger et al., 2021; Neumann et al., 2024; Liao et al., 2023). While conservative models calculate interatomic forces as the derivatives of total energy with respect to the positions of the atoms, conservvvveee models predict them directly, therefore breaking this constraint. Although this practice can lead to more computationally efficient neural networks, the use of such models in practical atomistic simulations has not yet
'Laboratory of Computational Science and Modeling, IMX, École Polytechnique Fédérale de Lausanne, 1015 Lau- sanne, Switzerland. Correspondence to: Michele Ceriotti <michele.ceriee cecbhetleccriotttie.. @pff.cch>.
Proceedings of the 4nnd International Conference Machine Learning. Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).

---

The dark side of the forces
been studied in detail.
emerged, with various goals (sampling different thermody- namic ensembles (Andersen, 1980), improving sampling efficiency, observing rare events (Laio & Parrinello, 2002), or accounting for nuclear quantum effects (Chandler & Wolynes, 1981), etc.). In this work, we focus on some of the simplest variants:
This work compares conservative and non-conservative machine-learned interatomic potentials. Following brief review of the most common types of potentials and their applications, we discuss the theoretical implications of us- ing non-conservative forces to drive atomistic simulations, highlighting potential shortcomings. We then demonstrate the impact of these conceptual problems in several case stud- ies, highlighting the pitfalls of using non-conservative inter- atomic forces in production, and demonstrating how, con- trary to the case of geometric symmetry breaking, it is diffi- cult to monitor and correct the impact of non-conservative models. Instead, it might be preferable to supplement con- servative model with direct force predictions, using them to accelerate simulations that have well-defined, energy- conserving forces the ground truth.
MD in the NVE ensemble. The number of atoms N. volume V and total energy (kinetic plus potential) are kept con- stant, describing the behavior of an isolated system. Usually, this is accomplished by using the velocity Verlet integra- tor (Verict, 1967) with a short time step At.
MD in the NVT ensemble. Here, the goal is to simulate system at constant temperature T. This achieved by modifying the dynamics with so-called thermostats. Local thermostats apply velocity corrections on each atom inde- pendently (Bussi & Parrinello, 2007). While they correctly sample the desired constant-temperann be used to evaluate static, equilibrium properties, this comes at the cost of disrupting the dynamics of the system. Global thermostats (Bussi et al., 2007) address this shortcoming by modifying the atomic velocities in a concerted manner, making it possible to investigate time-dependent proper- ties (Bussi & Parrinello, 2008).
ensemble and can
2. Background and related work
2.1. Interatomic potentials and atomistic simulations
An interatomic potential V describes the potential energy between the N atoms of a molecule or material as a function of their positions (r. J and chemical nature {a.]
V(rr, M))
(L)
Many other applications of interatomic potentials (such as lattice dynamics, phonon calculations, etc.) require the cal- culation of forces and their higher-order derivatives. How- ever, for simplicity, this work focuses on the most common use cases, discussed above.
While some applications, such as Monte Carlo simulations or energy difference calculations (for transition state de- fect energies, material stability, etc.) need only the values of V, the most popular uses of interatomic potentials require the evaluation of interatomic forces, defined as the negative derivatives of V with respect to the atomic positions:
2.2. Physical symmetries
Interatomic potentials obey a number of physical symme- tries and constraints.
f = OV (fr,, 1 Ōr,
The most notable applications which make use of inter- atomic forces are:
E(3)-invariance. Interatomic potentials are invariant un- der the transformations of the Euclidean group in three dimensions E(3), including translations, rotations, and re- flections. Mathematically, given group element g E E(33) that acts on all positions r,
Geometry optimization. This technique consists of find- ing one more minima of the potential energy surface V to identify the preferred structure of microscopic system. Both local and global optimization algorithms are employed, and the vast majority use the forces fj both during optimiza- tion and as a stopping criterion. Similar methods are used to identify saddle points of V, which are associated with the energy barriers that determine the rate of chemical reac- tions (Henkelman et al., 2000).
V(lg rs, at (n a((
Permutation invariance. Interatomic potentials are in- variant with respect to permutations of atom indices, i.c..
V(.., rs, 0i rj. Fj. a.. -)- V(.,, (..., ry. -r. r.,...... (4)
Molecular dynamics. MD aims at simulating the behav- ior of microscopic system by solving its classical equations of motion numerically. Using a time step At, the simplest forms of molecular dynamics propagate a discretized ver- sion of Hamilton's 's equations (Verlet, 1967). Over more than 50 years, many variants of molecular dynamics have
Energy conservation. Interatomic forces are conservative, i.c., their mechanical work W over closed loop is zero:
N f dr, 0.
W
(5)
1=1
2

---

The dark side of the forces
This is direct consequence of (2). The term conserva- tive alludes to the fact that an isolated physical system in which only conservative forces act obeys the principle of conservation of mechanical energy.
with the same asymptotic computational cost as the en- ergy prediction) with automatic differentiation (Griewank & Walther, 2008). Nevertheless, differentiation incurs a computational overhead, typically 2-3x for inference and 3x for training; the exact theoretical factors are discussed further in Appendix B. This overhead can be avoided by directly predicting forces during the forward pass. How- ever, by removing the relationship between forces and en- ergy defined in Eq. (2), such models do not enforce energy conservation. The possibility of performing a direct, non- conservative, evaluation of the forces was realized early (Li et al., 2015), but used only sparsely until recently, when it was applied to equivariant graph neural networks such as GemNet (Gasteiger et al., 2021) and several more recent universal machine learning interatomic potentials including ORB (Neumann al., 2024) and Equiformer (Liao et al., 2023), which are trained on large datasets and promise broad applicability for chemical modeling. Some of these architec- tures also perform direct evaluation of stresses, which leads to the violation of conservation of enthalpy (see App. G.3).
2.3. Interatomic potentials via graph neural networks
For many years, most machine learning interatomic poten- tials made use of feature- based models (Behler & Parrinello, 2007; Bartók et al., 2010) trained on energies and forces from quantum mechanical calculations. This combination of physically inspired descriptors (Musil et al., 2021; Langer et al., 2022) and classical machine learning methods (linear and kernel models, as well as shallow neural networks) de- livers good accuracy and computational efficiency on small datasets, usually created for a specific application in chem- istry and/or materials science.
However, recent efforts to build much larger and more com- prehensive datasets (Deng et al., 2023; Schmidt et al., 2024; Eastman et al., 2023; Tran et al., 2023; Chanussot et al., 2021), have made it apparent that graph neural networks exhibit superior accuracy and scalability compared to earlier models (Batzner et al., 2022), as they respect the permu- tation symmetry of interatomic potentials while delivering rich and flexible representations of the atomic structures. For a review of graph neural networks, we redirect the inter- ested reader to Zhou et al. (2020), while their applications to interatomic potentials is detailed in Duval et al. (2023).
The impact of lack of energy conservation for practical simulations and property prediction has recently been in- vestigated in related works: Eissler et al. (2025) probe the limits of unconstrained architectures and find that lack of energy conservation becomes more pronounced for larger target systems. Loew et al. (2024) find that direct-force mod- els perform poorly at the prediction of phonon properties, which require second derivatives of the potential energy sur- face. Fu et al. (2025) observe a better correlation between test set error and downstream predictive performance, for instance for phonons, for models that can perform stable MD simulations, i.e., that are able to conserve energy.
2.4. Breaking physical constraints
Several recent architectures do not enforce all physical con- straints listed in 2.2. aiming to increase expressivity, ability to scale to large datasets, and computational efficiency. Effi- cient inference is crucial to enable the length and time scales required for practical molecular dynamics simulations, as forces must be evaluated at every time step.
2.5. Use of forces for training
Machine learning models of V are typically trained jointly on energy and force labels. The relative weight of these labels in training, or in the extreme case, whether to train exclusively on energy or on forces, has been discussed exten- sively, both for interatomic potentials and diffusion models for atomistic systems (Chmiela et al., 2018: Christensen & von Lilienfeld, 2020; Wang et al., 2024; Ren et al., 2024).
Rotationally unconstrained models. In the case of rota- tional symmetry this approach, combined with training- time rotational augmentation, has been shown to yield accurate and efficient models (Pozdnyakov & Ceriotti, 2023; Neu- mann et al., 2024: Qu & Krishnapriyan, 2024) with neg- ligible, or casily controllable, impact on physical observ- ables (Langer et al., 2024).
One important consideration is that focusing on energies or forces emphasizes different components of a potential en- ergy surface, with consequences that are not directly visible when inspecting training and validation accuracies. To a first approximation, the distortions observed in NVT MD simulations - a common strategy to build training sets can be interpreted as collection of quasi-harmonic oscillators of different frequency w, with high- greeffreyency modes associ- ated with short-range covalent bonding, and low-frequency ones associated with collective motions, which are often the most relevant for applications. The statistical mechanics of
Non-conservative models. The conservative property for forces in Eq. (5) is satisfied if and only if there exists function (usually named the potential energy) of which the forces are the spatial derivatives, see Eq. (2). Therefore, a machine learning model that calculates interatomic forces as derivatives of the potential energy according to Eq. (2) will be trivially conservative. The required derivatives, the gradient of the scalar V, can be obtained efficiently (i.e.,
3

---

The dark side of the forces
harmonic oscillator imply that (V) x 1, while (f2) x w² (see Appendix A). In other terms, the largest contributions to the forces from a thermally sampled dataset come from high- frequency modes that are hard to integrate and may lead to instabilities in the dynamics, while the contribution from the slow modes, which are hard to sample, but usually the most relevant, are under-emphasized. The potential energy provides a more balanced representation of the different molecular time and length scales.
conservative component. À can also be computed only for entries associated with an atom i, pair of atoms ij, providing a finer-grained assessment of the violation of Eq. 6, as seen in Figure 1.
This local symmetry breaking can also be measured by inte- grating the work done by the forces over a closed loop, that is bound to be zero (within integration error) for conserva- tive force field (Eq. 5). An explicit test of non-conservative behavior can also be implemented by monitoring the total energy in an NVE MD simulation, or equivalently con- served quantity that keeps track of the heat flux associated with the thermostat (Bussi & Parrinello, 2007) in an NVT simulation. To find the power P (energy per unit time) in- jected by the non-conservative forces during a molecular dynamics trajectory, it is sufficient calculate the average rate of change in the conserved quantity C' during a section of the trajectory. P AC/AI.
The question of relative energy and force weights during training is related to, but distinct from, the question of learn- ing energy and forces consistently, i.e., enforcing that the predicted forces are exactly the gradient of the predicted energy. Models targeting mean-field energies, as in coarse- graining (Wang al., 2019) and centroid MD (Musil al., 2022), usually train exclusively estimates of the mean forces, but use conservative formulation as these forces are defined as derivatives of a thermodynamic potential that is difficult to evaluate explicitly. Disregarding the connec- tion with reference energies may also be beneficial when the two targets are subtly inconsistent, because of conver- gence issues, the use of heterogeneous calculation settings, or numerical techniques like Fermi smearing (Marzari et al., 1999), which can be required for metallic systems. This may contribute to the empirical observation that non-consistent training is preferable in datasets like OC20 (Chanussot et al., 2021) and OC22 (Tran et al., 2023) that contain a high frac- tion of metallic configurations.
3.2. Side-effects of direct force prediction
Predicting forces the derivatives of a translationally in- variant potential ensures that the total net force acting on all atoms is zero, and that for a potential that is rotation- ally invariant the torque on an isolated molecule is zero. A direct force model - irrespective of whether it is E(3) invari- ant - does not have the same guarantees. These spurious effects are easy to remedy, by subtracting the total force and torque from each prediction. This technique is used by ORB (Neumann et al., 2024), and we also adopt it here.
3. Theory
There is another non-trivial consequence of using a direct- force prediction architecture. Whereas the potential en- ergy is usually estimated as the sum of atomic contributions but is a global property, forces are atom-centered. When predicting them derivatives, many atomic environments contribute the force f; on each atom. When predicting directly, only the i-atom centered environment contributes to f,. Hence, direct force models can be expected to be affected more directly by the geometric degeneracies of low- body-order atom-centered descriptors (Pozdnyakov et al., 2020) and do not benefit from the same extended effec- tive interaction range as conservative forces (Artrith et al., 201 1). We discuss these effects in Appendix C, presenting some empirical evidence that direct force models require larger range to match the force accuracy of comparable conservative model.
Some of the consequences of using non-conservative models descend directly from their mathematical formulation and can be used to foresee their impact on typical applications.
3.1. Measuring non-conservative behavior
If a vector field f is the derivative of a smooth scalar function V, then its Jacobian J (the Hessian of V) contains the mixed second derivatives of V. and must therefore be symmetric:
Hb.a Jıa,59 0r arga
JJ3,ta
In order to quantitatively capture the amount of non- conservation in specific force prediction from trained model, it is then possible to compare the Frobenius norm (or any other matrix norm) of the antisymmetric component of the Jacobian to that of the Jacobian itself:
3.3. Effects on geometry optimization
In order to assess the stability of materials or molecules at low temperature, a common approximation is to search for minimum-energy configurations. This can be achieved by minimizing the potential energy V as a function of the atomic positions with most of the widely used algorithms relying (in some cases exclusively) on the value of the gradi-
x |Jantilie
||]||
where Janti = (J -J-J X then defines a metric going from 0 for conservative forces to for forces that have
4

---

The dark side of the forces
ent. The lack of a consistent potential energy is problematic for most optimization schemes: those which require an plicit evaluation of the objective function (e.g., to perform line searches) cannot be used; those that just "follow the forces", relying on the vanishing magnitude of the force as a stopping criterion, can fail because non-conservative forces can keep driving indefinitely in the same direction, e.g., following closed loops with negative total work.
An alternative approach is to train both a conservative and non-conservative predictor and adding force-matching loss term, or simply training both on the same forces labels. We will discuss this approach in further detail in Section 4.8. It is important to note that, in both strategies, energy conser- vation can be trained both on labeled and unlabeled data.
4. Results
3.4. Effects on molecular dynamics
Having introduced the subject of non-conservattie force fields and discussed the potential pitfalls that might be in- curred when using them in practice, we will now examine their effect on a range of applications, using liquid water as the main, paradigmatic example.
is not uncommon to observe violations of energy conser- vation in MD simulations, since finite-timestep integrators violate the exact correspondence between the change in kinetic and potential energy along a trajectory. This is usu- ally accepted, because (1) in well-designed simulations this leads to small fluctuations, and not to run-away effect; (2) the notion of a shadow Hamiltonian (Hairer et al., 2006) ensures that simulations reach a steady state that is "sta- tistically close" that generated by an exact integrator: (3) thermostatting techniques can relatively easily control small integration errors, so that structural and dynamical observables are not affected significantly (Morrone et al., 1). The fact that no underlying Hamiltonian can be defined for the dynamics generated by a non-conservativy force field suggests that, in this case, artifacts might be more pronounced and harder to correct. For example, the sym- plectic nature of Hamiltonian dynamics is no longer valid (this is easy to see, for example, from Eq. 5.2, Chapter 6 of Hairer et al. (2006)), and the theorem of equipartition of energy (whose proof is also based on the existence of a Hamiltonian (Pathria, 2017)) does not apply. As we shall see, this is what we observe empirically.
4.1. The models
In order to substantiate our empirical observations, we per- form our experiments on multiple models. Our main ex- amples rely on the rotationally unconstrained PET architec- ture (Pozdnyakov & Ceriotti, 2023), trained from scratch the bulk water dataset of Ref. 22 using both conser- vative (PET) and non-conservative (PET-NC) architecture. Additionally, we train "PET-M" to predict both conservative and non-conservative forces. To assess the implications of a direct prediction of forces in the most favorable possi- ble context, we primarily use the best-performing. custom- trained models. We also show some results for the non- conservative ORB-v2 model (Neumann et al., 2024), which is currently state-of-the-ar for several materials prediction benchmarks. Even though ORB is not trained on this spe- cific dataset, and is therefore at a clear disadvantage, it provides an indication of the relevance of the issues we discuss. In the appendices, we also discuss several other architectures, including a "legacy" SOAP-BPNN architec- ture, as well pre-trained foundation models, MACE-MP- 0 (Batatia et al., 2023), SevenNet (Park et al., 2024), and EquiformerV2 (Liao et al., 2023), which we apply to few diverse materials in Appendix G. A table describing all employed models, along with full details on the different architectures, can be found in Appendix D.
3.5. Learning conservative behavior
The standard approach to making symmetries more eas- ily learnable from data is to employ data augmentation at training time: A random element of the underlying sym- metry group, for instance rotations, is selected and applied for every sample mini-batch. This approach has been successfully employed in the domain of MLIPs, as well as in range of other applications including computer vi- sion (Quiroga et al., 2020). However, it is only applicable to explicit geometric symmetries, and therefore not to en- ergy conservation, which not a symmetry with respect to inputs, but rather a symmetry of derivatives, as discussed in Section 3.1. Nevertheless, we briefly consider different schemes to promote energy conservation during training.
4.2. Accuracy
In terms of sheer accuracy, see Table 1, our tests confirm that forces provide very useful information to train an in- teratomic potential, in particular for a dataset containing relatively large configurations. Using forces in the training of conservative model dramatically improves the accu- racy of energy predictions with just minor degradation in the accuracy for f with respect to a model trained only on forces. A non-conservative force model exhibits about 30% higher force error than a conservative architecture, and including a separate energy head leads to lower error than a model trained just on indicating that the sharing of
One approach is to include the measure of Jacobian symme- try X, Equation (7), as a term in the loss function. This faces severe practical issues: With automatic differentiation, com- puting J requires multiple (3NV in the absence of sparsity or stochastic approximations) evaluations of the potential.
5

---

The dark side of the forces
weights within the architecture is beneficial. We also show results for the PET-M hybrid architecture, which makes both conservative and non-conservative force predictions. Its accuracy is less than 10% worse than the best models for either architectures. As we will discuss in Section 4.8, this is an excellent way to exploit con-conservvtive forces simulations. Before doing so, however, we will assess the behavior of purely con-conservative models.
102
10
10a
10a
10c
10
// evev
10
10
10-2 2
10
Jv
10-3
1g 1/112
10
Ag
4.3. Non-conservative behavior
4
10
10
The asymmetry the Jacobian the most direct, point- wise measure of non-conservative behavior. Different non- conservative models show widely different values of X 0.015 for ORB, 0.017 for Equiformer, 0.032 for SOAP- BPNN-NC and 0.004 for PET-NC, computed on few wa- ter structures from the test set of Cheng et al. (2019). As we shall see, the magnitude of À correlates qualitatively with the stability of models in simulations. The sym- metry of J applies separately to each pair of atoms, and so it possible to extract further insights by computing the norm of the antisymmetric part of each block Jij function
Fij A
Figure I. Comparison of the norm of each block of the Jacobian, Jij. and of its antisymmetric component, for different pairs of atoms, as function of their distance Tij. computed for randomly selected bulk water structure from the set.
Non-conservative behavior can also be demonstrated by computing the work along closed path (see Appendix E). Given that the choice of the path is arbitrary, think it is more relevant to quantify the practical implications of this issue in terms of an energy drift in molecular simulations.
solved for different atomic pairs and plotted of the interatomic distance (Figure 1). The asymmetry is also present for "on-site" blocks, i.e., swapping only the Cartesian coordinates used the derivatives for a given atom; the relative magnitude Jasy small compared to the magnitude of the Jacobian; the asymmetric component between atoms i and decays with the interatomic distance more slowly than the absolute magnitude of J (which has been used as a measure of the interactions between pairs of atoms (Herbold & Behler, 2022)), and in the intermedi- ate regime around 6Â it becomes comparable size so that the pair-resolved Jacobian asymmetry Xi approaches I for large interatomic distances. This latter observation has important implications when applying these models in simulations, the impact of non-conservative behavior different atomic-scale processes is not uniform, and it tends to be larger- - in a relative sense - for collective processes involving long-range correlations.
4.4. Constant-energy molecular dynamics
consider the of non-conservative MLIPss in the context of MD simulations. Also in this case, we focus our analysis on the best-performing models for liquid water, PET and PET-NC. A more thorough comparison, including several foundation models and a few homogeneous and het- erogeneous material structures, is discussed Appendix G, make
Let
and consistently corroborates the observations here. constant-pressut simulations, which, if performed with direct-stress models, break conservation of enthalpy and therefore lead to drift in the volume of the simulation, are shown in App. G.3 using the general-purpose PET-MAD potential (Mazitov et al., 2025).
Given that non-conservative models lack well-defined con- served quantity by construction, rely indirect measure- ments of the sampled ensemble. We consider in particular the kinetic temperature T 2K 2K/(3Nkpp where N is the number of particles considered. This is just rescaling of the kinetic energy K; its ensemble average should corre- target temperature for NVT trajectories (300 in these tests), and to value in its vicinity for NVE trajecto- ries initialized from a thermally equilibrated configuration. As sensitive indicator of the dynamical behavior of the system, we compute the Fourier transform of the velocity- velocity correlation function, 6ve (w). Its peaks are closely related to the density vibrational modes and to infrared and Raman spectra, and its w - 0 limit is proportional to the diffusion coefficient.
Table 1. Test errors (energies in meV per atom, forces in meV/A) of PET models trained the bulk water dataset. In this table, and all subsequent tables, rows corresponding to non-conservative models are highlighted with gray background.
spond
ARCH.
TYPE
TRAINING
MAE(V)
MAE(F)
PET
V
4.7
PET
C
f
18.6
PET
C
V.F
0.55
19.4
PET
NC
f
24.3
PET
NC
V.
1.42
24.8
C
20.2
PET-M
V,f
0.59
NC
26.7
6

---

The dark side of the forces
800
ORB. NVE PET. NVE PET-NC, NVE PET-M, NVE MTS-8
Table 2. Mean kinetic temperature excess (AT) - (T) - 300 K. and atom-type resolved temperatures (ATo) and (ATit). for ORB, the conservative (C), non-conservative (NC) and multiple time- stepping (MTS) PET models (which evaluates conservative forces once every 8 steps). White-noise Langevin (WN) and stochastic velocity rescaling (SVR) thermostats are also compared.
700
600
500
400
THRM.
TYPE
r/fs
(T)/K
(Ta)/K
(To)/K
ORB
300
WN
NC
1000
51.0(6)
60.4(5)
33(1)
WN
NC
100
4.2(2)
5.9(3)
0.9(3)
10000
20000
30000
40000
WN
NC
10
0.4(1)
0.6(1)
0.1(1)
step
SVR
NC
10
1.0(1)
36.2(8)
-70(2)
Figure 2. Time series for the kinetic temperature along a NVE MD trajectory, for ORB, the conservative and non-conservative PET models, and for PET-M model using a multiple time-stepping (MTS) algorithm that evaluates conservative forces every steps.
PET
WN
C
100
0.1(2)
0.0(2)
0.3(3)
WN
NC
1000
12.8(5)
11.2(7)
16.2(5)
WN
NC
100
1.4(2)
1.3(2)
1.6(3)
WN
NC
10
0.1(1)
0.0(1)
0.2(1)
SVR
C
10
0.1(1)
-0.4(3)
1.0(7)
The failure of the non non-conservative model in NVE dynamics apparent in Figure 2. Whereas for a conservative potential the kinetic temperature fluctuates around the initial value, the spurious work associated with non-conservative forces leads to large drift T:: To put it on a human scale, this un- physical drift corresponds to rate of heating of about 7'000 billion degrees per second for the custom-trained PET-NC model, and another 10 times larger for the general-purpose ORB model. This spurious energy flow is a clear signature of consennsevvttive behavior, and makes direct-force mod- els entirely useless for constant-energy simulations. This runaway increase of the kinetic energy can be mitigated whenever auxiliary model is available to evaluate the potential energy - by adjusting the velocities to enforce en- engy conservation artificially (see Appendix G.4). Similar to what we will discuss for constant-temperraurre simulations, the trajectories are still affected by large artefacts.
SVR
NC
10
0.3(1)
-4.4(3)
9.9(6)
SVR
M-8
10
0.0(1)
-0.1(4)
0.1(9)
feres with atomic motion) are needed, even at 100 fs the average temperature is significantly off the target value. We discuss how these deviations in the equilibrium temperature affect structural properties of water in Appendix G.5. The upshot is that for accurate models and strong thermostatting, the effects are small but noticeable. Furthermore, the strong Langevin thermostatting reduces the sampling efficiency, and so even with respectable I ns trajectory, many simple structural averages are not fully converged.
4.6. Sampling efficiency and time-dependent properties
Aggressive Langevin dynamics is bound to dramatically change time-dependent properties, and in particular to re- duce the diffusion coefficient - so that longer trajectories are needed to collect statistically independent atomic con- figurations. This slow-down is apparent when looking at the velocity correlation spectra (Figure 3). In the weak coupling regime (WN, 1000 fs) there is a (small) increase in diffusion coefficient relative the reference, because of the unphysically higher temperature, while the high-frequency peaks corresponding to stretching and bending are only weakly perturbed. Stronger couplings alter the dynamics dramatically, and reduce the diffusion coefficient (and hence the efficiency in sampling slow, collective motion) by fac- tor of about 1.5 (7 100 fs) and 5 (r 10 fs), negat- ing the inference speed-up of non-conservative model while making it impossible to accurately evaluate any time- dependent property.
4.5. Equilibrium properties in the NVT ensemble
The very use of a finite time step in the integration of MD trajectories causes energy fluctuations, and it is not uncom- mon use advanced simulation schemes that violate ergy conservation (Kühne et al., 2007; Mazzola & Sorella, 2017; Morrone et al., 201 I; Laio & Parrinello, 2002), for instance because they use approximations that yield forces contaminated by a stochastic noise. In these cases, using a thermostat can counterbalance the energy error, and obtain stable trajectories that yield configurations and equilibrium average properties close the correct NVT ensemble despite the drift of the conserved quantity that generalizes the total simulations.
energy for constant-temperaate
Judging by the average temperature (T') (Table 2), it is rela- tively easy to control the non-conservative behavior using a white-noise (WN) Langevin thermostat. However, strong couplings (the time scale over which the thermostat inter-
A potential solution applied often in similar cases, includ- ing to control the artifacts non-invariant predictions of V (Langer al., 2024) is to resort a global thermostat
7

---

The dark side of the forces
10
400
PET-C PET-C, fine-tuned PET-NC
Val. force MAE MAE/ /mev/A
(w((w) /arb. units
200
10-1
100
10-2
C, SVR, T= fs
50
NC. WN. T=
NC. WN. 1=100
NC. WN, T= fs
10-3
NC, SVR, T=10 fs
20
10a
102
10a
10
20
100
200
/ cm
GPU hours
Figure 3. Velocity power spectrum (w), ). for different PET mod- els and thermostat types: conservative (C) and concooservvatie (NC) model using white-noise Langevin (WN) and stochastic ve- locity rescaling (SVR) thermostats.
Figure 4. Training curves for PET-C and PET-NC models, and the conservative fine-tuning of hybrid PET-M model initialized (epoch marked with arrow) from the potential energy head of the PET-NC model.
on gradients, to avoid the complication of using inconsis- tent energies, as discussed in Section 3.3, comparing the FIRE (Bitzek al., 2006) algorithm that similar in spirit to zero-temperature MD and LBFGS (Liu & Nocedal, 1989) quasi-Newton algorithm that uses approxima- tion of the Hessian to accelerate convergence. Comparing different models on the task of optimizing a water snapshot from a MD configuration (Figure 7) shows that inaccurate non-conservative models, such as SOAP-BPNN-NC. fail catastrophically at geometry optimization, while more accu- rate models, such PET-NC and ORB, can reach locally stable configuration, especially using FIRE. We note, how- ever, that non-conservative models are less stable when used with a basssia--based method, with large fluctuations in the residual force that make hard to define a stopping criterion. On practical level, non-conservative forces are bound to make geometry optimization more fragile, and to require careful choice of the minimization algorithm and its con- vergence parameters, as observe in Appendix F when comparing different general-purpose models.
that only acts on the total kinetic energy rather than on in- dividual particle momenta, achieving efficient temperature control without disrupting dynamics. We use the stochas- tic velocity rescaling (SVR) method (Bussi et al., 2007), which indeed brings the average temperature to within 1% of the target, without dramatically altering (w) even with a strong = 10 fs coupling. However, a global thermostat cannot help when con-conservative terms act differently on the various degrees of freedom. This is evident in how the temperature O and H atoms, computed separately, devi- ates by up 10% from the target, which is reflected in loss of structure measured by g(r), and in overestimation of the diffusion coefficient.
Some of the more sophisticated thermostats used to stabilize other types MD approximationss such as carefully tuned generalized Langevin equations (Ceriotti et al., 2009; Mor- rone et al., 2011) can also be used to enforce more aggres- sive temperature control with reduced dynamical disruption in conservative molecular dynamics, but they still mod- ify the natural dynamical properties significantly, and they can fail catastrophically when used with con-conservative forces, as shown in Appendix G.6).
4.8. Non-conservative forces accelerators
While have observed that conservative MLIPs are bet- ter suited for practical simulations, we suggest that hybrid models, which additionally support direct, non-conservative, force predictions, can be used for faster inference and train- ing. Such models can be obtained by training both force heads jointl,, as demonstrated the PET-M model, or, more efficiently, by first training non-conservative model and then fine-tuning its energy head to yield conservative forces. As shown in Figure 4 and further discussed in Appendix H, conservative fine tuning leads to the accuracy and physi- cal correctness of conservative models at highly reduced training time.
These experiments show clearly that while is possible to mitigate the runaway temperature increase associated with the lack of energy conservation, doing so in way that does not disrupt structural and/or dynamical observables is highly nontrivial or even impossible.
4.7. Geometry optimization
MD performed at very low temperature can be regarded as form of geometry optimization. Our observations from MD suggest that sufficiently accurate non non-conservative models should be able to reach reasonable, low-energy structures. We restrict ourselves to optimization algorithms based only
In simulations, one can then use the conservative forces of
8

---

The dark side of the forces
such a hybrid model for validation, error monitoring and correction, and the direct forces for faster inference. A good example is to use multiple time-stepping (MTS) tech- niques (Tuckerman et al., 1992) for molecular dynamics, where the non-conservative forces are used to integrate the equations of motion, and the conservative forces are applied every M steps as a correction. This reduces the theoretical overhead of a conservative trajectory from a factor of F x 2 to one of 1 1+ (F-1)) 1)/M. The results using this technique in Table 2, Figure 2 and in Appendix I are essentially indistin- guishable from fully conservative ones, using M = 8, which leads to a small, approximately 20% slowdown compared to a direct-force, non-conservative trajectory. Appendix contains further explanation of the MTS technique, as well as more detailed results for MTS simulations using models trained on the water and OC20 datasets. The technique can be also successfully used for constant- preesurt-peesure simulations, as shown in Appendix G.3.
validation set accuracy, but are less reliable in describing the slow, collective structural rearrangements that are often the key drivers of the most relevant microscopic processes,
Given that the target forces are conservative, accurate mod- els usually exhibit less pronounced non-conservative behav- ior. As a consequence, one can expect that, as the field moves to larger training datasets and more expressive mod- els, some of the pathological effects we observe will become less severe. Our findings, however, suggest that the best way to exploit the speed-up afforded by direct prediction of the forces is not to replace conservative models, but to aug- ment them with a non-conservative head. This can be used to accelerate training by first training a con-conservative model and then fine-tuning its energy head to yield accurate conservative forces through differentiation. The resulting "multi-force" models can also be used to speed up many different types of simulations, by alternating conservative and non-conservative evaluations, avoiding the narrower applicability and inherent instability associated with relying exclusively on direct force predictions. This insight enables the efficient training of the next generation of universal machine-leaming interatomic potentials while retaining the physical correctness required for practical simulations.
5. Discussion
Chemical and materials modeling at the forefront of de- velopment in the applications of machine learning sci- ence. The field has long been advocating for the use of physically informed architectural constraints, but there are indications that its bitter-lesson moment is coming, with the realization that deploying physics-agnostic models at scale provides better outcomes than exploiting physical priors. It appears that this is the case for some of the geometric sym- metry constraints (Pozdnyakov & Ceriotti, 2023; Neumann et al., 2024; Qu & Krishnapriyan, 2024), and growing number of architectures disregard the physical connection between the interatomic potential and the corresponding forces (Gasteiger et al., 2021; Neumann et al., 2024; Liao et al., 2023). With respect to this latter constraint, our study paints a nuanced picture. Atomistic simulations rely on the assumption that forces are the exact derivatives of the potential, and small deviations from this constraints lead to instabilities. Non-conservative behavior also results in molecular dynamics trajectories exhibiting a spontaneous drift away from the desired thermodynamic conditions. Con- trolling this effect with thermostats requires careful tuning, and disrupts both time-dependent properties and the sam- pling efficiency of the trajectory, negating the computational advantage of direct-force architecture.
Software and Data
Code and data required to reproduce the re- sults in this work are available on Zenodo at wwwpodd.orgg/eeeeee wwtg:record//177889909 An example of how to perform multiple-time-step dynamics with conservative and direct forces can be found at https://atomistee cockbookcoooook..oooooo org/ wwwllle//pettetee naatt-mmmdncp/pe-mmadade html, and an example of conservative fine-tuning at
cookbook org/ exampless intttttfiettttnnn//eeeeeee foat--ff-ree ft-ne. html. More details on the software used are available in Ap- pendix J.
kttps://atomist
Acknowledgements
The authors would like to thank Federico Grasselli and Niklas Schmitz for stimulating discussion, and Rafael Gomez-Bombarelli for a coffee-break conversation which inspired us to look into this problem. ML and MC acknowl- edge funding from the European Research Council (ERC) under the European Union's 's Horizon 2020 research and innovation programme Grant No. 101001 890-FIAMMA. FB and MC acknowledge support from the NCCR MAR- VEL, funded by the Swiss National Science Foundation (SNSF, grant number 182892) and from the Swiss Platform for Advanced Scientific Computing (PASC).
Contrary to the case of rotational symmetry that is easy to monitor and correct at inference time (Langer et al., 2024), and learn through data augmentation, assessing non- conservative behavior requires the explicit evaluation of the Jacobian both as diagnostics and as additional loss term. Fur- thermore, energy and forces are complementary targets, and disregarding the former may lead to potentials that appear stable, resilient to dataset inconsistencies, and with good
9

---

The dark side of the forces
Impact Statement
Kozinsky, B. E(3)-qquivarian graph neural networks for data-efficient and accurate interatomic potentials. Nature Communications, 13(1):2453, May 2022. ISSN 2041- 1723. doi: 10. 100880/444400 467-022-29939-55
This paper presents work whose goal is to advance the ap- plication of machine learning to simulations of molecules and materials, which can be used for many purposes. As potential impacts hard to predict and wide-ranging, we do not believe it is necessary to discuss them here in detail.
Behler, J. Four Generations of nighh-Dimensional Neural Network Potentials. Chemical Reviews, 121(16):10037- 10072, August 2021. ISSN 0009-2665, 1520-6890. doi: 10. 10.102//asscchemree0000883
References
Behler, J. and Parrinello, M. Generalized Neural-Network Representation of High-Dimensional Potential-Energy Surfaces. Physical Review Letters, 98144)144440, April 2007. ISSN 0031-9007. doi: 10.1 103/PhysRevl.en.99 146401.
Allen, M. P. and Tildesley, D. J. Computer Simulation of Liquids, volume 1. Oxford University Press, Novem- doi: 10. 1093/0//////
ber 2017. ISBN
9780198803195000 011.118888811150010000000.
Andersen, H. C. Molecular dynamics simulations at con- stant pressure and/or temperature. The Journal of Chem- ical Physics, 72(4):2384-233)) 1980. ISSN 00219606. doi: 10. 10.1003/1.00 .439486.
Bigi, F., Chong, S., Kristiadi, A., and Ceriotti, M. Flashmd: long-stride, universal prediction of molecular dynamics. arXiv preprint arXiv:2505. 19350, 2025.
Bitzek, E., Koskinen, P., Gähler, F., Moseler, M., and Gumb- sch, P. Structural Relaxation Made Simple. Physical Re- view Letters, 1(11):170201 October 2006. ISSN 0031- 9007, 1079-7114. doi: 10.1 10.1105Phhyy
Artrith, N., Morawietz, T., and Behler, J. High-dimensional neural-network potentials for multicomponent systems: Applications zinc oxide. Physical Review 83(15):
153101, April 2011. ISSN 1098-0121. doi: 10.1 DtssRevB.88.1131000
Bussi, G. and Parrinello, M. Accurate sampling using Langevin dynamics. Physical Review E. 75(5):56707, 2007. doi: 10.1 1..10//PhysReee..5.5500000
Barroso-Luque, L.,, Shuaibi, M., Fu, X., Wood, B. M., Dzamba, M., Gao, M., Rizvi, A., Zitnick, C. L.., and Ulissi, Z. W. Open materials 2024 (omat24) inor- ganic materials dataset and models. arXiv preprint 2rXiv:2410.12771 2024.
Bussi, G. and Parrinello, M. Stochastic thermostats: Com- parison of local and global schemes. Computer Physics Communications, 179(1-3):2--29, July 2008. ISSN 00104655. doi: 10. w0.10066j..pc22008.00 00.000.
Bartók, A. P., Payne, M. C., Kondor, R., and Csányi, G. Gaussian Approximation Potentials: The Accuracy of Quantum Mechanics, without the Electrons. Physical Review Letters, 104(13)) April 2010. ISSN 0031- 9007. doi: 10.1100/hysReeeeett.. 1...114....444. 136403.
Bussi, G., Donadio, D., and Parrinello, M. Canonical sam- pling through velocity rescaling. Journal of Chemical Physics, 126(1):1410) 2007.
Bussi, G., Zykova- Zykova-Timan, T., and Parrinello, M. Isothermal- isobaric molecular dynamics using stochastic velocity rescaling. The Journal of Chemical Physics, 130(7): 074101, February 2009. ISSN 0021-9606. doi: 10.1063/ 1. 1.3073889.
Bartók, A. P., Kondor, R., and Csányi, G. On repre- senting chemical environments. Physical Review B, 87(18):184115 May 2013. ISSN 1098-0121. doi: 10.1 1..100//hyssevvB.S0 1....8e2.......11....
Batatia, I., Kovacs, D. P., Simm, N. C., Ortner, C., and Csanyi, G. MACE: Higher order equivariant message passing neural networks for fast and accurate force fields.
Ceriotti, M., Bussi, G., and Parrinello, M. Langevin equa- tion with colored noise for constant-temperature molec- ular dynamics simulations. Physical Review Letters, 102(2):020601, January 2009. ISSN 00319007. doi: 10. 10.1103/hhsRReveet..o0 102.000000.......0000.
In Oh, A. H., Agarwal, A., Belgrave, D., and Cho, K. (eds.). Advances in Neural Information Processing Sys- tems, 2022.
Ceriotti, M., Bussi, G., and Parrinello, M. Colored-noise thermostats à la Carte. Journal of Chemical Theory and Computation, 170):1170-1180,, April 2010. ISSN 15499626. doi: 10. 10.1001/c1900566
Batatia, I., Benner, P., Chiang, Y., Elena, A. M., Kovács, D. P., Riebesell, J., Advincula, X. R., Asta, M., Avaylon, M., Baldwin, W. J., et al. A foundation model for atomistic materials chemistry, arXiv preprint arXiv:2401.0m 00090,0096.0000 2023.
Chandler, D. and Wolynes, P. G. Exploiting the isomor- phism between quantum theory and classical statistical mechanics polyatomic fluids. Journal of Chemical Physics, 74(7):4078 4075, 1981.
Batzner, S., Musaelian, A., Sun, L., Geiger, M., Mailoa, P. Kornbluth, M., Molinari, N., Smidt, T. E,, and
10

---

The dark side of the forces
Chanussot, L.,, Das, A., Goyal, S., Lavril, T., Shuaibi, M., Riviere, M., Tran, K., (leras-Domings J., Ho, C.,, Hu, W., et al. Open catalyst 2020 (oc20) dataset and community challenges. Acs Catalysis, ((11))6059-6077 2021.
Griewank, A. and Walther, A. Evaluating Derivatives.
Society for Industrial and Applied Mathematics, Jan- uary 2008. ISBN 77808987165977 9780898717761. doi: 10. 10.1197/119780089777776666 URL https: /doi.
Cheng, B., Engel, E. A., Behler, J., Dellago, C., and Ceriotti, M. Ab initio thermodynamics of liquid and solid water. Proceedings of the National Academy of Sciences of the United States of America, 116(4): 10-1115, January 2019. ISSN 10916490. doi: 10. ww..as../pns...a.. 18111///p1a...11117116
Hairer, E.,, Hochbruck, M., Iserles, A., and Lubich, C. Geo- metric numerical integration. Oberwolfach Reports, 3(1): 805-882. 2006.
Henkelman, G., Uberuaga, B. P., and Jónsson, H. climb- ing image nudged elastic band method for finding saddle points and minimum energy paths. Journal of Chemical Physics, (1(2))99901, 2000.
Chmiela, S., Sauceda, H. E., Müller, K.-R., and Tkatchenko, A. Towards exact molecular dynamics simulations with machine-learned force fields. Nature Communicationss 9: 3887, 2018. doi: 10.1038/s41467001 8-0662--.
Herbold, M. and Behler, J. Hessian-bass assessment of atomic forces for training machine learning interatomic potentials. The Journal of Chemical Physics, 156(11): 114106, March 2022. ISSN 0021-9606, 1089-7690. doi: 2.11063/5000880
Christensen, A. S. and von Lilienfeld, O. A. On the role of gradients for machine learning of molecular energies and forces. Machine Learning: Science and Technology, 1(4): 045018, 2020. doi: 10.
Hjorth Larsen, A., Jørgen Mortensen, J., Blomqvist, J., Castelli, I. E., Christensen, R., Dutak, M., Friis, J., Groves. M. N., Hammer, B., Hargus, C., Hermes, E. D., Jennings, P. C., Bjerre Jensen, P., Kermode, J., Kitchin, J. R., Leonhard Kolsbjerg, E., Kubal, J., Kaasbjerg, K., Lysgaard, S., Bergmann Maronsson, J., Maxson, T., Olsen, T., Pastewka, I.., Peterson, A., Rostgaard, C., Schiøtz, J., Schütt, O., Strange, M., Thygesen, Vegge, T. Vilhelmsen, L., Walter, M., Zeng, Z., and Jacobsen, K. W. The atomic simulation environment a Python library for working with atoms. Journal of Physics: Condensed Malter, 29(27):273002, July 2017. ISSN 0953-8984, 1361-648X. doi: 10. (088/1361-6488) aa680c.
Deng, B., Zhong, P., Jun, K., Riebesell, J., Han, K., Bartel, C. J., and Ceder, G. Chgnet as a pretrained universal neural network potential for charge-informed atomistic modelling. Nature Machine Intelligence, 5(9):1031-104) 2023.
Duval, A., Mathis, S. V., Joshi, C. K., Schmidt, V., Miret, S., Malliaros, F. D., Cohen, T., Lio, P., Bengio, Y., and Bronstein, M. A hitchhiker's guide to geometric gnns for 3d atomic systems. arXiv preprint a:Xiv2311..0511, 2023.
S.,
Eastman, P., Behara, P. K., Dotson, D. L., Galvelis, R., Herr, E., Horton, T., Mao, Y.. Chodera, D., Pritchard, B. P., Wang, Y., De Fabritiis, G., and Markland, T. E. SPICE, A Dataset of Drug-like Molecules and Peptides for Training Machine Learning Potentials. Scientific Data, 10(1):11, January 2023. ISSN 2052-4463. doi: 10. $41597-002-011882--
Kapil, V., Vandelon J., and Ceriotti, M. Accu- rate molecular dynamics and nuclear quantum effects at low cost by multiple steps in real and imaginary time: Using density functional theory to accelerate wave- function methods. Journal of Chemical Physics, 144 (5)005111,, February 2016. ISSN 00219606. doi: 10. 14491114
Eissler, M., Korjakow, T., Ganscha, S., Unke, O. T., Müller, K.-R., and Gugler, S. How simple can you go? an off-the- shelf transformer approach to molecular dynamies, 2025. URL httph:tp/////. www.aaiiiivc.coo org/ abs 25033. 01431.
Kühne, T. D., Krack, M., Mohamed. F. R., and Parrinello, M. Efficient and Accurate Sar-Parrinello-like Approach Born- Oppe-Oppeneeimer Molecular Dynamics. Physical Review Letlers, 98(6)66401, February 2007. ISSN 0031- 9007. doi: 10. 28.11003Physseeeeett......66
Fu, X., Wood, B. M., Barroso-Luque, L., Levine, D. S., Gao, M., Dzamba, M., and Zimnick, C. I.. Learning smooth and expressive interatomic potentials for physical property prediction, 2025. URL httptt/s:://aaa //arxivv///xvvva//// orgoorgggooo abs/ 2502. 12147.
Laio, A. and Parrinello, M. Escaping free-energy minima. Proceedings of the National Academy of Sciences, 99(20): 12562-12566. 2002.
Gasteiger, J., Becker, F., and Günnemann, S. Gemnet: Uni- versal directional graph neural networks for molecules. Advances in Neural Information Processing Systems, 34: 6790-6802, 2021.
Langer, M. F., Goeßmann, A., and Rupp. M. Represen- tations of molecules and materials for interpolation of quantum-mechanice simulations via machine learning.

---

The dark side of the forces
npj Computational Materials, 8(1):41, December 2022. ISSN 2057-3960. doi: 10. 10.1038//4115444022--772122
Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H., and Teller, E. Equation of State Calcula- tions by Fast Computing Machines. Journal of Chemical Physics, 21(6):1087-1092 1953.
Langer, M. F., Pozdnyakov, S. N., and Ceriotti, M. Prob- ing the effects of broken symmetries in machine learn- ing. Machine Learning: Science and Technology, 5 (4):04LT01, December 2024. ISSN 2632-2153. doi: 10. 00.1088/2632-2155/akke
Morrone, J. A., Markland, T. E., Ceriotti, M., and Berne, B. J. Efficient multiple time scale molecular dynam- ics: Using colored noise thermostats to stabilize reso- nances. Journal of Chemical Physics, 134(():14103, Jan- uary 2011. ISSN 00219606. doi: 10.1063/1.3518366
Li, Z., Kermode, J. R., and De Vita, A. Molecular dynamics with on-the-fly machine learning
teantum-mechane forces. Physical Review Letters, 1 1444),, March 2015. ISSN 10797114. doi: 10.1 10.1100/hyysRevLtt.... 096405.
Musil, F., Grisafi, A., Bartók, A. P., Ortner, C.,, Csányi, G., and Ceriotti, M. Physics-Inspired Structural Representa- tions for Molecules and Materials. Chemical Reviews, 121 (16):9759-98) August 2021. ISSN 0009-2665, 1520- 6890. doi: 10.1021/acs. wwwmrev.1c0000
Liao, Y.-I.,, Wood, B., Das, A., and Smidt, TT. Equiformer22: Improved equivariant transformer for scaling to higher- degree representations. arXiv preprint 2.fi::200.120555 2023.
Musil, F., Zaporozhets, I., Noé, F., Clementi, C., and Kapil, V. Quantum dynamics using path integral coarse- graining. The Journal of Chemical Physics. 157(18): 181 102, November 2022. ISSN 0021-9606, 1089-7690. doi: 10.1063/5.0120386
Litman, Y., Kapil, V., Feldman, Y. M. Y., Tisi, D., Begušićć T., Fidanyan, K., Fraux, G., Higer, J., Kellner, M., Li, T. E., P6s, E. S., Stocco, E., Trenins, G., Hirshberg. B., Rossi, M., and Ceriotti, M. I-PI 3.0: A flexible and efficient framework for advanced atomistic simulations. The Journal of Chemical Physics, (6)00224, August 2024. ISSN 0021-9606, 1089-7690. doi: 10.1063/5. 0215869.
Neumann, M., Gin, J., Rhodes, B., Bennett, S., Li, Z., Choubisa, H., Hussey, A., and Godwin, J. Orb: A fast, scalable neural network potential. arXiv preprint arXi::4410.2257. 2024.
Nigam, J., Pozdnyakov, S., Fraux, G., and Ceriotti, M. Uni- fied theory of atom-centered representations and message- passing machine-learning schemes. The Journal of Chem- ical Physics, 156(20):204115, May 2022. ISSN 0021- 9606, 1089-7690. doi: 10.1063/5.008704
Liu, D. C. and Nocedal, J. On the limited memory BFGS method for large scale optimization. Mathematical Pro- gramming, 45(1-3):503-5273 528, August 1989 ISSN 0025- 5610, 1436-4646.
doi: 10. 1007//FF018889166
Park, Y., Kim, J., Hwang, S., and Han, S. Scalable parallel algorithm for graph neural network interatomic potentials in molecular dynamics simulations. Journal of Chemical Theory and Computation, 2024.
Loew, A., Sun, D., Wang, H.-C., Botti, S., and Marques, M. A. L. Universal machine learning interatomic potentials are ready for phonons, 2024. URL https::////aa.. /arviv. 16551.
Pathria, R. K. Statistical Mechanics: International Series of
Monographs in Natural Philosophy, volume 45. Elsevier, 2017.
Marzari, N., Vanderbilt, D., De Vita, A., and Payne, M. C. Thermal Contraction and Disordering of the AI(1 Surface. Physical Review Letters, 82(16):3296- 3299, April 1999. ISSN 0031-9007, 1079-7114. doi: 10.1 10.100hyyssee..eee 82.3296.
Pozdnyakov, S. and Ceriotti, M. Smooth, exact rotational symmetrization for deep learning point clouds. In Advances in Neural Information Processing Systems, vol- ume 36, pp. 79469-79501. Curran Associates, Inc., 2023.
Mazitov, A., Bigi, F., Kellner, M., Pegolo, P., Tisi, D., Fraux, G., Pozdnyakov, S., Loche, P., and Ceriotti, M. PET- MAD, universal interatomic potential for advanced ma- terials modeling. arXiv preprint arXiv:2503. 14118, 2025.
Pozdnyakov, S. N., Willatt, M. J., Bartók, A. P., Ortner, C., Csányi, G., and Ceriotti, M. Incompleteness of Atomic Structure Representations. Physical Review Letters, 125: 166001, 2020. doi: 10. 1a..1o/Phyyseeeettt.........
Mazzola, G. and Sorella, S. Accelerating ab initio Molecular Dynamics and Probing the Weak Dispersive Forces in Dense Liquid Hydrogen. Physical Review Letters, 118 (1):015703, January 2017. ISSN 0031-9007, 1079-7114. doi: 10. 103/PhysRevL.ett. 18.005
Qu, E. and Krishnapriyan, A. S. The importance of being scalable: Improving the speed and accuracy of neural network interatomic potentials across chemical domains. arXiv preprint arXiv:2410. 20.24/62.99... 2024.
12

---

The dark side of the forces
Quiroga, F., Ronchetti, F., Lanzarini, L.,, and Baríviera, A. F. Revisiting data augmentation for rotational invariance in convolutional neural networks. In Modelling and Sim- ulation in Management Sciences: Proceedings of the International Conference on Modelling and Simulation in Management Sciences (MS-18), pp. 127-141. Springer, 2020.
Wang, J., Olsson, S., Wehmeyer, C., Pérez, A., Charron, N. E., De Fabritiis, G., Noé, F., and Clementi, C. Machine Learning of Coarse-Grained Molecular Dynamies Force Fields. ACS Central Science, 5(5):755-761 May 2019. ISSN 2374-7943, 2374-7951. doi: 10.1021 8b00913.
Wang. Y.. Wang, L., Shen, Y. Wang, Y., Yuan, H., Wu, Y., and Gu, Q. Protein conformation genera- tion via force-guided SE(3) diffusion models. In Pro- ceedings of the 41st International Conference on Ma- chine Learning, pp. 56835-56859, PMLR, 2024, doi: notavailable. URL https: / /proceedings mmlr. press/ /23/s/v235//// wttp///ng22.....o /Fess//2255wwag22vvv..
Ramakrishnan, R., Dral, P., Rupp. M., and Von Lilienfeld, O. A. Quantum chemistry structures and properties of 134 kilo molecules. Scientific data, 1())11-7,2044
Ren, Y., Zheng, D., Liu, C., Jin, P., Shi, Y., Huang. L.., He, J., Luo, S., Qin, T., and Liu, T.-Y. Physical con- sistency bridges heterogeneous data in molecular multi- task learning. In Advances in Neural Information Pro- cessing Systems 38 (NeurIPS 2024), Vancouver, Canada, Dec 10-Dec 15, pp. not available. not available, 2024. doi: notavailable. URL https: setpps:/openreeeee http://pppeeeeeeeeeeeee n forum? 1i=:dddtaGGFFeeee
html.
Zhou, J., Cui, G., Hu, S., Zhang, Z., Yang, C., Liu, Z., Wang, L., Li, C., and Sun, M. Graph neural networks: A review of methods and applications. Al open, 1:57-81, 2020.
Schmidt, J., Cerqueira, T. F., Romero, A. H., Loew, A., Jäger, F,, Wang. H.-C., Botti, S., and Marques, M. A. Improving machine-learning models in materials science through large datasets. Materials Today Physics, 48: 101560, 2024.
Tran, R., Lan, J., Shuaibi, M., Wood, B. M., Goyal, S., Das, A., Heras-Domingg J., Kolluru, A., Rizvi, A., Shoghi, N., et al. The open catalyst 2022 (oc22) dataset and challenges for oxide electrocatalysts, ACS Catalysis, 13 (5):3066-3084, 2023.
Tuckerman, M. Statistical Mechanics and Molecular Simu- lations. Oxford University Press, 2008
Tuckerman, M., Berne, B. J., and Martyna, G. Reversible multiple time scale molecular dynamics. The Journal of Chemical Physics, 97(3): 1990, 1992. ISSN 00219606. doi: N0.1063//.463133
Unke, O. T., Chmiela, S., Sauceda, H. E., Gastegger, M., Poltavsky, I., Schütt, K. T., Tkatchenko, A., and Müller, K.-R. Machine learning force fields. Chemical Reviews, 121(16):10142-1111 10186, 2021. doi: 10.1021/acs.chemrev 0c01111.
Verlet, L. Computer "Experiments" on Classical Flu- ids. Thermodynamical Properties of Lennard-Jones Molecules. Physical Review, -59, July 1967. ISSN 0031-899X. doi: 10.1103/hh:Ree.155..9
Villar, S., Hogg, D. W., Storey-Fisher, K., Yao, W., and Blum-Smith, B. Scalars are universal: Equivariant ma- chine learning, structured like classical physics. Advances in Neural Information Processing Systems, 34:28848- 28863, 2021. doi: 10. 20..4880///xx..216...66
13

---

The dark side of the forces
A. Forces and potentials as training targets
Including both forces and the potential energy V as targets when training a ML potential (either separately, or jointly for a conservative model) requires choosing a weighting factor to combine the errors into single loss function. Independently from the relative weight, is interesting consider how these two targets affect the accuracy of model for different types of molecular displacements. To investigate this aspect, we can approximate material in the vicinity of stable structure, i.e., local minimum of the potential energy V. as a quadratic form, which can be written in the basis of the eigenvectors of the (mass-scaled) Hessian to give a superimposition of harmonic modes,
1 - Σ mws Te: V(g)= = Σ Vil(n)) -
(
where m is the atomic mass (the expression generalizes the case of multiple atomic types), wi are normal mode frequencies and q): the displacements along the eigenvectors.
If we now consider how configurations are sampled at a constant temperature T (running short MD trajectories is a common strategy to generate training sets for MLIPs), one sees that each harmonic mode is distributed as p(gx) oc exp(-Vi(gs)/kBB)) = exp(-mu@gl /2kB)) As a consequence, one can easily compute the expectation values
TTT w
(Vi(g))) x TT,, ( x Tu
These textbook results highlight the following facts: (1) fww-freuencyy normal modes are those associated with the largest structural deformations - and therefore, often, with phase transitions and important molecular rearrangements; (2) thermal excitations affect all normal modes equally in terms of potential energy contributions; (3) the largest force contributions come from high-frequency dlw-displaccmmnt) vibrations.
Thus, when training on forces using an L loss, molecular modes associated with high-frequency vibrations are over- emphasized. For example, if one considers the water dataset we use in this work (Cheng et al., 2019), the total force acting on each water molecule has root mean square of 0.97 ev/,, while the residual "ooeraamo forces (that predominantly short-ranged and associated with high-frequency molecular vibrations) are four times larger, 3.93 eV/Å. This very crude analysis highlights the non-trivial implications of using forces as (direct or indirect) training targets.
B. Theoretical computational cost
Conservative forces are generally computed from potential energies by backward propagation of gradients. The vast majority of operations in neural networks (and nearly all those that take up significant computational time) are binary operations which can be expressed as the computation of f(x, y) starting from and y. During the backward step corresponding to such operation, /V/ and 0V/Oy must be found from 0V/0f. In the case of matrix multiplication, the forward calculation of [z, y). the backward calculation of 0v/O= 0x and the backward calculation of 0V/0y each consist of a matrix multiplication with the same computational cost. Since matrix multiplications can be assumed to be the most costly components of neural networks, this generally implies that backward gradient computations are around twice as expensive as the corresponding forward function evaluation.
However, in the case of backward force evaluation, operations where is an internal representation and y is weight can save the /Vy calculation. In simple multi-layer perceptron, where all linear layers correspond to this type of computation, this would yield a backward propagation that is roughly as expensive as the forward pass. This is not the case in transformers, as the attention mechanism involves comparatively expensive operations where both x and y are internal representations, and one can expect the backward propagation of gradients to be somewhere between 1 x and 2x the cost of the forward evaluation.
C. Range of pack-propagated and direct force models
The use of cutoff to restrict the range of interactions is ubiquitous in the construction of physics-based potentials, and is also an integral part of descriptor- ML potentials (Musil et al., 2021). It is often argued that models that incorporate correlations between at least two neighbors of each central atom achieve an effective interaction range of twice the cutoff distance (Artrith et al., 2011) 1) (see Figure 5a). A similar effect also applies to message- -aassggg architectures (Nigam et al., 2022).
14

---

The dark side of the forces
b
k
cut
T cut
Figure 5. schematic representation of the interactions a (ML) interatomic potential. (a) Contributions an i-centered prediction of energy or force from threc or more neighbors, labeled j and k. (b) Contributions from a j-centered energy prediction to the force on the i-th atom, computed by pack-propagation.
One can see clearly that this extension of the interaction range beyond the cutoff (or the receptive radius of the GNN) does not apply to the case in which forces are predicted directly. Considering for simplicity the case three-body model in which the potential contribution from each atomic environment can be written sum of a function of its distances with two neighbors j and k and the distance between the neighbors Tyk
V1= = Σ J.k
urigg Tik, Tyk);
(10)
one sees that the dependency on interatomic distances greater than the cutoff due the relative position of atoms other than the central atom. Thus, direct model that predicts a similar 3-body force
f
Σ 3.k
f(rig, Tak, Tyk)s
Tak,
(1I)
any other functional form limited the neighbors of the i-th atom, only contains information atoms within the cutoff. The dependency of f, on the coordinates of far-away atoms is consequence of the fact that the total energy is built as sum over multiple centers. It is only through the terms of the form 0V, that occur naturally when evaluating forces through back-propagation, than the force depends on the position of the neighbor's ss neighbors. More generally, in a message-passing implementation, backpropagation ensures that force predictions benefit from an effective receptive radius that is twice that of atom-centered energy predictions.
A further concern for direct force models is that least for low-body-order models atom-centered descriptors can be shown to have low resolving power, with pairs of distinct atomic environments having precisely the same representa- tion (Pozdnyakov et al., 2020). For all known degeneracies, descriptors centered on other atoms allows distinguishing the structure as whole, and as consequence the total energy and interatomic forces can be still differentiated. This would not be the case for direct force predictions, which would fail completely to differentiate degenerate pairs when using atom-centered low-body-order models, and would be more sensitive to numerical instabilities for higher-order models.
In practice we find consistent evidence of the practical impact of these considerations. As an example, Table reports the accuracy of PET-C and PET-NC models, using 2 and 3 message passing layers. It can be seen that the accuracy of the direct force model benefits much more from the increase in the receptive radius.
Model
2 message- peessage-passing layers
message-passingg layers
PET-C PET-NC
20.8
18.6
32.8
24.3
Table 3. Test-set force MAE, in units of me V/A, for conservative and non-conservative force-only PET models using and 3 message- passing layers, respectively. Upon increasing the receptive radius, the accuracy of the model improves by around 30%, while the accuracy of the C model only improves by around 10%.
15

---

The dark side of the forces
D. Models and architectures
Even though main text is focused on custom-trained models, and emphasize the most accurate non-conservative model we could obtain, we also want to provide an overview of the behavior of less-performant custom-trainee model, and of several publicly available general purppss models. A comprehensive list of the models we consider is reported in Table 4.
Model
Description
ORB-v2
Non-equivariant, non-conservative model, trained on the Alexandria and MPtrj
datasets.
Equiformer
Equivariant, non-conservative model, trained on the Alexandria and MPuj datasets.
MACE-MP-0
Equivariant, conservative model, trained the MPujj dataset.
SevenNet
Equivariant, conservative model, trained the MPtrj dataset. A re-implementation of the PET architecture, trained on the bulk water dataset.
PET-C
PET-NC
A modified PET architecture, trained on the bulk water dataset.
SOAP-BPNN-C
A SOAP-BPNN architecture, trained on the bulk water dataset. A modified COAP-BPNT architecture, trained on the bulk water dataset to directly predict forees.
SOAP-BPNN-NC
Table 4. Models used in the present work.
should be noted that:
The ORB model (orb-v2) is more accurate than MACE-MP-0 and SevenNet, as the former is pre-trained on the Alexandria (Schmidt et al., 2024) dataset and then fine-tuned on MPtrj, while the latter two are only trained on MPtrj. Despite its higher accuracy, ORB yields problematic physical behavior as discussed in this work.
The PET-NC and SOAP-BPNN-NC models are simply obtained from the respective conservative models by changing the output head to predict atomic forces directly. In the case of SOAP-BPNN-NO an equivariant vector representation is generated internally thanks to formalism in Villar et al. (2021). In both cases, due to the marginal increase in number of parameters in the force head, the con-conservative models have slightly more parameters than their conservative counterparts. Within the calculators for these non-conservative models, implemented the net force removal suggested in Neumann al. (2024).
The architectures of these models are further described here:
ORB: rotationally unconstrained and non-conservative architecture, presented in Neumann et al. (2024).
MACE: A rotationally invariant and conservative architecture, presented in Batatia et al. (2022).
SevenNet: The SevenNet model (Park et al., 2024) makes use of the NequlP (Batzner et al., 2022) architecture, which is rotationally invariant and conservative.
PET-C and PET-NC: A re-implementation of the architecture in Pozdnyakov & Ceriotti (2023), which is rotationally unconstrained and conservative. The non-conservative version changes the final head to predict forces instead of energies.
COAP-BPNNC and using SOAP (Bartók et al., 2013) descriptors. This architecture is rotationally invariant and conservative. The non- conservative version makes use of the formalism in Villar al. (2021) to predict forces (a vector) from scalar internal representation.
Behler-Pariinello neural network architecture (Behler & Parrinello, 2007),
D.I. Timings of general-purpose models
Table 5 shows the timings for the four eeneral-pupposs models tested this work, compared with the PET models we train here. The large version MACE-MP was used in this table and throughout this work. The small version of EquiformerV2 was used in this table and throughout this work, except to calculate work loops, where the large version was used.
16

---

The dark side of the forces
Table 5. Timings (ms) of general-purpose purpose models on a bulk water NVE simulation with 64 water molecules (192 atoms), Nvidia H100 GPU. Note that timings refer to overall evaluation time for a single configuration, performed using the ASE calculator distributed with each model, and therefore may contain overheads that are not directly associated with the model evaluation (neighhor list construction, etc.). The timing of the water-focused PET-C and PET-NC are also reported, for reference.
MODEL
TIMING PER STEP
TIMING PER STEP PER ATOM
MACE (C)
26.9
0.140
SEVENNET (C)
52.8
0.275
ORB (NC)
11.9
0.062
EQUIFORMER (NC)
1580
8.230
PET (C)
19.4
0.101
PET (NC)
8.58
0.047
ARCHITECTURE TYPE
TRAINING
MAB(V)
MAE(f)
TIMING (TR.)
TIMING (EV.)
PET
V
4.7
1025.6*
5.48
0.0264
PET
C
1.26
18.6
15.30
0.0713
PET
C
V.F
0.55
15.31
0.0716
PET
NC
24.3
5.55
0.0224
PET
NC
V,F
1.42
24.8
5.63
0.0269
C
20.2
0.0715
PET-M
V,F
0.59
15.43
NC
26.7
0.0265
PET-M-FT***
C
V,f
0.50
20.0
56.42**
0.0714 14
SOAP-BPNN
2.16
177.0*
3.57
0.1065
SOAP-BPNI
C
1.89
40.6
36.10
0.6394
SOAP-BPNN
C
V,
1.38
41.4
36.20
0.6367
SOAP-BPNI
NC
112.2
5.81
0.1515
SOAP- BPNN
NC
V,
3.20
111.9
6.34
0.1674
Table 6. Test errors (energies in meV per atom, forces in VeÀ) training and evaluation timings for models trained on the bulk water dataset. Training timings correspond the time to compute a single epoch (in seconds) on H100 GPUs with total batch size 64. Evaluation timings correspond to the average time per atom (in ms) for energy and/or force evaluations for single structures across the test set.
*The force emorr of energy-only models computed by evaluating forces derivatives of the energies, despite the fact that explicit
force training took place.
** A linear fit executed minimize training errors energies for the force-only models, in order calculate the best constant shift for the fictitious energy of which the forces are the derivatives. **eTrained single GPU opposed four, and for single day opposed two.
17

---

The dark side of the forces
D.2. Accuracy and timings of water models
Similarly to Table 1, we evaluate the accuracy of the SOAP-BPNN architecture on the same bulk water dataset under different training conditions. The results in Table 6 show, again, that the lack of energy conservation without corresponding data augmentation strategy seems to hurt the accuracy of the models.
In general, although the non-conservative models trained in this work on the bulk water dataset (with little more than 250000 targets) seem to show worse accuracy. training on larger datasets has shown that non-conservative models can be competitive in accuracy with conservative models. This is not only because energy conservation can then be effectively learned, but also because non-conservative models, by virtue of being faster, can train for larger number of epochs a given computational budget. Training duration be the limiting factor accuracy large datasets.
The timings of the PET models are fully consistent with the theoretical cost analysis in Appendix B. In contrast, the SOAP-BPNN models rely on the SOAP atomic descriptors cwm/met met http://etttteereeeeeeeeeee /ftnooocr/eeetoic...... Within this library for atomistic descriptors, three imple- mentation details account for the SOAP-BPNN timings: 1) although the models are trained and evaluated on GPU, the feature calculation executed on CPU: 2) feature calculation is parallelized across different structures, but not different chemical environments within the same structure (effectively meaning that no evaluation time parallelization is present); 3) the equivariant calculation of forces makes it necessary evaluate additional features with angular momentum quantum number of 1.
implemented in
E. conncoonsevvttive work
Since non-conservative models do not obey equation (5), we evaluate indicative magnitudes of the work over a closed loop for the non-conservative models considered this study. These are shown in Figure
The closed path corresponds the rotation and deformation of single water molecule within a liquid water structure, while keeping all the other atoms fixed. Figure shows the cumulative work of the models considered in this study. Even though the cumulative work curves are very similar, the conservative models results in zero overall work on the closed path; meanwhile, the non-conservative models leads to an overall non-zero work along the path.
Figure Cumulative work along a closed path for all models considered the study. While the total work for the conservative models zero up machine precision, the non-conservative models exhibit 132 VeV, SOAP-BPNN-NC: -410 meV). The first figure from the left shows the overall path, while the other two zoom on the last part.
total work (ORB: meV, Equiformer: -241 meV, PET-NC:
0.10
0.4
PET-NC
PET-C
Cumulative work ev
SOAP-NC
ORB
SevenNet
-10-
Equiformer
MACE-MPO
0.05
0.2
-20-
0.0
0.00
-30
0.2
-0.05
40
-0.4
-0.10
500
1000
980
990
1000
990
995
1000
Structure number
Structure number
Structure number
18

---

The dark side of the forces
F. Geometry optimization
F.I Optimization trajectories
We compare the behavior of different models when quenching a liquid- snapshot, optimizing it towards the nearest potential energy minimum, (Figure 7) showing both the convergence in terms of the force modulus, and the trajectory of configurations latent-space projection built on geometric descriptors. The latent-space plots obtained by computing SOAP descriptors for all configurations in all the trajectories, averaged atomic centers in each structures, and projected on the axes of highest variability using simple Principal Component Analysis (PCA). The qualitative features of the latent space plots are insensitive to the details of the SOAP descriptors used.
10
10
SOAP-BPNN-NC SOAP-BPNN
PET-NC
SOAP-BPNN-NC SOAP-BPNN
PET-NC
PET
PET
10
ORB
100
ORB
/vv/A
10
PCA[1]
[f| /eee
10
PCA[1]
10
PCA[2]
10-2
PCA(2
10-3
10-3
10-4
10
10a
1
10
10a
102
10a
step
step
Figure 7. Force modulus function of step along FIRE (left) and BFGS (right) optimization trajectories of liquid water snapshot, for different MLIPs. The inset shows a latent-space projection of the trajectory in configuration space; the triangle indicates the starting configuration, and the squares the final one
We discuss first the results for the FIRE algorithm (left panel Figure 7). Despite the very different level of accuracy, the conservative SOAP-BPNN model and PET converge to similar configuration (with PET forces saturating at about 10-3 eV/Á, due to numerical precision). The non-conservative PET-NC model also converges to similar structure, indicating that sufficiently accurate non-conservative model can indeed be used with gradient-based structural optimization algorithm. The non-con non-conservative SOAP-BPNN model, however, displays a catastrophic mode of failure, with the force never decreasing below 0.1 eV/Á, and the trajectory drifting off in a different direction without ever reaching stable state. ORB converges to a very different structure than the other models, which might be due to the different reference DFT functional used for training.
The BFGS optimization trajectories (right panel in Figure 7) very similar to those obtained with FIRE, except for SOAP-BPNN-NC whose catastrophic failure is apparen in the trajectory going in completely different direction, and in the force modulus never converging below 0.1 eV/À. Thanks to the second-order nature of LBFGSS convergence is much faster, until PET-based models reach their precision limit. ORB shows an interesting behavior, as it reaches similar configuration as with FIRE, but the force modulus exhibits large fluctuations, and would be hard to determine a clear threshold to establish convergence. On a practical level, this experiment shows that non-conservative forces make geometry optimization more fragile- with first-order methods being somewhat more stable although slower and require careful choice of the minimization algorithm and its convergence parameters.
F.2. Failure rates
In Table 7, geometry optimization is attempted with range of conservative and non-conservative models. For each model, three cases are considered: geometry optimization of gas-phase water molecules, starting from the experimental geometry, randomly displacing the coordinates with a standard deviation of 0.5 Å, and relaxing the geometry; 2) geometry optimization of bulk water structures from the test set of (Cheng et al., 2019), 3) geometry optimization of molecules chosen at random from the QM9 dataset (Ramakrishnan et al., 2014), randomly displacing the coordinates with standard deviation of 0.5 À before relaxing the structures. Geometry optimization is performed with the L-BFGS (Liu & Nocedal, 1989) algorithm as implemented ASE (Hjorth Larsen et al., 2017). Optimization runs that do not converge within 1000 optimization steps are considered as failed.
19

---

The dark side of the forces
Table 7. Percentage success in geometry optimization.
MODEL
(1)
(1)
QM9
MPTRJ
PRB-LOW-PRECIIION (NC)
3
0
10
ORB (NC)
69
1
76
SEVENNET (C)
92
MACE (C)
94
83
94
99
PET-NC (NC)
75
PET-C (C)
83
58
BAAP-NPNN- (NC)
79
SOAP-BPNN-C (C)
91
59
*This an ORB model used with default settings, which lower precision of matrix multiplications. Given results here, we deactivated this setting for all ORB results shown in the rest of this work.
Non-conservative models consistently show lower rates of success in geometry optimization. It should be noted that striet convergence criteria were used (Emax=1e-5A for molecules and fmax=1 1--4A for bulk systems).
20

---

The dark side of the forces
G. Additional molecular dynamics results
G.I. NVE temperature profiles
We repeat the experiments Sec. 4.4 using all the potentials we tested. As shown Figure 8, all conservative models yield stable kinetic temperature profile (there spread in the average temperature that compatible with the nature of constant-energy trajectories) and all the non-conservative models show noticeable temperature drift. Equiformer displays drift comparable to PET-NC, which is remarkable for ceneral-purpose model (but is too expensive be used in practical simulations). SOAP-BPNN-NC shows extremely large drift, consistent with the fact it is very inaccurate model.
800
PET-NC SOAP-NC Equiformer
PET-C SOAP-C SevenNet MACE-MPO
700
600
500
400
300
10000
20000
30000
40000
step
Figure thorough comparison kinetic temperature profiles for NVE trajectories performed with different conservative (right column in the legend) and non-conservative (left column in the legend) ML potentials.
Graphene
Amorphous carbon
FCC aluminium
Benzene on graphene
Alanine dipeptide in water
600
400
ORB
200
SevenNet
10000
20000
10000 20000
10000 step
20000
10000 step
20000
10000 step
20000
stop
step
Figure 9. Kinetic temperature profiles for NVE MD trajectories, using three universal models, for three different materials, as indicated. Simulations for graphene and amorphous carbon were performed with time step of fs, and those for aluminum with time step of reflecting the different timescale of atomic motion.
We also perform NVE MD trajectories on a wider range of chemical systems using the general-purpon potentials MACE-MP, SevenNet, Equiformer and ORB. Fig. 9 shows that the non-conservative models ORB and Equiformer exhibit strong kinetic energy gain during the simulation of graphene, amorphous carbon and FCC aluminium, with ORB exhibiting an explosive behavior for graphene (possibly because it is outside the training set). The conservative models exhibit stable temperature fluctuations, even for amorphous carbon where the high distortion of the structure leads to higher initial energy and broader kinetic energy fluctuations. The order of magnitude of the drift is similar for systems with different degrees of chemical heterogeneity.
21

---

The dark side of the forces
G.2. NVT
We repeated the bulk water experiments in Sec. 4.5 with wider range of models and thermostats. The results summarized in Table 8.
Table Average kinetic temperatures (in K) for different models and thermostatting schemes. For SVR, we show also the separate kinetic lemperature computed only for bydrogen and oxygen atoms. All thermostat time constants are reported in units of fs.
THERMOSTAI TYPE
WN
WN
WN
SVR
SVR 10
SVR
THERMOSTAT
10
100
1000
10
ATOMS
ALL
ALL
ALL
H
O
ORB (NC)
300.4...1)
(04.2(0.2)
351 351.0/o.6
301.0.o.1
336.2.0.8)
229.5
MACE (C)
300. (a.1)
300.4(0.2)
300. (0..4
300.1
299. ()
302.3(.1)
SEVENNET (C)
300.0(o.1)
300.1
299.9(0.5)
300.1...
299.9(0.4)
300.5
PET-NC (NC)
300.1
301.4
312. (1.31
300.3
295.6,u
309 (0.69
PET-C (C)
300.1
300. (o.2)
300.3
300.1.o.1
299.6(0.3)
301.0(o.7)
SAAP-BPN--C (NC)
301.9
340.7 (
1.2 10
309. (0..7.o.u
265.5(0.2)
399.5
SOAP- CAAP-BPNN-C(O (C)
299.9(0.1)
299.8(0.4)
300.5g0.5)
(0..1
298 (
302.4(1.5
G.3. Direct stresses: NPH and NPT
In order test the effects of direct stress prediction, which can be employed in constant-pressan simulations, we run MD in the NPH and NPT ensembles. In this case, means disentangle the effects of non-conservative stresses from those of non-conservative forces, run the "non-conservative" (NC) models using non-conservative stresses and conservative forces. For these simulations, employ the PET-MAD general-purpon potential (Mazitov et al., 2025), whose version 1.1 is trained to provide both conservative and non-conservative stresses.
15
PET-MAD-C, NPH PET-MAD-NC. NPH PET-MAD-C, NPT PET-MAD-NC, NPT PET-MAD-M, NPT MTS-8
Cell edge A
14
13
12
2
4
8
10
Time / ps
Figure 10. Length of the cubic cell edge, as a function of time, during NPH and NPT simulations of liquid water using the PET-MAD model, in conservative mode, non-conservative mode, and using multiple time stepping with M 8. The Lesss-Zykova-Parrinng harostat (Bussi et al., 2009) was used in all simulations.
Fig. 10 shows how the conservative runs produce reasonable cell volumes, both the constant-enthalpy NPH ensemble and the cathermal-isobarri NPT ensemble. The concconservvttve NPH run shows rapid and unphysical expansion of the cell, which reflected in the rapid growth in the enthalpy of the simulation (which instead conserved very well by the conservative NPH run). This "explosion" of the cell is prevented by the thermostat in the NPT ensemble, although the exact value of the cell parameter is inconsistent with the conservative NPT baseline. Finally, we show MTS (multiple time stepping, see App. I) run using M 8, which correctly reproduces the cell volume of the conservative NPT simulation while only evaluating the more expensive conservative stresses every time steps. This simple example illustrates how all the considerations that were made for non conserrvvve forces in the context of NVE/NVT simulations are also applicable to NPH/NPT simulations when using direct stress predictions.
22

---

The dark side of the forces
600
500
Temperature /
400
300
200
PET-MAD-C
PET-MAD-NG
PET-MAD- C.
O
PET-MAD-NN NC. PET-MAD-NE NCC
PET-MAD C, H
100
0
2
4
6
8
10
Time / ps
Figure I1. Total and atom-lype-resolved kinetic energy profiles for NVE simulations of liquid water using the conservative and non- conservative modes the PET-MAD (Mazitov et al., 2025) potential. non-conservative case, energy conservation is enforced by velocity rescaling.
G.4. Energy conservation by velocity rescaling
The easiest way to stabilize NVE molecular dynamics using direct forces is rescale the velocities (or, equivalently, the momenta) at each step to enforce conservation of energy, proposed in Bigi et al. (2025). While this leads to an improvement, in the sense that NVE trajectories are guaranteed to be stable, many fundamental issues remain. As an example, Fig. 11 compares traditional NVE MD with a conservative model to NVE MD with a non-conservative model, where the latter is stabilized by velocity rescaling. The conservative model shows kinetic temperatures close to the temperature at which the starting structure was equilibrated, both for the overall simulation and for individual atomic types. In contrast, the non-conservative run departs from the equilibration temperature and, more importantly, shows much different average kinetic energies for atoms belonging different chemical species. This behavior is in violation of the principle of equipartition of energy, and it manifestation of the fact that the target ensemble (NVE, in this case) being sampled correctly. This is very similar to what we observe with global thermostats, and is analogous manifestation of the fact that enforcing the correct values of global quantities does not guarantee that all the parts of a simulation will equilibrate correctly, in the absence of a well-defined Hamiltonian.
23

---

The dark side of the forces
G.S. Structural correlations
Even though the kinetic temperature can be used diagnostic of the ensemble error, does not guarantee that the trajectory yields structural properties consistent with the correct NVT ensemble. We compute different kinds of structural conrelation functions: the O-O and H-H pair correlation functions goo(r) and (Hи(r) reports on the probability of finding two oxygen or two hydrogen atoms at a given distance and the orientational correlation function Cuu (r) reports the average scalar product between the orientation vector of two water molecules at a given distance. Results for the PET models (Figure 12) show that structural correlations are consistent for all thermostatting schemes for the conservative runs. The 10 white-noise Langevin runs exhibit larger error bars, consistent with the poor sampling efficiency, All non-conservative PET-NC runs show statistically significant discrepancies with respect to the conservative reference. Although it is not possible rigorously disentangle sampling errors due non-conservative forces, and the overall difference in the force values between PET and PET-NC, we hypothesize that the differences seen in the LooongLangevin--e Langevin-coupling limit are due to differences in the underlying forcefield, while the differences for the weaker-coupling runs are associated with the deviation in the sampling temperature.
and global SVR thermostat correlations with correct sampling. The large difference in reference correlations with respect the custom-trained PET models is due the difference in reference DFT energeties: the PBEsol functional used MPTraj is known to yield severely overstructured water, while the revPBEO-D3 functional used in (Cheng et al., 2019) was shown give good agreement with experimental data. For the purpose of this work, however, this difference is not important, and the main takeaway from Figure 13 is that for ORB there are large discrepancies between different thermostatting schemes, which is indicative of the artifacts induced by non-conservative forces, and the difficulty in mitigating them using thermostats.
For the ORB models (Figure 13) we consider the high-coupling limit as the most reliable reference of the expected
24

---

The dark side of the forces
Reference (PET-C, SVR)
A-ref. PET-C
A-ref. PET-NC
0.10
SVR == 10 fs
WN T= 1 ps
WN T= 100 fs
0.05
WN 10
goo(r)
0.00
00.05
2
5
2
0.05
(mn(()
0.00
-0.05
1
2
4
5
2
3
4
5
2
3
5
0.04
0.4
0.02
Cou(f)
0.2
0.00
0.02
0.0
-0.04
3
//A A
r/ A
r/
Figure 12. Structural correlations for simulations of room-temperature water using different PET models. From top to hottom, the rows report the O-O correlation function, H-H correlation function, the orientational correlation function. The left column reference value (SVR thermostat for the conservative PET model), the middle column the differences with respect to this reference for PET-C using different thermostats, and the right column differences to the reference using the non-conservative PET model and different thermostats, Shaded encompass two standard emrrr around the values.
25

---

The dark side of the forces
(ORB,
fs)
Á-ref.
0.75
SVR
3
0.50
WN T= 100 fs
WN T ps
goolr)
2
0.25
1
0,00
-0.25
0.1
9malr)
-0.1
2
0.1
0.4
0.0
Cuul))
0.2
-02
0.0
-0.3
3
5
3
5
r/
r/A
Figure 13. Structural correlations for simulations of room room-temperanee temperatrrter water using the ORB model and different thermostatting scheme. From top to bottom, the rows report the O-O correlation function, the H-H correlation function, the orientational correlation function. The left column shows the reference value sstrogg-coullnng differences with respect to this reference for ORB trajectories using different thermostats. Shaded areas encompass two standard emoor around the mean values. Note the much larger range of the discrepancies with respect the PET-NC runs in Figure 12.
10 is white-noise Langevin thermostatting) and the right column the
26

---

The dark side of the forces
G.6. Generalized Langevin equations
The generalized Langevin equation (GLE) thermostat (Ceriotti et al., 2010) local stochastic thermostat that extends traditional Langevin dynamics by introducing history-dependen friction term. The possibility of tuning the memory kemel of the friction and the noise allows tuming the effect of the thermostat, by effectively coupling different modes of the physical system to different friction timescales. Thanks to this construction, near-optimal damping can be achieved for large range of characteristic frequencies, as opposed to the traditional Langevin thermostat, which only achieves optimal damping (and therefore sampling) for narrow frequency range.
Among the many different applications of GLE thermostats, that of Ref. 50 appears particularly suited to contrast the temperature drift observed with the direct prediction of forces: by applying high-pass filter the noise, allows to thermostat aggressively the fast degrees of freedom of the system, without disrupting the long-time dynamics and the sampling of the diffusive degrees of freedom. We use the same parameters that have been used successfully to control temperature drift caused by the integration errors for aqueous systems in Ref. 50, using 1/70- 83.33 ps, 1/76o = 10fs, and wy = 300 ps where all GLE parameters refer to their definition in Morrone et al. (2011).
We probe the behavior of conservative and non-conservative models under #LE-thermostatted dynamics. Table shows that this GLE, that could successfully control integration errors for an empirical water model, fails catastropically when used with non-conservative models. This can be explained by the observation that the near-diffusive modes, which are those that show the most non-conservative behavior according the analysis in Figure 1, are almost not damped all by this thermostat, Even though one could in principle design low-pass filter, or use more aggressive parameters for this high-pass thermostat, the analysis performed in Ref. 50 indicates that such GLE would slow down long-time sampling, exactly as it is the case for a traditional Langevin thermostat. By reducing the sampling efficiency, this effect negates the advantages of direct force prediction.
Table Average kinetic temperatures for different models and atomic species, for the molecular dynamics of water using GLE thermostat.
THERMOSTAT TYPE
GLE
GLE
GLE
ATOMS
II
ORB (NC)
(181.5(5..)
1099.2(4.))
1362.5
MACE (C)
303.1
301.7a.0
306.0(166)
SEVENNET (C)
304.1
302.6(1.))
307.3(2.4)
PET-NC (NC)
524. (4.6)
495.5
585.5
PET-C (C)
301.8 (012)
301.6(1.2)
302.2(1.6)
SOAP-BPNN-NC (NC)
10
2.1 10
5.9 10?
SOAP-BPNN-C (C)
301.3(1.6)
301.0(.4)
301.9(.8)
27

---

The dark side of the forces
H. Conservative fine-tuning
Figure 14 and Figure 15 show that a non-conservative model can be trained to its asymptotic accuracy with 2-4 times smaller investment of computational resources. This is consistent with previous observations (Gasteiger et al., 2021), and one of the main reasons for which direct prediction of the forces has been proposed in previous work. More importantly the figures also show training curves corresponding non-conservatively), and whose V head was then "fine tuned" using energy and back-propagated forces (while the non- conservative force head was also trained to avoid its degradation). This "PET-M-FT" model achieves the accuracy of a conservative model trained "from scratch" in about 1/3rd of the total training time, demonstrating that this simple fine-tuning strategy makes it possible to recover most of the speed-up afforded by direct-force force training, while having the convenience and accuracy of conservative model. might also be possible (although technically more complicated) alternate gradient steps using only the direct head and steps using the conservative forces achieve a similar effect.
model that has been initially trained with direct V, f heads (i.e.,
Figure 14. Training curves (showing validation MAE force as function of the GPU time expenditure NVIDIA H100) for conservative PET-C model, for the non-conservative head a PET-NC model, and for the fine-tuning of a conservative model initialized from the potential energy head of the PET-NC model. Training conducted the liquid water dataset.
400
PET-C PET-C, fine-tuned
Val. force MAE/ / /mev/A
PET-NC
200
100
50
20
10
20
50
100
200
GPU hours
Figure 15. Training (showing validation MAE force error function of the GPU time expenditure NVIDIA H100) for conservative PET-C model, for the non-conservative head of a PET-NC model, and for the fine-tuning of conservative model initialized from the potential energy head of the PET-NC model. Training was conducted on small subset of OC20 (Chanussot al., 2021) (training on the first 15000 samples of the 200k S2EF training set, and using the next 5000 for validation).
200
Val. force MAE mev/A
PET-C
PET-NC
PET-M, fine-tuned
100
70
2
5
10
20
50
100
GPU hours
It should be noted that models trained with this fine-tuning procedure will offer accurate conservative and non-conservative forces, making them suitable for the multiple time stepping strategy discussed in Section 4.8. The overall procedure would allow the resulting models to be trained and evaluated at the computational cost non-con conseevserveeive models, while achieving the levels of accuracy and energy conservation of conservative models.
28

---

The dark side of the forces
Type
Relative speed (0 on Sc-Co-Al alloy)
Relative speed (water)
Relative speed (theoretical)
NC
1.0
1.0
C
1.93
1.92
F 1.925
MTS, M 4
1.50
1.45
1.48
MTS, M/ 8
1.19
1.18
1.24
MTS, M = 16
1.05
1.04
1.12
Table 10. Empirically observed slowdown conservative MD simulations multiple time-stepping with different strides M, well simulation where the conservative forces are evaluated every step, compared with the theoretical cost + F/M with F F= 1.925. Timings are given for selected systems from the experiments Section 4.4 and Appendix G, using runs NVE ensemble.
I. Multiple time stepping
Multiple time stepping (MTS) is based on a Trotter splitting of the Liouville operator (Tuckerman et al., 1992), and can be applied any time one can decompose potential into a "slowly varying" and expensive compute part Vstow and a "quickly varying", inexpensive part Visa- The dynamics is propagated using Via and small discrete time step A.. Every M steps, one applies correction based on Valow. In the case of a Velocity- Verlet integrator (Verlet, 1967), a MTS integrator amounts to:
p - p + M falow At/2
p - p + fins At/2
q a 1ttpAtnn m
Repeat M times
(12)
p p fra At/2
p f p Mf At/2
where q and p indicate the vectors containing positions and momenta of all particles, and m the masses. Forces have to be recomputed if they are needed after position update. In the limit of ideal time scale separation, the resulting dynamics (and the sampled ensemble for constant-temperatuue simulations) is consistent with Vslow at much-reduced cost: Letting F denote the slowdown of "slow" over "fast" forces, and assuming that an evaluation of "slow" forces also yields the "fast" ones, MTS reduces the theoretical cost of well-behaved molecular dynamics from therefore requires a trade-off between accuracy and speed. Morrone et al. (2011) discusses the failure modes of MTS high M. and shows how to use thermostats for mitigation.
(F 1)/M. Choosing M
Here we use the implementation in i-PI (Kapil et al., 2016; Litman et al., 2024), and use the non-conservative forces as the "fast" forces, fist fNc. and the difference between conservative and non-conservative forces as a slow correction, flow = fc - fNC. Our current implementation does not reuse computation between the conservative and non-conservative forces, and so the practical overhead given by 1 + F/M, with an empirical factor of 2 (see Appendix B). As seen in Table 10, empirical timings match the theoretical ratios rather closely.
NVE MD trajectories based the MTS algorithms stable up to time step factor of M 8, while trajectories with M = 16 become unstable and show large kinetic temperature drift (Figure 16). Note that the cause of drift here is different than for non-con non-conservative models, and is associated with resonances between the natural dynamics of the fast degrees of freedom and the integration errors due to imperfect time scale separation. Running at M = 8 recovers almost in full the speed-up associated with the use of a non-conservative force head, while allowing to run stable constant-energy trajectories. Using a SVR thermostat makes possible to sample the NVT ensemble, with excellent agreement with the reference PET-C trajectories (Figure 17). The unstable MTS-16 trajectories, on the other hand, show large structural anomalies, that might be in part corrected by application of suitably tuned colored-noise thermostat (Morrone al., 2011). Given that the computational savings would be negligible, however, there is little reason to M greater than when using non-conservative forces in the fast step.
29

---

The dark side of the forces
800
PET-C, NVE PET-NC, NVE PET-M, NVE MTS-4
700
PETM. NVE MTS-8 PET-M, NVE MTS-16
600
500
400
300
10000
20000
30000
40000
step
Figure 16. A comparison of kinetic temperature profiles for NVE MD trajectories of bulk water performed with the PET-M model and different orders MTS integration. The step numbers refer to the "fast" evaluation.
Reference (PET-C, SVR fs)
A-ref. PET-M, MTS
MTS, M =
2
0.4
MTS, M 16
MTS, M 8
goo(r)
0.2
0.0
-0.2
3
gre(())
-0.1
2
2
0.10
0.4
0.2
0.00
0.05
0.0
-0.10
3
5
3
5
r/
r/A
Figure 17. Structural correlations for simulations of room-temmpertuure water using PET-M model, using multiple time-step integrators alternte the evaluation of non-conservative forces and (every M steps) the conservative model. From top bottom, the report the O-O correlation function, the H-H correlation function, and the orientational correlation function. All simulations use an efficient and gentle global SVR thermostat, with
fs coupling time. The left column shows the reference value (full PET-C MD) and the right column the differences with respect to this reference for MTS trajectories using different M values. Shaded areas encompass two standard errors around the mean values.
30

---

The dark side of the forces
We also investigate the behavior of MTS for the more generic PET potential trained on the OC20 dataset in Appendix H. We the model after conservative fine-tuning, which is able predict both conservative and non-conservative forces with good accuracy. The simulations were run system of single oxygen atom adsorbed on an alloy surface, and their kinetic energy profiles are shown in Figure 18. Also in this case, well-behaved molecular dynamics can be performed up to a high multiple time-step factor M 8, a value that retains nearly all the speed-up afforded by non-conservative forces.
400
PET-C, NVE PET-NC, NVE PET-M, NVE MTS-4 PET-M, NVE MTS-8 PET-M, NVE MTS-16
350
300
250
200
150
0
20000
40000
60000 step
80000 100000
Figure 18. A comparison of kinetic temperature profiles for NVE MD trajectories of oxygen atom adsorbed alloy surface, performed with the various heads of a model which was produced by "conservative fine-tuning" on subset of the OC20 dataset. The step numbers refer the "fast" evaluation.
J. Software
Molecular dynamics simulations have been performed using the i-PI software (Litman et al., 2024), and geometry optimiza- tions using ASE (Hjorth Larsen et al., 2017).
ORB, MACE-MP-0, Equiformer, SevenNet calculations were performed using the publicly available repositories, www.mmmaleeeereeaeaaaaaaa models/ (orb-moded) models 0.4.0), https: (mace-torch 0.3. 8), https://github. com Cow./AIR-CCC///C/c//cecaa ffircar (sevenn 0.10.1).
respectively https: / /gtttub..tttttgg attp://eeeeotiiaaaaaaataaaaaaaaa //github. wwt.///////tttt/meeemee (fairche--core 1.3.0) 0), and https::pp/ pettps//gitthue Sttp:a/.././/O...SS/SSSSSS
31


Reviews:
[{"review_id": "w1p9LrgXND", "reviewer": "Reviewer_K2Jj", "summary": "This paper critically evaluates the implications of using non-conservative (NC) force models for learning machine learning interatomic potentials (MLIPs). This paper investigates the general accuracy of the NC models' forces, the effect of using NC force models on measured values of specific properties of interest in atomistic simulations, and geometric optimization. \n\nWhile direct force prediction using NC force models provides a faster method for running atomistic simulations, the paper shows that even the best NC force models result in unstable simulations. The authors show that for a bulk water simulation, the NC forcefields result in a drift equivalent heating rate of about 7000 billion degrees per second. The authors also argue that some well-known correction procedures do not alleviate these issues. Using local thermostats to control the instability/run-away behaviors of NC force model-based simulation results can lead to slower convergence as it leads to a smaller diffusion coefficient, making it impossible to measure any time-dependent properties. Similarly, global thermostats can lead to the introduction of artificial artifacts as the non-conservative terms have different behaviors on different degrees of freedom. Through this experiment and other geometric optimization experiments, the authors convincingly show the unstable behavior introduced by using NC force models for atomistic simulations and geometric optimization. \n\nFinally, the authors also provide a multi-step integration procedure for using hybrid models to effectively reap the benefits of the training and inference speed of NC force models while also enjoying the stability of conservation force models.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the Reviewer for their comments and questions. To provide more context, the use of multiple time-stepping to reduce the cost of a simulation while matching the accuracy of a “slow” potential is well-established in the molecular dynamics community. It was originally introduced to reduce the cost of long-range electrostatics (which is more computationally demanding, and harder to parallelize, than short-range components in a physically-motivated model) [Tuckerman 1992], and it has been used extensively together with different levels of electronic-structure theory [Marsalek 2016] and with machine-learning potentials [Rossi 2020]. It has also been recently applied to a very similar problem - namely correcting for symmetry breaking in unconstrained machine-learning models [Langer 2024]. In all these cases, and in that of non-conservative forces, MTS acts as a mitigation strategy, with smaller initial errors allowing for a larger stride M. \n\n\nNormalizing the cost of non-conservative MD to 1, and letting $F$ denote the slowdown for conservative force predictions, $2 \\leq F \\leq 3$, MTS reduces the theoretical cost of well-behaved molecular dynamics from $F$ to $1+(F-1)/M$. The evaluation of the conservative force is reduced by a factor $1/M$ and there is one common evaluation every $M$ steps. Note that in our current practical implementation, we do not yet re-use data for the common evaluation, and therefore incur a cost of $1+F/M$. We find that for the PET model employed in this work, $F\\approx2$.\n\nIn all our MTS experiments, the empirical timings match the theoretical ratios rather closely:\n\n| Type        | Relative speed (O on Sc-Co-Al alloy) | Relative speed (water) |\n| ----------- | --------------------- | -------------------- |\n| NC         |          1.0           |  1.0                    |\n| C            |          1.93        |  1.92                    |\n| MTS, M=4  |          1.50                  | 1.45                     |\n| MTS, M=8  |          1.19                  | 1.18                    |\n| MTS, M=16 |          1.05                  | 1.04                     |\n\nFor $M=8$ and $F\\approx2$, the simulation is only $\\approx20$% slower than using exclusively non-conservative forces. Note that, even in the most unfavorable case (where only $M=2$ is stable), the theoretical computational overhead due to the use of conservative forces could be halved with an efficient implementation of MTS.\n\nEven though these general arguments suggest that MTS should be broadly applicable, we do agree with the Reviewer that further evidence that “hybrid” models are useful in general would be very valuable. To this end, we also consider the case of the diverse OC20 dataset (containing molecular fragments adsorbed on different catalytic surfaces) [Chanussot 2020].\n\nFirst, we performed a “conservative fine-tuning” exercise - reducing the training time by a factor of more than 3 while achieving the same final accuracy as the conservative model. This shows that the idea of pre-training with direct forces and then turning on backpropagation to fine-tune, therefore obtaining a hybrid model inexpensively, is far from specific to the water system. (Incidentally, this training method that we propose has already been picked up by the community to train state-of-the-art models, for example in [Fu 2025].)\n\nSubsequently, we use this model to perform MD on a representative structure (an oxygen atom adsorbed on an alloy surface). Even though the direct-force model shows a much worse temperature drift than observed in the experiments on water (consistent with the lower accuracy of the fit for this more broadly applicable potential) MTS still succeeds in stabilizing NVE trajectories up to M=8 (https://imgur.com/a/PY96n58). This provides a concrete example that MTS also works for less accurate, general-purpose direct-force models.\n\nTo summarize, our work not only highlights the potential problems related to using non-conservative forces, but it also provides simple and effective solutions in the form of conservative fine-tuning (at training time) and MTS (at inference time), which allow practitioners to retain nearly all the advantages of using direct forces while avoiding their physical pitfalls.\n\n[Tuckerman 1992] https://doi.org/10.1063/1.463137\n\n[Marsalek 2016] https://doi.org/10.1063/1.4941093\n\n[Rossi 2020] https://doi.org/10.1021/acs.jctc.0c00362\n\n[Langer 2024] https://doi.org/10.1088/2632-2153/ad86a0\n\n[Chanussot 2020] https://doi.org/10.1021/acscatal.0c04525\n\n[Fu 2025] https://arxiv.org/abs/2502.12147", "claims_and_evidence": "The authors claim that the NC force models do not implicitly learn force conservation and that the direct use of Non-conservative forces for atomistic simulation can lead to unstable simulations. The experimental setup and the evidence provided by the authors convincingly prove these claims.\n\nThe authors also claim that multi-step integrators with a hybrid multi-headed conservative and non-conservative force model can produce fast and stable simulations. While the authors show this claim to hold in one experiment, it is unclear if this still holds for more general systems and if the optimal M (shown to be 8 in the experiment considered) still leads to a non-trivial speed up in the simulations.", "methods_and_evaluation_criteria": "The proposed metric of using lambda to understand the local deviations from the conservative fields and the measurement of the average kinetic energy and the velocity-velocity correlation function serve as appropriate and relevant measures for understanding the stability of the considered force models. Using the bulk water model provides a large and simple enough system to test the various problems about the stability of the NC force models.", "theoretical_claims": "N/A", "experimental_designs_or_analyses": "The experimental design and analysis is sound.", "relation_to_broader_scientific_literature": "As stated in the paper, many new works in MILPs are starting to look into direct force predictions without enforcing energy conservation. Generally, it has been seen that using physics-agnostic models trained at scale over enough data can lead to better models as the physical constraints can be learned from the data. This paper convincingly shows that ignoring energy conservation in learning Force models can lead to unstable atomistic simulations, and these problems cannot be easily fixed at inference time. This urges the field to either train better hybrid models using multi-head forces or develop better inference time methods to make the simulations more stable.", "essential_references_not_discussed": "The authors have provided a detailed overview of the relevant related works.", "other_strengths_and_weaknesses": "Strengths\n- The paper brings up an essential issue in the field of MILPs.\n- The proposed metrics and experiments provide a holistic picture of the problems associated with using Non-Conservative force models for simulations. \n- The paper also highlights some problems with using current methods for making inference time corrections to improve the stability of atomistic simulations run using NC force models.\n\nWeaknesses\n- The paper does not provide ample evidence to support using hybrid models to fix some of the highlighted problems.", "other_comments_or_suggestions": "N/A", "questions_for_authors": "Please refer to the previous sections.", "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": "", "extra_scores": {}}, {"review_id": "pai3mpMIeP", "reviewer": "Reviewer_NcLF", "summary": "The paper presents a systemmatic study of non-conservative force models for atomistic machine learning. Traditionally, forces are computed as the derivatives of potential energies to enforce symmetries and conservation laws. Several recent studies, however, directly predicted forces and learned the conservation laws implicitly during training to achieve a better trade-off between accuracy and computational efficiency. In this paper, the authors conducted microscopic simulations to investigate the non-conservative force fields. Experimental results show that these models suffer from several fundamental issues, from ill-defined convergence of geometry optimization to instability. In particular, energy conservation is harder to learn compared to rotational symmetry. The authors argue that the optimal way to use non-conserative force models is to integrate it with a conservative model, reducing the computation cost instead of completely eliminating the additional cost of backpropagation.", "questions": "", "limitations": "", "rating": 5, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for pointing out additional relevant literature, which we have added to the manuscript. We also added the recent related work by Eissler et al. (arXiv:2503.01431), in which a non-conservative and non-equivariant transformer model is trained on a large dataset and used to study the degree of energy non-conservation in MD simulations.\n\nWe chose water as a benchmark because it is an important system, ubiquitous in chemistry and biophysics, and because it is challenging due to the combination of high-frequency vibrations of the covalent O-H bond and the long-term dynamics of the hydrogen-bond network. Its relative compositional simplicity and homogeneity also means we can train a very accurate direct-force model, which allows us to show that non-conservative forces can produce unacceptable sampling artifacts even in this most favorable scenario.\n\nHowever, we agree that additional testing on more diverse molecular and materials systems is desirable. We had already included three examples (graphene, amorphous carbon, aluminum) in the appendices. Following the recommendation of this Reviewer (as well as Reviewer J2Ug), we performed additional tests on two heterogeneous systems: a benzene molecule adsorbed on a graphene surface, and alanine dipeptide solvated in water. We find that the non-conservative universal ORB and Equiformer models display severe lack of energy conservation in these cases, comparable to or worse than in the other systems we considered (https://imgur.com/a/CV3YmI3).\n\nRegarding the use of additional thermostats, we deliberately did not test Nosé-Hoover thermostats, which are non-ergodic in their original form and which, when using the “chain” extension [Martyna 1992] which improves (but does not fully fix [Patra 2014]) ergodicity by producing chaotic thermostat dynamics, amount essentially to a complicated strategy to generate poor-quality random numbers. We cannot think of a single case (besides perhaps the theoretical study of classical chaotic dynamics) in which N-H thermostats cannot be replaced by a stochastic thermostat. We agree with the Reviewer that it would be ideal to explore an “adaptive” thermostat to better substantiate our claim that thermostatting is not a valid fix for non-conservative dynamics. To that effect, we ran some tests using the “GLE-RESPA” [Morrone 2011] thermostats that use a bespoke generalized Langevin equation, tuned to stabilize noisy dynamics without disrupting the sampling of diffusive dynamics. Being local and adaptive, one could expect this approach to be preferable to a global (SVR) or plain Langevin thermostat for stabilizing direct force predictions. Our tests show that instead this GLE thermostat (which is less aggressive for slow degrees of freedom and hence does not slow down diffusion) is less effective at controlling the temperature drift induced by direct-prediction forces. The table below shows the atom-type resolved temperature for a simulation of bulk water at a target temperature of 300K, similar to table 2 in the manuscript. Results closer to 300K are better.\n\n| Model      | GLE (all atoms)         | GLE (H)           | GLE (O)           |\n|-----------------------|-------------------|-------------------|-------------------|\n| ORB (NC)              | 1181.5 ± 5.2      | 1099.2 ± 4.6      | 1362.5 ± 7.7      |\n| MACE (C)              | 303.1 ± 1.2       | 301.7 ± 1.1       | 306.0 ± 1.6       |\n| SevenNet (C)          | 304.1 ± 0.8       | 302.6 ± 1.3       | 307.3 ± 2.4       |\n| PET-NC (NC)           | 524.1 ± 4.6       | 495.5 ± 4.0       | 585.5 ± 5.9       |\n| PET-C (C)             | 301.8 ± 1.2       | 301.6 ± 1.2       | 302.2 ± 1.6       |\n| SOAP-BPNN-NC (NC)     | 1.6 × 10⁸         | 2.1 × 10⁸         | 5.9 × 10⁷         |\n| SOAP-BPNN-C (C)       | 301.3 ± 1.6       | 301.0 ± 1.4       | 301.9 ± 1.8       |\n\nThis is consistent with our observation that non-conservative effects are more pronounced for long-range correlations and slow dynamics and reinforces our claim that thermostats can mitigate, but not solve, the artifacts induced by non-conservative forces. Finally, we would argue that the MTS strategy we propose provides a much better and more broadly applicable fix to the shortcomings of non-conservative forces compared to thermostats.\n\n[Martyna 1992] https://doi.org/10.1063/1.463940\n\n[Patra 2014] https://doi.org/10.1103/PhysRevE.90.043304\n\n[Morrone 2011] https://doi.org/10.1063/1.3518369", "claims_and_evidence": "Most of the claims are supported by clear and convincing evidence. However, the following claim might be problematic to some extent:\n\n1. The paper made claims on the limitations of non-conservative force fields, whereas the empirical study focuses primarily on bulk water as the benchmark system. While water is a good test case, it remains unclear whether the same effects would be observed in other materials.", "methods_and_evaluation_criteria": "1. The authors conduct theoretical analysis and empirical validation on several benchmarks, including geometry optimization and molecular dynamics in NVE or NVT ensembles\n\n1. The choice of bulk water as a benchmark dataset is appropriate, as it is widely studied and sensitive to force inaccuracies. \n\n1. The authors use energy drift, Jacobian asymmetry, and kinetic temperature stability as the evaluation metrics which can directly capture the deficiencies of non-conservative models in simulations.", "theoretical_claims": "I checked the mathematical derivation for the symmetry of Jacobians and it is correct.", "experimental_designs_or_analyses": "I checked the soundness and validity of the experimental design and analysis in the MD simulations and geometry optimization tests. The experiments are well-designed to test the claims made in the paper:\n\n1. The paper compares multiple models, conservative or non-conservative, ensuring that findings are not tied to a single architecture.\n\n1. NVE simulations effectively demonstrate the unphysical energy drift in non-conservative force models. NVT simulations evaluate whether thermostats can alleviate this issue. The use of both Langevin and SVR thermostats is appropriate for evaluating different correction methods.\n\n1. The use of $\\lambda$ to quantify deviations from the conserved energy is a novel and effective metric.\n\n1. The manuscript provides the test results for different optimization methods (FIRE vs. L-BFGS) to demonstrate the impact of non-conservative forces on convergence.", "relation_to_broader_scientific_literature": "The paper reveals a deep insight into the limitations of non-conservative force fields which have been broadly studied in a lot of scientific literature for efficient force predictions. The paper points out that physics constraints should not be thoroughly removed, which provides a promising direction for the research in machine learning force fields.", "essential_references_not_discussed": "All the relevant and essential references are discussed. There are two recent papers that also discuss the limitations of non-conservative force models, though they may not be perfectly within the scope of this paper.\n\n[1] Loew, A. et al. Universal Machine Learning Interatomic Potentials are Ready for Phonons. Preprint at https://arxiv.org/abs/2412.16551v1.\n\n[2] Fu, X. et al. Learning Smooth and Expressive Interatomic Potentials for Physical Property Prediction. Preprint at https://arxiv.org/abs/2502.12147.", "other_strengths_and_weaknesses": "Strengths:\n\n1. The paper discusses the limitations of non-conservative force fields and points out that the best way to use non-conservative fore models is to integrate them into a hybrid model that follows physical constraints but also additionally supports direct, non-conservative,\nforce predictions. Therefore, instead of simply investigating the limitations, the paper provides a constructive pathway to mitigate the shortcomings of non-conservative force models, making them more applicable to real-world simulations.\n\nWeaknesses:\n\n1. The authors observe that thermostats cannot fully correct non-conservative behavior in NVT simulations, whereas only Langevin and SVR thermostats are tested. Exploring other thermostats (Nose-Hoover, adaptive thermostats, etc.) might provide a more complete picture.", "other_comments_or_suggestions": "See Strengths and Weaknesses.", "questions_for_authors": "1. Why is the water dataset chosen for benchmarking instead of other systems? Is there any particular chemical insight that we can tell from benchmarking on the water dataset?", "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": "", "extra_scores": {}}, {"review_id": "MnFIcMxj0c", "reviewer": "Reviewer_J2Ug", "summary": "This paper investigates the implications of using machine learning models that predict non-conservative forces for atomistic simulations. Non-conservative models predict interatomic forces directly, rather than computing them as the derivative of a potential energy, which offers computational advantages but violates the fundamental principle of energy conservation. The authors compare conservative and non-conservative models and assess their performance in various simulation scenarios, including constant-energy (NVE) and constant-temperature (NVT) molecular dynamics (MD), as well as geometry optimization. They demonstrate that non-conservative forces lead to significant issues, such as ill-defined convergence in geometry optimization, energy drift in NVE MD, and artifacts in NVT MD, even with strong thermostatting. Further, the authors propose a hybrid approach that combines a conservative model with a non-conservative force predictor, which combines the benefits of direct force prediction (increased efficiency) with the stability and physical correctness of conservative forces.\n\n## update after rebuttal\n\nPlease see comments below.", "questions": "", "limitations": "", "rating": 5, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the suggestion to investigate a solvated peptide, and have done so: we studied alanine dipeptide solvated in water, and additionally a benzene molecule adsorbed on a graphene surface, beyond the homogeneous cases of graphene, amorphous carbon and aluminium that we had already discussed in Appendix F. In all cases, the observation that non-conservative models like Equiformer and ORB display lack of energy conservation in NVE MD holds. The effect is comparable to or larger than that seen for homogeneous systems, and much larger than the drift observed for the more accurate single-purpose PET-NC model we trained for bulk water (https://imgur.com/a/CV3YmI3). \n\nAs we argue with Reviewer NcLF, our choice of a comparatively simple system, for which very accurate models can be trained, is aimed at giving non-conservative predictions the most favorable treatment possible, and these additional tests confirm that the artifacts become even more severe when using general-purpose, universal models that are less accurate and entail more pronounced violation of energy conservation.\n\nWe also thank the Reviewer for pointing out some typographical errors, which we have fixed.", "claims_and_evidence": "The central claim is that while non-conservative force models offer computational benefits, their use in atomistic simulations leads to fundamental problems due to the violation of energy conservation. This claim is well-supported by both theoretical arguments and empirical evidence.", "methods_and_evaluation_criteria": "The proposed methods and evaluation criteria are appropriate for the problem at hand and applied thoroughly.", "theoretical_claims": "The theoretical claims are correct and well-explained. The derivation of the non-conservation metric (λ) is mathematically sound, and other statements in the paper are based on well-known facts and established theory.", "experimental_designs_or_analyses": "The experimental design and analyses are all sound and well-executed. It would have been interesting to also investigate less homogenous systems beyond bulk water (e.g., a solvated peptide) to investigate how general trends transfer between different systems (e.g., are non-conservative forces a bigger issue for heterogenous systems?), but this is a minor nitpick.", "relation_to_broader_scientific_literature": "The paper clearly places the research in the context of current literature.", "essential_references_not_discussed": "The paper cites all essential references as far as I can judge.", "other_strengths_and_weaknesses": "**Strengths**\n\n+ The paper is very well-written and the arguments are presented logically and are easy to follow.\n+ The paper addresses an important issue in the field of machine learning potentials, as non-conservative models are becoming increasingly popular and many practitioners seem unaware of the downsides of such models.\n+ The authors provide a thorough investigation of the problem, covering both theoretical and practical aspects.\n+ The authors propose a hybrid approach that combines the advantages of conservative (accurate and reliable) and non-conservative (fast) models.\n+ The supplementary materials are very well constructed and provide many additional details to support the claims.\n\n**Weaknesses:**\n\n- The analysis is primarily focused on bulk water. While water is definitely a relevant system and I agree with its use as a paradigmatic example, it would be interesting to see how the findings generalize to other chemical systems. However, this is a very minor point (and the supplementary information also already contains some results for geometry optimizations of other systems).", "other_comments_or_suggestions": "I found two small typos in the supplementary information:\n* Line 912 on p.17 (\"without never\" should be \"without ever\")\n* Line 1028 on p.19 (\"langevin\" should be \"Langevin\")", "questions_for_authors": "I have no important questions for the 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": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: The authors present a rigorous analysis of several nonconservative machine-learned force fields, showing how the deviations from energy conservation can lead to severely undesirable behavior. The reviewers all agreed the paper was high-quality, valuable work, an appreciated the inclusion of additional experiments beyond bulk water. The recommendation to accept was unanimous, with many reviewers strongly in favor, and I concur.
Decision: Accept (oral)