Paper deep dive
Implicit Machine Learning Force Fields Accelerate Molecular Dynamics Simulations
Johannes Maeß, Leon Werner, J. Thorben Frank, Winfried Ripken, Martin Michajlow, Joshua Futterer, Klaus-Robert Müller, Stefan Chmiela
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 91%
Last extracted: 8/3/2026, 2:41:17 AM
Summary
The paper introduces Implicit Machine Learning Force Fields (I-MLFFs), a method that replaces explicit neural network stacks with self-consistent fixed-point equations. This approach allows for the reuse of intermediate representations across molecular dynamics (MD) timesteps, significantly reducing computational cost and memory footprint by 2-5 times compared to explicit models, while maintaining high accuracy and full atomistic resolution. The method is demonstrated across various Graph Neural Network architectures (SchNet, PaiNN, SO3net) and leverages regularization and warm-starting techniques to accelerate convergence.
Entities (10)
Relation Signals (9)
Implicit Machine Learning Force Fields → uses → Fixed-Point Equation
confidence 95% · We introduce implicit machine learning force fields (I-MLFFs), which replace explicit stacks of neural network layers with self-consistent fixed-point equations.
Implicit Machine Learning Force Fields → accelerates → Molecular Dynamics
confidence 94% · Implicit Machine Learning Force Fields Accelerate Molecular Dynamics Simulations
Implicit Machine Learning Force Fields → reduces → Compute Footprint
confidence 92% · Each yields a two- to five-fold reduction in compute and memory footprint.
Warm starting → enables → Implicit Machine Learning Force Fields
confidence 91% · enables intermediate representations to be reused across successive timesteps, thereby warm-starting force evaluation.
Implicit Machine Learning Force Fields → demonstratedon → SO3net
confidence 90% · We demonstrate this across three major classes of graph neural networks: ... SO(3)-equivariant spherical-tensor architectures.
Implicit Machine Learning Force Fields → demonstratedon → PaiNN
confidence 90% · We demonstrate this across three major classes of graph neural networks: ... PaiNN architecture
Implicit Machine Learning Force Fields → demonstratedon → SchNet
confidence 90% · We demonstrate this across three major classes of graph neural networks: ... SchNet architecture
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:We introduce implicit machine learning force fields (I-MLFFs), which replace explicit stacks of neural network layers with self-consistent fixed-point equations. In molecular simulations, this formulation enables intermediate representations to be reused across successive timesteps, thereby warm-starting force evaluation. The resulting models effectively combine the computational footprint of a shallow, single-layer MLFF with the representational capacity and accuracy of a deep neural network. Our approach unlocks architecture-agnostic efficiency gains that are inaccessible when force prediction and trajectory integration are considered separately. We demonstrate this across three major classes of graph neural networks: invariant, equivariant Cartesian tensor, and SO(3)-equivariant spherical-tensor architectures. Each yields a two- to five-fold reduction in compute and memory footprint. Crucially, these gains are achieved while retaining full atomistic resolution and the original integration timestep, avoiding spatial or temporal coarse graining. Our contribution therefore advances the scaling frontier of quantum-mechanically faithful molecular simulation, enabling longer trajectories and larger atomistic systems within fixed GPU memory and compute budgets, and thereby opening access to new insights across biomolecular and material systems.
Tags
Links
- Source: https://arxiv.org/abs/2607.29158v1
- Canonical: https://arxiv.org/abs/2607.29158v1
Trouble viewing inline? Open PDF directly →
Full Text
106,231 characters extracted from source content.
Expand or collapse full text
Implicit Machine Learning Force Fields Accelerate Molecular Dynamics Simulations Johannes Maeß, 1, 2 Leon Werner, 1, 2 J. Thorben Frank, 1, 2 Winfried Ripken, 1, 2 Martin Michajlow, 1, 2 Joshua Futterer, 1, 2 Klaus-Robert M ̈uller, 1, 2, 3, 4,∗ and Stefan Chmiela 1, 2,∗ 1 BIFOLD, Berlin, Germany. 2 Machine Learning Group, TU Berlin, Berlin, Germany. 3 Department of Artificial Intelligence, Korea University, Seoul, South Korea. 4 MPI for Informatics, Saarbr ̈ucken, Germany. We introduce implicit machine learning force fields (I-MLFFs), which replace explicit stacks of neural network layers with self-consistent fixed-point equations. In molecular simulations, this for- mulation enables intermediate representations to be reused across successive timesteps, thereby warm-starting force evaluation. The resulting models effectively combine the computational foot- print of a shallow, single-layer MLFF with the representational capacity and accuracy of a deep neural network. Our approach unlocks architecture-agnostic efficiency gains that are inaccessible when force prediction and trajectory integration are considered separately. We demonstrate this across three major classes of graph neural networks: invariant, equivariant Cartesian tensor, and SO(3)-equivariant spherical-tensor architectures. Each yields a two- to five-fold reduction in compute and memory footprint. Crucially, these gains are achieved while retaining full atomistic resolution and the original integration timestep, avoiding spatial or temporal coarse graining. Our contribu- tion therefore advances the scaling frontier of quantum-mechanically faithful molecular simulation, enabling longer trajectories and larger atomistic systems within fixed GPU memory and compute budgets, and thereby opening access to new insights across biomolecular and material systems. I. INTRODUCTION A major obstacle in accurate molecular dynamics (MD) simulations [1] of physical systems, such as large proteins in aqueous environments, is the need for mil- lions of computationally expensive quantum mechanical calculations to compute observables that may then be directly compared with experimental results. Over the past decade, machine learning force field (MLFF) devel- opment has been oriented toward this application. Mod- ern MLFFs are several orders of magnitude faster than ab-initio methods [2], while maintaining the accuracy of large-basis set density functional theory (DFT) [3– 12] and coupled-cluster reference calculations [13] for well-defined inference tasks.Remarkably, the fastest MLFFs available today achieve linear scaling [14]. While these models can achieve up to 10 6 simulation steps per day [15], they have still not reached the necessary effi- ciency for studies at desired practical scales. Many ex- perimentally relevant phenomena, such as protein fold- ing, occur over timescales that span microseconds to mil- liseconds, requiring billions to trillions of molecular dy- namics (MD) simulation steps [16]. Consequently, clas- sical FFs [17, 18] or coarse graining [19–23] remain the methods of choice for studying phenomena at the largest scales, even when the high accuracy of MLFFs would be preferable. Although continued advances in MLFF ar- chitectures and implementations will reduce the cost of individual force evaluations, conventional inference still repeats the full computation at every integration step. ∗ Corresponding authors. We therefore shift the focus from optimizing the MLFF in isolation to exploiting efficiencies across the MD sim- ulation loop as a whole. In atomistic simulations, the classical Newtonian equa- tions of motion are solved numerically to evolve nuclear positions over time. The computational bottleneck arises from the inherent need for fine temporal resolution, typ- ically on the order of one femtosecond (10 −15 s), to cap- ture high-frequency dynamics of light atoms accurately. This constraint reflects the Nyquist–Shannon sampling limit. Extending the timestep can lead to instabilities and unrealistic system behavior in the simulation [24]. However, at room temperature (∼300 K), the change in atomic positions during a single simulation step is typically on the order of only 0.001 ̊ A to 0.01 ̊ A (less than ∼1 % of a typical covalent bond length in organic molecules), leaving the global structure of the system mostly unchanged. As a result, successive evaluations of the MLFF along a simulation trajectory are highly similar. This calls for novel approaches that exploit this inherent redundancy to avoid starting the calculations from scratch at every step. To address this challenge, we propose using the prin- ciple of implicit modeling [25–28] to amortize inference across simulation time by systematically building on prior computations. Implicit models treat inference as the solution of an underlying dynamical system, effec- tively a continuous-time counterpart to traditional fixed- depth neural network (N) architectures. Deep Equilib- rium Models [27] are a specific formulation of implicit models in which inference is formulated as a convergent, self-consistent procedure whose limit corresponds to an equilibrium (fixed point) of a differential equation de- fined by a single learned layer [29–32]. In our setting arXiv:2607.29158v1 [cs.LG] 31 Jul 2026 2 h t * ( ) 3rd+ MD step: Linear extrapolation 1 iteration 1st MD step: Many iterations to converge 2nd MD step: Reuse of fixed-point 2-5 iterations h (t 2 ) (0) h (t 3 ) (0) h(t 1 ) * h(t 3 ) * f f f f f f h t * ( ) h (t 1 )(t 1 ) (0)(0) * h(t 2 ) t 3 t 2 t 1 Time h (0) u (0) h RF h (k) u (k) u * h * u h * Time t 3 t 2 t 1 (6) Warm- starts (3) Force prediction (4) Backward solver (1) Energy prediction (2) Forward solver Layers stored in memory Layers computed PaiNN SO3net SchNet Implicit (ours) Explicit Layers stored in memory 1 2 3 4 1 2 3 4 (c) Algorithmic flowchart of implicit MLFF (b) Fixed-point extrapolation accelerates force prediction (a) Continuity of fixed-points (d) Implicit MLFFs are 2-5x faster & smaller FIG. 1. Temporal continuity in MD enables computational reuse in implicit MLFFs. (a) Learned fixed-point representations h ∗ (t) of the implicit MLFF evolve smoothly in representation space along a C 2 H 4 torsion coordinate (black arrow ), with larger changes near bond breaking events. (b) This temporal continuity enables reuse and extrapolation of previous fixed points (red dashed arrows ,), reducing the number of solver iterations required in subsequent MD steps. The first step (bottom) requires many applications of f to converge, whereas later steps can be warm-started from previous fixed points using fixed-point reuse (middle) or linear extrapolation (top). (c) Algorithmic flowchart of an implicit MLFF. A forward fixed-point solve h ∗ = f (h ∗ , R) ( ↶ ↶ ↶ ) yields E = f E (h ∗ ). Forces, F = −∇ R E, are computed in the backward pass by solving a second fixed-point problem for u ∗ ( ↷ ↷ ↷ ). Both solves are warm-started from preceding MD steps. Corresponding equation numbers are indicated. (d) The computational footprint of implicit models is 2–5× smaller than explicit models with matching force accuracy, averaged across MD17 and MD22 datasets at 300 K. Implicit models using linear extrapolation need only 1–2 layer evaluations to predict forces, with implicit differentiation requiring only one layer in memory. In contrast, explicit models have to evaluate, store, and differentiate all their 3–5 layers in every MD step. this equilibrium defines a representation of the molecu- lar graph from which energies and atomic forces are pre- dicted. Because the solution is defined as a limit rather than a fixed sequence of computations, the solver can be initialized from arbitrary states. We use this flexibility to exploit the temporal coherence of simulation trajectories: intermediates from prior steps provide effective initializa- tions for subsequent solves. This enables warm-started inference, substantially reducing the number of iterations needed to converge the latent molecular graph represen- tation. Once the fixed point is reached, we further implic- itly differentiate the energy, avoiding the need to unroll the solver to obtain conservative forces. Together, this re- duces the inference cost per simulation step significantly, below that of traditional explicit MLFFs, bringing them closer to the practical applicability of classical FFs with- out compromising their characteristic accuracy. These gains are achieved without relaxing any aspect of physical accuracy, which sets apart our approach from many existing strategies for extending accessible simu- lation timescales, including prior work that exploits the temporal continuity of MD simulations [33, 34]. Coarse- grained methods, for example, accelerate simulations by sacrificing atomistic resolution [19–23]. Other ap- proaches predict molecular evolution by taking large ef- fective time steps [35–37], which can skip over high- frequency dynamics, such as bond vibrations or rapid hydrogen-bond rearrangements. A unique aspect of our solution is that it preserves energy conservation, a fun- damental property of isolated classical and quantum me- chanical systems [38]. In contrast, direct force-prediction methods [39–45] can introduce non-conservative force components that contribute to thermodynamic drift and runaway instabilities in molecular simulations [46, 47]. This limits them to short simulation times that are in- sufficient for realistic MD applications.While small energy-conservation violations can be corrected with thermostats, they disrupt the sampling efficiency of the trajectory [47]. I. RESULTS Many state-of-the-art MLFFs are built on Graph Neural Networks (GNNs), which encode the molecu- lar geometry into latent atomic representations by it- eratively passing messages along interatomic edges [2, 48]. Most explicit GNN-based MLFFs can be reformu- lated as implicit models with only minimal architectural changes. To demonstrate the general and architecture- independent applicability of our approach, we apply it to three GNN architectures that are considered repre- sentative of the main types of MLFFs: the invariant continuous-filter convolutional SchNet architecture [49], the Cartesian tensor message-passing PaiNN architec- ture [4], and a generic SO(3)-equivariant spherical-tensor MPNN, representing scalar, vectorial, and higher-order tensor feature spaces, respectively. The latter serves as a typical template for a broad class of contemporary equiv- ariant MLFF [7, 9–12, 42, 43, 50–53]. GNNs encode the molecular graph x := [R, Z], consisting of the nuclei po- sitions R ∈ R N×3 and atomic numbers Z ∈ N N + of N 3 atoms, into a latent representation h ∗ (x)∈ R N×F . From this a readout function f E infers the energy and forces: E = f E (h ∗ ) and F =−∇ R E.(1) Explicit models form the latent representation of the molecule by composing a stack of K interaction layers h (K) = f (K) ◦·◦f (1) (x), where each layer may have its own learned parameters. Notably, they offer no mecha- nism for adaptive compute or warm-starting, and all K layers need to be kept in memory in order to compute F with backpropagation. By contrast, implicit models iter- ate a single, learned interaction layer h (k+1) = f (h (k) , x) until the residual change in h is small, reaching a fixed point: h ∗ = f (h ∗ , x).(2) To guarantee the existence of a fixed point, f must map a compact convex set to itself [54], which is achieved through normalization. We further inject R and Z at every iteration to ensure that the fixed point remains conditioned on the molecular structure (Sec. IV A). By the Implicit Function Theorem (IFT) [55], equation (2) locally defines a differentiable mapping x7→ h ∗ (x). This enables efficient force computation by implicit differen- tiation (cf. Sec. IV B), expressing F entirely in terms of derivatives evaluated at h ∗ : F =− ∂f E ∂h ∗ ⊤ I− ∂f ∂h ∗ −1 | z (u ∗ ) ⊤ ∂f ∂R .(3) Because the gradients depend only on h ∗ , and not on the path taken to reach this fixed point, the forward solver can be chosen freely and does not require unrolling. Only the final evaluation of f must be retained in memory, reducing memory usage by up to a factor of K relative to an explicit K-layer network (cf. SI Fig. S6). This allows larger systems to be simulated on the same fixed-memory GPU. To avoid explicit matrix inversion, we reformulate the fixed-point adjoint u ∗ in equation (3) as the solution of a second fixed-point equation (Sec. IV B): u ∗ = ∂f ∂h ∗ ⊤ u ∗ + ∂f E ∂h ∗ .(4) We monitor convergence atom-wise using the residual max a∈atoms ∥□ (k−1) a −□ (k) a ∥ 2 /∥□ (k−1) a ∥ 2 for h and u, and find that a threshold of 10 −2 is sufficient for accurate force prediction (Sec. I A). Fast convergence of the fixed- point iteration is encouraged by a regularization scheme during training (Sec. IV C) and further accelerated in MD simulations by warmstarting the fixed-point solver with initial guesses h (0) ≈ h ∗ and u (0) ≈ u ∗ (which we de- velop in Sec. I B). The full algorithm for implicit energy- conservative FFs is illustrated in Fig. 1C and written out in Algorithm 1 (SI Sec. S7). Further intuition is provided in SI Sec. S11 through the mean-field Curie–Weiss Ising model, which reduces to a simple scalar implicit equation that is easier to interpret. A. Implicit MLFFs: Accuracy & Regularization Implicit MLFFs are governed by a trade-off between accuracy and cost: forces are only accurate once the latent fixed point is sufficiently converged, so accuracy is inherently coupled to the number of solver iterations required per prediction. In this section we study this trade-off and ask whether training-time regularization re- duces solver cost to the low, MD-compatible budgets that would make implicit models competitive. We begin by evaluating the accuracy of the implicit approach when fully converged. Implicit models reuse one parameter- ized layer across iterations, whereas explicit models grow parametric capacity with depth, which can potentially translate into higher accuracy. This leads to the ques- tion of whether parameter sharing fundamentally limits the implicit approach. To answer it, we compare a fully converged implicit model to explicit models with 1–8 lay- ers on the MD17 benchmark [38] (Fig. 2A), using PaiNN as the backbone architecture (see Sec. IV C for train- ing details). We also include a parameter-tied explicit variant (sharing parameters across layers) to decouple depth from parameter count in our analysis. Remarkably, the prediction performance of the implicit model exceeds both explicit baselines up to four layers, a commonly used depth for MLFFs [9, 10, 49, 56–58]. At greater depths, the explicit models do eventually achieve higher accu- racies, suggesting that the additional structure and con- straints introduced to make the implicit formulation well- defined may act as a limiting factor. Nevertheless, the strong performance at practically relevant depths makes implicit models a parameter-efficient alternative that pre- serves accuracy where it matters most. However, fully converging implicit models to machine precision, as in the previous test, is impractical: the re- quired iterations can exceed explicit depth and erase run- time gains. Implicit models are only useful if they yield accurate representations in a few iterations. Therefore, our second experiment quantifies the convergence toler- ances ε needed for the latent fixed-point embedding h ∗ to reach a target force-prediction accuracy (within ±2 % of the fully converged implicit F-MAE). We compare performance on the MD17-aspirin dataset across three representative MLFF architectures, SchNet, PaiNN, and SO3net (Fig. 2B). A key focus is the implicit backward pass used to compute forces efficiently (equation (4)), since the IFT relies on a sufficiently converged fixed point.If the forward iteration is stopped too early, gradients may become inconsistent and any computa- tional savings from early stopping could be offset. In practice, we observe that sensitivity is approximately symmetric with respect to the forward and backward tolerances. Restricting both to the same value there- fore yields efficient and accurate forces across all three architectures, a simplification that reduces the search space without sacrificing performance. Encouragingly, all models reach the target force accuracy at compar- atively loose convergence thresholds, with tolerances in 4 Residual 10 -1 10 -2 1 2 3 4 5 6 7 8 9 10 Iterations 10 -3 True trajectory Previous five Next fixed-point h t * ( ) t Warmstar No (6.4) Reuse (2.1) Linear (1.12) Deg=2 (1.02) Deg=3 (1.07) Deg=4 (1.3) (on I-PaiNN UL ) 4 5 6 7 8 9 10 Iteration 0.6 0.7 0.8 0.9 1.0 10 -2 10 -1 10 0 10 -3 1 2 3 4 5 6 7 8 9 Residual per iteration Regularization (on I-PaiNN UL ) No Force MAE (kcal/mol/A) I-PaiNN UL 10 -2 10 -1 10 0 Backward tolerance MD17 Ethanol Malonaldehyde Benzene Uracil Toluene Salicylic Acid Naphtalene Paracetamol Aspirin Azobenzene (5.0) (5.0) (9.2) (6.7) (7.4) (7.7) (8.2) (5.8) (6.1) (9.0) Aspirin 10 -2 10 -1 10 0 10 -2 10 -1 10 0 Backward tolerance Forward tolerance MD22 Ac-Ala3-NHMe DHA AT-AT Stachyose AT-AT-CG-CG Buckyball-Cat. (6.7) (8.7) (7.9) (6.5) (8.7) (9.0) (6.1) (7.4) (3.9) Architecture Force error of explicit relative to implicit models Explicit Explicit (tied) Mean, 95% conf. Layers of explicit MLFF ( Memory consumption) 2 2.5 3 1 1.5 0.75 1 2 3 4 5 6 7 8 Implicit Force MAE with 1 layer in memory (converged) I-PaiNN on MD17 (a) Accuracy-memory advantage (b) Required tolerance to converge F-MAE up to 2% (e) Sketch of warmstarts (d) Acceleration with warmstarts (c) Acceleration with regularization losses Aspirin 10 8 10 7 10 2 320 100 32 10 8 3.2 10 10 7 10 6 10 5 1 0.32 10 6 10 -1 10 5 10 4 10 3 10 4 10 10 8 0.032 10 3 10 -1 10 -2 1 No Doublewalled-N. (8.5) I-PaiNN UL I-SO3net UL I-SchNet UL FIG. 2. Understanding and implementing implicit MLFFs. (a) Implicit MLFFs with one layer match or exceed the accuracy of explicit models with up to four layers on MD17. We report average force MAEs of I-PaiNN and show the mean with 95 % confidence interval of the relative error increase (E− I)/I between explicit (E) and implicit (I) models. (b) Implicit models remain accurate with approximate solver convergence. Across three architectures (left) and 16 systems (right subplot), forward and backward solves require comparable tolerances of at most ∼10 −2 to converge the force error up to 2 %. Encouragingly, the necessary avg. iteration count (shown in brackets in the legends and in the marker size) does not generally increase with the molecule size. (c) Regularization can accelerate convergence without sacrificing accuracy. All four regularization losses reduce iteration counts, but only itc (γ=0.4) and jac preserve force accuracy. We show test MAEs versus the average iterations required to reach a 10 −2 residual, with the inset indicating convergence to this threshold (red line). (d) Warmstarts yield near one-layer implicit MLFFs. Linear extrapolation (equation (6)) accelerates convergence by 6× over no warmstart and 2× over fixed-point reuse, while preserving high accuracy. High-order extrapolation can further reduce iteration counts but may overfit. (e) Warmstart construction. Adams–Bashforth-inspired warmstarts extrapolate from the current fixed point using finite differences of previous fixed points to estimate derivatives. the range 10 −2 –10 −1 . These results suggest that early stopping is feasible across implicit formulations of dif- ferent architectures and largely insensitive to parameter settings. Whether a single acceleration technique can transfer broadly depends on how consistently these con- vergence characteristics hold across systems. To assess this, we repeat the test across molecules spanning differ- ent sizes, compositions, and conformational flexibilities, from the MD17 [38] and MD22 [59] benchmarks. Here, we focus on our implicit PaiNN variant (I-PaiNN), the more efficient of the two equivariant networks tested in the previous experiment (Fig. 2B-I). We find that for- ward and backward convergence criteria remain closely aligned across diverse systems. Convergence tolerance is weakly correlated with system flexibility (with the excep- tion of naphthalene and toluene): small, comparatively rigid aromatic molecules (e.g., salicylic acid, paraceta- mol, aspirin) tend to lie toward the upper right (looser tolerances), while multi-torsional and chain-like systems (e.g., Ac-Ala3-NHMe, stachyose, DHA) cluster some- what lower. Notably, we find no evidence that larger molecules require tighter convergence thresholds. For ex- ample, the longer DNA segment containing two Adenine– Thymine and two Cytosine–Guanine base pairs (AT-AT- CG-CG) converges in fewer iterations than the shorter DNA duplex segment AT-AT. These results suggest that the implicit approach remains robust to system complex- ity. Given that early stopping is feasible and consistent, we next ask whether training-time regularization can push 5 solver cost down to the low iteration counts MD re- quires. Because fixed points are easier to compute when the forward dynamics are simple, we consider three reg- ularizing losses to improve the stability of our implicit models [60] and identify which one is the most effec- tive: Jacobian regularization (jac) [61] adds a stochasti- cally estimated penalty on the Jacobian of f , effectively shrinking its spectral radius and encouraging contraction around the fixed point. Iterate correction (itc) explicitly pulls all earlier solver iterations h (k) toward h ∗ , encour- aging them to converge quickly. Finally, in truncated prediction (trunc), energies and forces are derived from the unrolled computation path of early representations h (1) , h (2) . The latter two regularization techniques were adapted from the approach introduced in Ref. [62] (see Sec. IV C). Fig. 2C quantifies the resulting trade-off be- tween convergence speed (average number of iterations to reduce the embedding residual below 10 −2 , on diverse conformations, cf. SI Sec. S9) and force accuracy (F- MAE). The unregularized implicit model (black dot) re- quires 10 iterations to reach the residual threshold, yield- ing an F-MAE of ∼0.62 kcal/mol/ ̊ A on the aspirin sys- tem. Increasing trunc strength accelerates convergence, but systematically degrades force accuracy. Applying itc with γ=1 yields even larger speed-ups, yet strong con- traction penalties over-regularize the dynamics and lead to substantially less accurate models. The variant itc with γ=0.4 focuses the convergence penalty on late it- erations, preserving accuracy. We find that jac provides the most favorable Pareto trade-off: it reduces the re- quired iterations from 10 to∼6 while preserving accuracy (∼0.6 kcal/mol/ ̊ A at a regularization of strength 0.32). We speculate that the Jacobian penalty improves local conditioning near the solution, accelerating convergence, while leaving the relative position of the fixed points and the function class represented by the model largely un- changed. Guided by these insights, models throughout the paper are trained using jac at a moderate strength of 0.32 to stabilize solver dynamics, and with itc (γ=0.4) at a mild strength of 10 3 to guarantee convergence within the 10 iterations allotted during training (Sec. IV C). Taken together, these experiments show that training- time regularization helps but plateaus: even efficient reg- ularizers require ∼6 calls of an interaction layer, short of the strict compute and accuracy demands required for MD. B. Warm-Starting Implicit MLFFs To improve this performance plateau, we exploit a source of structure unavailable to generic solvers: the temporal coherence along MD trajectories. The fixed point h ∗ (t) = h ∗ (R(t)) varies smoothly along the tra- jectory because the atomic positions R(t) evolve con- tinuously under Newton’s equations and the mapping R 7→ h ∗ (R) is differentiable by the IFT. Because MD uses small integration timesteps, consecutive simulation states h ∗ (t) , h ∗ (t+∆t) ,... are therefore close. This moti- vates the central objective of this work: warm-starting both the fixed-point and adjoint solves to accelerate the force prediction. In the simplest case, the fixed point from the previous step provides a natural initial guess for the next solve h (0) (t+∆t) :=h ∗ (t) , a zeroth-order (constant) approximation of h ∗ (t) used by prior work [34, 62]. We instead exploit the smoothness of the trajectory more carefully: to increase the accuracy of the warmstart fur- ther, we propose a first-order approximation (i.e., a linear extrapolation) T 1 across the MD time step ∆t: h ∗ (t + ∆t)≈ T 1 (t + ∆t) = h ∗ (t) + dh ∗ dt t ∆t.(5) We adopt a backward finite difference approximation for the derivative, dh ∗ dt t ∆t ≈ h ∗ (t) − h ∗ (t−∆t) , and initialize the fixed-point iteration in equation (2) at t + ∆t with: h (0) (t+∆t) := 2h ∗ (t) − h ∗ (t−∆t) .(6) However, to drive efficient MD with energy-conservative forces, not only the energy prediction in the forward pass but also the backward pass used for force derivation needs to be accelerated. A key methodological contribution of this work is the analogous warm-starting of the itera- tive solver used for implicit differentiation in equation (4) (Sec. IV B). This linear predictor generalizes to higher- order polynomial extrapolators using additional prede- cessors, similar to Adams–Bashforth multi-step integra- tors (SI Sec. S8). We note that higher orders can reduce iterations further but eventually risk overfitting to the local trajectory (Fig. 2E). This scheme incurs negligible overhead as it only relies on inexpensive combinations of already computed fixed points (SI Sec. S6). Fig. 2D evaluates the effectiveness of warmstarting simulations on the precomputed MD17-aspirin trajec- tory, averaged over both the forward and adjoint solve, with six different initializations: a baseline without warmstarts, the constant predictor (reusing the predeces- sor), the linear predictor (equation (6)), and polynomial predictors of degree k=2–4, using the k+1 predecessors, h ∗ (t−∆t) to h ∗ (t−(k+1)∆t) . We solve to the 10 −2 residual threshold established above on a broad set of conforma- tions as described in SI Sec. S9. We find that warm- starting consistently and substantially reduces iteration cost compared with the no-warm-start baseline. Notably, linear extrapolation converges after only 1.12 iterations on average, a striking six-fold speedup over the no-warm- start reference and an approximate two-fold improvement over simply reusing the preceding solution. We observe that most conformations converge after a single iteration, effectively reducing inference cost to that of a one-layer neural network. Additional iterations are invoked adap- tively only for more challenging configurations, such as those encountered in high-energy regions (cf. Fig. 4). 6 Mean Iterations (line) MD Stability (color) Tolerance 10 -3 10 -2 Tolerance 10 -3 10 -2 10 -1 Temperature Drift (K/fs) 10 -3 I-PaiNN Ac-Ala3-NHME, NVE without thermostat UL Lorem ipsum Iterations 1 2 3 4 5 6 7 8 1 2 3 1 2 3 MD17 .37 .33 .25 .2 .8 .76 MD22 .39 .26 .24 .6 .37 .32 1.38 1.27 1.25 . Layers/ Implicit F. MAE ( ) kcal mol A (a) Accuracy-compute advantage with warmstarts (c) Stability using thermostats (b) Temperature drift per tolerance PaiNN SO3net SchNet I-MLFF wins at low cost Implicit Force MAE (1-8 iter.) MD22 MD17 Mean, 95% conf. Force error of explicit relative to implicit models 2 2.5 1 1.5 FIG. 3. Accuracy–cost tradeoffs and stability of implicit MLFFs in MD simulations. (a) Implicit MLFFs outperform explicit baselines at matched computational cost, measured by interaction-layer calls averaged over MD17/MD22 trajectories. With linear fixed-point extrapolation from the two previous MD steps and solver tolerance 10 −2 , they reach most of their accuracy after one iteration and typically converge within two calls. We report mean force MAEs over 20,000 test samples and the relative error increase, (E− I)/I, between explicit (E) and implicit (I) models at equal average compute. (b) Energy drift in NVE simulations per solver tolerance. Plotted is the mean temperature drift with 95 % confidence intervals over 15× 50 ps MD trajectories initialized at 300 K run in NVE (without a regulating thermostat). Both (b) & (c) simulate Ac-Ala3-NHMe using an I-PaiNN UL model. (c) Stability of 50 ps MD simulations across solver tolerances and thermostat strengths. Background color gives the fraction of stable simulations; black lines and annotations show the mean implicit iteration count. Stronger thermostats suppress energy drift and broaden the stable tolerance range. C. Fast Molecular Dynamics: Accuracy & Cost In the following we will assess whether warmstarts can give implicit models an accuracy–compute advan- tage against their explicit counterparts in true MD sim- ulations driven by the I-MLFF. To make the compari- son independent of the force-field architecture and im- plementation, we measure cost as the number of calls to the interaction layer f , averaged over the forward and backward pass. Other contributions, such as the energy- readout head f E , are generally negligible and compara- ble across all models. For implicit models, the number of layer calls is controlled by the fixed-point solver tol- erance. Throughout, the solvers are warm-started with an initial guess obtained by linear extrapolation from the two previous steps as in equation (6). Fig. 3 re- ports the relative dataset-averaged force MAEs of im- plicit vs. explicit variants as a function of this compute for PaiNN, SO3net, and SchNet on the MD17 and MD22 datasets. Implicit models are consistently more accurate in the low-cost regime of 1−2 layer calls, often by a large margin. At a single layer iteration, I-PaiNN reaches a dataset-averaged force MAE of 0.37/0.39 (MD17/MD22) versus 0.75/0.84 for its explicit counterpart; I-SO3net scores 0.25/0.60 versus 0.76/0.70 (all in kcal/mol/ ̊ A). Beyond∼3 layer calls, explicit models become more accu- rate, but their advantage stays modest (at most ∼25 %) even at large layer counts, where they carry many times the parameters of the implicit model. Where explicit models first cross over implicit accuracy is architecture- and dataset-dependent: it occurs as early as two layers (PaiNN and SO3net on MD22) but only at five layers for the PaiNN architecture on MD17. D. Stability of Molecular Dynamics To be practically useful in MD, MLFFs must not only be efficient, but also remain stable over long simulation timescales. Accordingly, we next assess the stability of the implicit model on the Ac-Ala3-NHMe molecule while varying the tolerance of both solvers (Fig. 3B). Stabil- ity depends on the solver tolerance because incomplete convergence in the forward pass violates the assumptions required for implicit differentiation in the backward pass. As a result, the backward pass can yield forces that are inconsistent with the computed energy, which breaks en- ergy conservation and manifests as a spurious tempera- ture drift. We therefore quantify stability by the mean temperature drift (per fs) in microcanonical (NVE) tra- jectories run without a thermostat, so that any system- atic drift directly reports non-conservative force incon- sistencies. We average over 15 trajectories of 50 ps, each initialized at ∼300 K from a relaxed MD22 conformation (SI Sec. S10). Observed drifts are close to zero on fine tolerances and show a gradual increase for looser con- figurations. In practice, small drifts can be mitigated by running with a Langevin thermostat in the canonical (NVT) ensemble. The thermostat drives the system to- ward the target temperature on a characteristic Langevin coupling time (LTC). We validate for a range of solver tolerances which ther- mostat coupling strength is sufficient to regulate the in- duced temperature drift and keep the simulation stable (Fig. 3C). Across 15 independent MD trajectories we re- port the fraction of runs remaining stable by complet- ing the full 50 ps horizon without an unphysical bond- breaking event (any bonded interatomic distance exceed- ing twice its mean in the MD22 reference). Because overly strong coupling can distort dynamics and poten- 7 T=100K 200K 300K 400K 500K 600K 1 3 5 7 8 6 4 2 Mean solver iterations No warmstart Fixed-point reuse Linear extrapolation 14 12 10 8 6 4 Increased Increased computation computation at higher energiesat higher energies Increased computation at higher energies 1.1 1. NoNo trainingtraining datadata No training data Training Training datadata Training data (c) Iterations increase with temperature (Converging to a fine tolerance of 10 -3 ) (b) Iteration count: No warmstart (a) Iteration count: Linear extrapolationAc-Ala3-NHMe dihedral angles FIG. 4. Extrapolation behavior and solver dynamics of implicit MLFFs. (a,b) Ramachandran plots for Ac-Ala3- NHMe from implicit MLFF (I-PaiNN UL ). Although training conformations are largely confined to φ<0 ◦ (black contours), an MD run at 400 K for 1 ns samples an additional free-energy minimum outside this region. This indicates that implicit models preserve the extrapolative sampling behavior of their explicit counterparts. Color intensity visualizes mean fixed-point iteration counts in 10 ◦ (φ,ψ) bins for warm-started (a) and independently initialized (b) trajectories. Warm-starting keeps the solver close to the one-iteration limit across most conformations, while adaptively increasing iterations in high-energy states that likely have rapid changes in momenta. (c) Temperature dependence of iteration counts with cold- and warm-started solves on I-PaiNN UL . Iterations increase with temperature, but warm starts with linear extrapolation reduce this growth. Plotted are iterations needed to reach a fine tolerance of 10 −3 . Reaching 10 −2 takes ∼1 iteration using both warmstarts, (SI Fig. S5). tially bias the integrator, we focus on weak-to-moderate thermostats. Under these conditions, the thermostat provides sufficient regulation to offset the drift produced by the implicit model even at a relatively loose toler- ance of 10 −2 , enabling substantially faster simulations (fewer solver iterations) without compromising practical stability in NVT. Without a thermostat, NVE runs re- main stable up to a tolerance of 3.2× 10 −3 , averaging 2.2 solver iterations per step. A moderate coupling of LTC=1 ps lowers this to 1.7, and a stronger LTC=100 fs to just 1.2 iterations on average. Again, this brings the computational cost close to the practical lower bound of one solver iteration per step, even in this realistic MD setting. The cost-matched baseline fails this test: a one-layer explicit model never reaches reliable simula- tion, passing only 87 % of runs even at LTC=100 fs and 67 % in NVE. Implicit models therefore retain stability at near-minimal per-step cost, whereas the explicit coun- terparts do not. E. Generalization and Adaptive Iteration Depth Next, we examine whether reformulating an MLFF as an implicit model preserves its ability to generalize be- yond the training distribution. To visualize this, we first reproduce a generalization test from Ref. [12] on the flex- ible peptide Ac-Ala3-NHMe: a 1 ns MD simulation with I-PaiNN UL at 400 K correctly explores previously unseen regions of configuration space, consistent with the corre- sponding explicit model. The black contours overlaid on the Ramachandran plots in Fig. 4a,b show the density of the MD22 training set aggregated onto the two backbone angles φ,ψ of the molecule. The MD simulation driven by the implicit model discovers the chemically correct minimum [12] on the right side of the plot (φ>0), which is absent in the training data. We then use the same MD trajectory to show how I-MLFFs dynamically allocate computation across chem- ical space, concentrating effort where it is needed most. Fig. 4a,b visualizes the average iteration count (with and without warmstarts) in (10 ◦ ) 2 bins of the Ramachandran plot. We find that the simulation becomes highly efficient using warmstarts (Fig. 4a), with only a small fraction of MD steps taking more than one iteration to solve. Both plots highlight a distinctive advantage of implicit mod- els: the number of solver iterations adapts for chemically rare or out-of-distribution conformations, such as tran- sition pathways, and converges particularly quickly for frequent conformations near common energy minima, es- pecially when using warmstarts. Fig. 4c generalizes this finding by showing the effect of temperature on the number of iterations required to reach a residual threshold of 10 −3 , chosen here instead of the standard 10 −2 to better resolve differences in itera- tion count. We average over 10× 500 fs MD trajectories of Ac-Ala3-NHMe each at 100 K to 600 K (SI Sec. S9). Simulations at temperatures of ≥300 K increase the it- eration count of the cold-started solver, likely because they sample less inside well-known minima and more in the high-iteration-count regions visible in Fig. 4b. Warm- started solvers show a similar rise in iteration count with temperature: at low temperatures, small atomic displace- ments enable near-perfect amortization and one-step con- vergence, whereas at higher temperatures the previous fixed point becomes an increasingly poor initial guess. Linear extrapolation mitigates this effect by accounting for the evolution of h ∗ , thereby maintaining accurate 8 warmstarts and low solver cost. SI Sec. S4 provides ad- ditional evidence that high potential energy correlates with high iteration counts in the MD of Fig. 4 (a). SI Sec. S2 extends this analysis to another physical observ- able and shows an increased modeling fidelity of implicit models in high-frequency vibrational spectra of ethanol. Finally, SI Sec. S1 examines the generalization proper- ties of implicit models in greater depth on a dataset of cumulene structures C n H 4 (n = 2,..., 9). Their simple linear carbon-chain geometry provides a controlled set- ting to show that implicit models can accurately capture long-range electronic effects requiring up to ten message- passing steps, the longest range accessible in this dataset. Our results indicate that, given a sufficiently strong train- ing signal for long-range interactions, implicit modeling can overcome the finite effective receptive field that lim- its how far explicit models propagate geometric informa- tion. Our results further show that, when trained only on shorter cumulenes (n≤6), implicit models extrapo- late accurately to longer chains in both energy and force profiles, unlike explicit models. Together, these results suggest that implicit depth not only extends the range over which geometric information can be integrated, but also improves extrapolation beyond the molecular sizes seen during training. I. DISCUSSION & CONCLUSION Standard MLFF inference is stateless: at every fem- tosecond timestep during an MD simulation, the latent molecular representation needed for force prediction is re- constructed from scratch, even though successive atomic configurations differ only marginally. The implicit MLFF framework introduced here recasts this computation as a self-consistent fixed-point problem, enabling each force evaluation to be warm-started from the preceding so- lution and thereby amortized along the trajectory. We further improve these initial guesses by polynomially ex- trapolating previous solutions along the simulation tra- jectory. Our results show that a single interaction layer, iterated to convergence, matches the force accuracy of conventional explicit MLFFs that typically use up to five distinct layers and substantially more parameters. This advance thus directly addresses a key bottleneck in scaling MLFFs to larger systems, where speed and GPU memory are often the limiting constraints. Beyond atomistic simulation, the implicit modeling framework developed here may provide a route to ac- celerating forward inference and derivative computations in other physics-based ML applications that rely on iterative simulation. Examples include neural solvers for quantum-mechanical and electronic-structure prob- lems, such as neural wavefunctions [63, 64], learned or accelerated self-consistent-field methods [65–68], La- grangian [69] and Hamiltonian neural networks [70] (demonstrated in SI Sec. S12), and fluid models in which pressures, constraints, or energy gradients are obtained through differentiation. Related ideas also appear more broadly in generative modeling: energy-based models, diffusion and score-based models [71–73], including their probability-flow ODE formulations [30], as well as con- tinuous normalizing flows [74, 75], all perform inference by integrating learned forces, scores, or vector fields and could thus profit from our framework. In summary, our results establish temporal continu- ity in molecular trajectories as a computational resource that can be exploited directly through implicit model design. Across three GNN architectures, two molecu- lar benchmarks, and systems spanning nearly two or- ders of magnitude in atom count, our I-MLFF framework translates into a consistent two- to five-fold simulation speedup and smaller memory footprint. Because it builds directly on existing GNN architectures, its improvements are orthogonal to (future) architectural advances, im- plementation optimizations, and improved training-data sampling strategies. IV. METHODS A. Constructing Implicit MLFFs In conventional GNN-based (explicit) MLFFs, the molecular representation h used to predict the potential energy and atomic forces (equation. 1) is constructed by successively refining atom-wise features. Starting from learned embeddings of the atomic types, h (0) := h Z , the representation is updated through a fixed sequence of distinct interaction layers, h (k) = Interact (k) R (h (k−1) ), each of which incorporates geometric information from the nuclear positions. In contrast, our implicit MLFF approach predicts these outputs from a self-consistent molecular representation, h ∗ = f (h ∗ , x), defined as the fixed point of a single learned interaction layer f that de- pends implicitly on x := [R, Z]. To convert an existing explicit MLFF into the implicit form, we introduce two small modifications that ensure the existence of such a fixed point while preserving the general structure of the original architecture: h (k) = f (h (k−1) , x) := Norm◦ Interact R h (k−1) + h Z . (7) First, we add a dependence on h Z into the interaction layer, in order to make the fixed point a function of the atomic numbers, and not only of the atomic posi- tions R. This is essential in the implicit setting, where the converged representation no longer retains informa- tion from the initial state h (0) used to start the solver. Second, we test two different normalizations similar to those being used in explicit models to enable depth scal- ing, to bound the interaction layer outputs. A sim- ple equivariant implementation of the root mean square norm [76] rescales h on each atom a and each of its fea- ture blocks b ∈ invariant, equivariant, independently 9 to unit length, Norm UL = h a,b /(∥h a,b ∥ 2 + η), with a small positive η = 10 −5 ensuring continuity. The more expressive equivariant merged layer norm [77] (Norm LN ) first subtracts the mean from the invariant features, and then divides all invariant and equivariant features per atom by their merged norm. It then scales each invari- ant and equivariant feature by a learned factor and adds a learned offset on each invariant feature. This increases the force accuracy of implicit models by around 20 % on most molecules (SI Table S3). The main results in Fig. 2a and Fig. 3a report accuracy using Norm LN models; qual- itative results using the simpler norm are denoted as e.g. I-PaiNN UL . SI Sec. S5 highlights accuracy differences between the two norms per molecule. Importantly, both normalizations make f a continuous self-map on a com- pact space, guaranteeing the existence of a fixed point by Brouwer’s fixed-point theorem [54]. B. Implicit Differentiation And Force Inference At each simulation step,we compute energy- conservative forces by implicitly differentiating the con- verged fixed point. Because implicit differentiation de- pends only on local derivative information at this fixed point, it avoids backpropagation through the full solver trajectory. The resulting forces are therefore independent of the warm start used to initialize the solver. Moreover, derivative computation requires retaining only the final application of f , so the memory cost matches that of a single-layer explicit model and does not grow with the number of fixed-point iterations used to obtain h ∗ . The implicit Jacobian dh ∗ dR in equation (3) is derived by differentiating the fixed-point equation h ∗ (x) = f (h ∗ (x), x) with respect to the atomic position R and collecting terms: dh ∗ dR = ∂f ∂h ∗ dh ∗ dR + ∂f ∂R ⇔ dh ∗ dR − ∂f ∂h ∗ dh ∗ dR = ∂f ∂R ⇔ I− ∂f ∂h ∗ dh ∗ dR = ∂f ∂R ⇔ dh ∗ dR = I− ∂f ∂h ∗ −1 ∂f ∂R In the definition of the interatomic forces in equation (1), the inverse term is absorbed into the fixed-point adjoint u ∗ : −F = ∂f E ∂h ∗ ⊤ dh ∗ dR = ∂f E ∂h ∗ ⊤ I− ∂f ∂h ∗ −1 | z (u ∗ ) ⊤ ∂f ∂R = (u ∗ ) ⊤ ∂f ∂R We then redefine u ∗ as a linear self-consistency equation that can be solved iteratively: u ∗ = I− ∂f ∂h ∗ −⊤ ∂f E ∂h ∗ ⇔ I− ∂f ∂h ∗ ⊤ u ∗ = ∂f E ∂h ∗ ⇔u ∗ − ∂f ∂h ∗ ⊤ u ∗ = ∂f E ∂h ∗ ⇔u ∗ = ∂f ∂h ∗ ⊤ u ∗ + ∂f E ∂h ∗ C. Training Implicit MLFFs We train all implicit models without warmstarts, using h (0) = h Z as initialization and unrolling ten iterations of the computation graph to h (10) . We find that this method efficiently approximates a converged fixed point if the residual between h (10) and h (9) is smaller than 0.01. Energy and force predictions are trained jointly with a mean-squared-error (MSE) loss over the n atoms of each system, weighting the force contribution by α =95 %: L E,F = (1−α)·|E−E True | 2 +α/n·||F− F True || 2 2 . (8) a. Regularization losses. We test three additional loss functions to encourage fast convergence of implicit solvers (see Sec. I A): Jacobian Regularization (jac) [61] promotes a smaller spectral radius ρ( ∂f ∂h ∗ ), thereby im- proving convergence in equations (2) and (4). In practice, the spectral radius is upper bounded by the Frobenius norm, which is estimated following Girard [78], L jac = ε ⊤ ∂f ∂h h ∗ 2 2 , with ε∼N (0,I d ). A grid search on the aspirin dataset in Sec. I A showed that applying Jacobian regularization with a moderate coefficient reduces iteration counts without compromis- ing force accuracy. Accordingly, we train all implicit models in this study with L jac with a moderate coeffi- cient of 0.32. Another loss, ‘truncated prediction’ (trunc), explicitly derives E (k) =f E (h (k) ) and F=−∇ R E (k) from the early representations h (1) and h (2) and minimizes their dis- tance to labels, akin to equation (8). Unlike fixed-point correction [62], we limit supervision to early iterates – each iterate’s force loss requires second-order backprop- agation, adding the compute and memory footprint of k interaction layers per training step. To avoid this over- head, we propose the lightweight ‘iterate correction’ (itc) which does not fit force labels but minimizes the Eu- clidean distance of all solver iterates to the final repre- sentation. In full, L γ itc = P k γ 10−k ||h (k) − h (10) || 2 2 . Here, we let gradients propagate only through h (k) , not h (10) , 10 and can exponentially downweight the loss on early iter- ates by setting γ<1. The grid search (Sec. I A) shows that both trunc and itc with γ=1 (weighting all iterates equally) reduce iteration counts at the cost of lower pre- dictive accuracy. Setting γ=0.4 focuses itc on the conver- gence of late iterates. We apply this variant with a mild loss coefficient of 10 4 on all following implicit models as it guarantees convergence of h (k) within the 10 iterations allotted during training without compromising accuracy. All losses are weighted by their coefficients and summed, giving L =L E,F + 0.32·L jac + 10 4 ·L γ=0.4 itc . b. Convergence criteria. All models are trained with the AdamW optimizer [79] until the loss L on the validation set shows no improvement for 500 epochs. The learning rate is initialized to 10 −3 and is halved whenever the validation loss has not improved for 250 epochs. D. MLFF Architectures We use the implementations of SchNet, PaiNN, and SO3net provided by SchNetPack [80, 81] for training, testing, and MD simulations. Our extensions comprise input injection and normalization for all three archi- tectures, the addition of the fixed-point solver, implicit differentiation, warm-start logic, and the regularization losses described above. The three analyzed architectures encode molecular ge- ometry at increasing fidelity. SchNet acts purely on in- variant features and views the geometry of the system only through the pair-wise distances between atoms [56]. PaiNN represents atomic neighborhoods by scalar (l = 0) and Cartesian equivariant vector (l = 1) features, while SO3net uses spherical-harmonic features up to degree L max =2. The latter is representative of a broad class of modern equivariant architectures [7–12]. All models are configured with a cutoff of 5 ̊ A, 50 radial basis functions, and the SiLU activation function [82]. E. Datasets We train explicit MLFFs of varying depth and implicit models with one interaction layer on each molecule of the MD17 and MD22 datasets. MD17 contains trajec- tories of 10 small molecules containing 9 to 24 atoms each. Reference energies and forces are calculated with DFT at the pbe+vdw-ts standard, driving simulations with a timestep of 0.5 fs, at a temperature of 500 K [38]. MD22 contains seven biomolecules and supra-molecules ranging from 42 up to 370 atoms. Each trajectory is sam- pled at a time resolution of 1 fs at temperatures ranging from 400 K to 500 K; energies and forces are calculated at the pbe+mbd standard [59]. We use the same random dataset splits across training runs of different models. On molecules of the MD17 dataset we use 950 training and 50 validation conformations and increase to 5000 and 1000 on the more complex systems in MD22. Because of the lower number of available conformations, we reduce training points to 3000 on the ‘buckyball catcher’ and ‘double-walled nanotube’ systems. Test errors through- out the paper refer to mean absolute force errors of the remaining conformations in kcal/mol/ ̊ A. The sequential nature of the MD trajectories in both datasets allows us to measure the precision achieved with warmstarts and a limited compute budget: Fig. 3 mea- sures this on 20,000 randomly chosen test conformations. First, their predecessor conformations in the reference trajectory are solved by the implicit model to a tolerance that is realistic in a production MD (10 −2 ). The prede- cessor fixed points are then used to predict the fixed point of the test conformation, analogously to our method dur- ing MD. Finally, a scan over possible iteration depths gauges the accuracy-speed tradeoff of the warmstarted implicit model. F. Code Availability The source code for training, testing, and running molecular dynamics simulations with I-MLFF can be found at https://github.com/johannesmaess/imlff. V. ACKNOWLEDGMENTS JM, LW, JTF, WR, M, KRM, and SC acknowl- edge support by the German Federal Ministry of Re- search, Technology and Space (BMFTR) under Grants BIFOLD24B, BIFOLD25B, 01IS18037A, 01IS18025A, and 01IS24087C. KRM was partly supported by the Insti- tute of Information & Communications Technology Plan- ning & Evaluation (IITP) grants funded by the Korea government (MSIT) (No.2019-0-00079, Artificial Intelli- gence Graduate School Program, Korea University and No. 2022-0-00984, Development of Artificial Intelligence Technology for Personalized Plug-and-Play Explanation and Verification of Explanation) and also by DFG and HFA. KRM was furthermore supported in parts by the Korea University Grant. Correspondence to KRM and SC. We thank Julia Henkel, Alin Banka, Tim Ebert, and Gregor Lied for continuous feedback on our research. JM thanks Ana Runji ́c, Saskia ̈ Ozt ̈urk, Rebekka Maeß, Aslıhan Y ̈uksel, Nicole Trappe, and Jan Vincent Szlang for support, design, and proofreading of the manuscript. 11 [1] M. Karplus and J. A. McCammon, Molecular dynam- ics simulations of biomolecules, Nat. Struct. Biol. 9, 646 (2002). [2] O. T. Unke, S. Chmiela, H. E. Sauceda, M. Gastegger, I. Poltavsky, K. T. Sch ̈utt, A. Tkatchenko, and K.-R. M ̈uller, Machine learning force fields, Chem. Rev. 121, 10142 (2021). [3] J. Gasteiger, J. Groß, and S. G ̈unnemann, Directional message passing for molecular graphs, in ICLR (2020). [4] K. Sch ̈utt, O. Unke, and M. Gastegger, Equivariant mes- sage passing for the prediction of tensorial properties and molecular spectra, in ICML (PMLR, 2021) p. 9377– 9388. [5] T. W. Ko, J. A. Finkler, S. Goedecker, and J. Behler, A fourth-generation high-dimensional neural network po- tential with accurate electrostatics including non-local charge transfer, Nat. commun. 12, 398 (2021). [6] J. Gasteiger, F. Becker, and S. G ̈unnemann, Gem- Net: Universal directional graph neural networks for molecules, in NeurIPS, Vol. 34 (2021) p. 6790–6802. [7] Y.-L. Liao and T. Smidt, Equiformer: Equivariant graph attention transformer for 3D atomistic graphs, in ICLR (2023). [8] Y. Wang, S. Li, X. He, M. Li, Z. Wang, N. Zheng, B. Shao, T.-Y. Liu, and T. Wang, ViSNet: an equivariant geometry-enhanced graph neural network with vector- scalar interactive message passing for molecules, arXiv preprint arXiv:2210.16518 (2023). [9] I. Batatia, D. P. Kov ́acs, G. Simm, C. Ortner, and G. Cs ́anyi, MACE: Higher order equivariant message passing neural networks for fast and accurate force fields, in NeurIPS, Vol. 35 (2022) p. 11423–11436. [10] S. Batzner, A. Musaelian, L. Sun, M. Geiger, J. P. Mailoa, M. Kornbluth, N. Molinari, T. E. Smidt, and B. Kozinsky, E(3)-equivariant graph neural networks for data-efficient and accurate interatomic potentials, Nat. Commun. 13, 2453 (2022). [11] A. Musaelian, S. Batzner, A. Johansson, L. Sun, C. J. Owen, M. Kornbluth, and B. Kozinsky, Learning local equivariant representations for large-scale atomistic dy- namics, Nat. Commun. 14, 579 (2023). [12] J. T. Frank, O. T. Unke, K.-R. M ̈uller, and S. Chmiela, A Euclidean transformer for fast and stable machine learned force fields, Nat. Commun. 15, 6539 (2024). [13] S. Chmiela, H. E. Sauceda, K.-R. M ̈uller, and A. Tkatchenko, Towards exact molecular dynamics simu- lations with machine-learned force fields, Nat. Commun. 9, 3887 (2018). [14] A.Kabylda,J.T.Frank,S.Su ́arez-Dou, A. Khabibrakhmanov, L. Medrano Sandonas, O. T. Unke, S. Chmiela, K.-R. M ̈uller, and A. Tkatchenko, Molecular simulations with a pretrained neural network and universal pairwise force fields, J. Am. Chem. Soc. 147, 33723 (2025). [15] D. P. Kov ́acs, J. H. Moore, N. J. Browning, I. Batatia, J. T. Horton, V. Kapil, W. C. Witt, I.-B. Magd ̆au, D. J. Cole, and G. Cs ́anyi, MACE-OFF23: Transferable ma- chine learning force fields for organic molecules, arXiv preprint arXiv:2312.15211 (2023). [16] M. E. Tuckerman, Ab initio molecular dynamics: basic concepts, current trends and novel applications, J. Phys.: Condens. Matter 14, R1297 (2002). [17] W. D. Cornell, P. Cieplak, C. I. Bayly, I. R. Gould, K. M. Merz, D. M. Ferguson, D. C. Spellmeyer, T. Fox, J. W. Caldwell, and P. A. Kollman, A second generation force field for the simulation of proteins, nucleic acids, and organic molecules, J. Am. Chem. Soc. 117, 5179 (1995). [18] A. D. MacKerell, D. Bashford, M. Bellott, R. L. Dunbrack, J. D. Evanseck, M. J. Field, S. Fischer, J. Gao, H. Guo, S. Ha, D. Joseph-McCarthy, L. Kuch- nir, K. Kuczera, F. T. K. Lau, C. Mattos, S. Mich- nick, T. Ngo, D. T. Nguyen, B. Prodhom, W. E. Rei- her, B. Roux, M. Schlenkrich, J. C. Smith, R. Stote, J. Straub, M. Watanabe, J. Wiorkiewicz-Kuczera, D. Yin, and M. Karplus, All-atom empirical potential for molecular modeling and dynamics studies of proteins, J. Phys. Chem. B 102, 3586 (1998). [19] J. Wang, S. Olsson, C. Wehmeyer, A. P ́erez, N. E. Charron, G. de Fabritiis, F. No ́e, and C. Clementi, Ma- chine learning of coarse-grained molecular dynamics force fields, ACS Cent. Sci. 5, 755 (2019). [20] B. E. Husic, N. E. Charron, D. Lemm, J. Wang, A. P ́erez, M. Majewski, A. Kr ̈amer, Y. Chen, S. Olsson, G. de Fab- ritiis, F. No ́e, and C. Clementi, Coarse graining molecu- lar dynamics with graph neural networks, J. Chem. Phys. 153, 194101 (2020). [21] M. Majewski, A. P ́erez, P. Th ̈olke, S. Doerr, N. E. Char- ron, T. Giorgino, B. E. Husic, C. Clementi, F. No ́e, and G. De Fabritiis, Machine learning coarse-grained poten- tials of protein thermodynamics, Nat. Commun. 14, 5739 (2023). [22] N. E. Charron, K. Bonneau, A. S. Pasos-Trejo, A. Gul- jas, Y. Chen, F. Musil, J. Venturin, D. Gusew, I. Za- porozhets, A. Kr ̈amer, C. Templeton, A. Kelkar, A. E. P. Durumeric, S. Olsson, A. P ́erez, M. Majewski, B. E. Hu- sic, A. Patel, G. De Fabritiis, F. No ́e, and C. Clementi, Navigating protein landscapes with a machine-learned transferable coarse-grained model, Nat. Chem. 17, 1284 (2025). [23] A. E. P. Durumeric, Y. Chen, A. S. Pasos-Trejo, F. No ́e, and C. Clementi, Learning data-efficient coarse-grained molecular dynamics from forces and noise, Nat. Com- mun. 17, 2493 (2026). [24] T. J. Lane, D. Shukla, K. A. Beauchamp, and V. S. Pande, To milliseconds and beyond: challenges in the simulation of protein folding, Curr. Opin. Struct. Biol. 23, 58 (2013). [25] P. Simard, M. Ottaway, and D. Ballard, Fixed point anal- ysis for recurrent networks, NeurIPS 1, 149 (1988). [26] J. Miller and M. Hardt, Stable recurrent models, in Inter- national Conference on Learning Representations (2019). [27] S. Bai, J. Z. Kolter, and V. Koltun, Deep equilibrium models, in NeurIPS, Vol. 32 (Curran Associates, Inc., 2019) p. 690–701. [28] E. Winston and J. Z. Kolter, Monotone operator equilib- rium networks, in NeurIPS, Vol. 33 (Curran Associates, Inc., 2020) p. 10718–10728. [29] Y. Lu, A. Zhong, Q. Li, and B. Dong, Beyond finite layer neural networks: Bridging deep architectures and numer- ical differential equations, in ICML (PMLR, 2018) p. 3276–3285. 12 [30] R. T. Q. Chen, Y. Rubanova, J. Bettencourt, and D. K. Duvenaud, Neural ordinary differential equations, in NeurIPS, Vol. 31 (Curran Associates, Inc., 2018) p. 6571–6583. [31] E. Haber and L. Ruthotto, Stable architectures for deep neural networks, Inverse Probl. 34, 014004 (2017). [32] L. Ruthotto and E. Haber, Deep neural networks moti- vated by partial differential equations, J. Math. Imaging. Vis. 62, 352 (2020). [33] L. L. Schaaf, I. Batatia, J. Tilly, and T. D. Barrett, BoostMD: Accelerated molecular sampling leveraging ml force field features, in NeurIPS 2024 Workshop on Data- driven and Differentiable Simulations, Surrogates, and Solvers (2024). [34] A. Burger, L. Thiede, A. Aspuru-Guzik, and N. Vijayku- mar, DEQuify your force field: Towards efficient simula- tions using deep equilibrium models, in AI for Accelerated Materials Design Workshop, ICLR (2025). [35] F. L. Thiemann, T. Resch ̈utzegger, M. Esposito, T. Tad- dese, J. D. Olarte-Plata, and F. Martelli, Force-free molecular dynamics through autoregressive equivariant networks, arXiv preprint arXiv:2503.23794 (2025). [36] F. Bigi, S. Chong, A. Kristiadi, and M. Ceriotti, FlashMD: long-stride, universal prediction of molecular dynamics, in NeurIPS (2026). [37] W. Ripken, M. Plainer, G. Lied, T. Frank, O. T. Unke, S. Chmiela, F. No ́e, and K.-R. M ̈uller, Learning hamilto- nian flow maps: Mean flow consistency for large-timestep molecular dynamics, arXiv preprint arXiv:2601.22123 (2026). [38] S. Chmiela, A. Tkatchenko, H. E. Sauceda, I. Poltavsky, K. T. Sch ̈utt, and K.-R. M ̈uller, Machine learning of ac- curate energy-conserving molecular force fields, Sci. Adv. 3, e1603015 (2017). [39] W. Hu, M. Shuaibi, A. Das, S. Goyal, A. Sriram, J. Leskovec, D. Parikh, and C. L. Zitnick, ForceNet: A graph neural network for large-scale quantum calcula- tions, arXiv preprint arXiv:2103.01436 (2021). [40] C. L. Zitnick, A. Das, A. Kolluru, J. Lan, M. Shuaibi, A. Sriram, Z. Ulissi, and B. Wood, Spherical channels for modeling atomic interactions, in NeurIPS, Vol. 35 (2022). [41] J. Gasteiger, M. Shuaibi, A. Sriram, S. G ̈unnemann, Z. Ulissi, C. L. Zitnick, and A. Das, GemNet-OC: Devel- oping graph neural networks for large and diverse molec- ular simulation datasets, Transactions on Machine Learn- ing Research 10.48550/arXiv.2204.02782 (2022). [42] S. Passaro and C. L. Zitnick, Reducing SO(3) convolu- tions to SO(2) for efficient equivariant GNNs, in ICML, Vol. 202 (PMLR, 2023) p. 27420–27438. [43] Y.-L. Liao, B. M. Wood, A. Das, and T. Smidt, EquiformerV2: Improved equivariant transformer for scaling to higher-degree representations, in ICLR (2024). [44] M. Neumann, J. Gin, B. Rhodes, S. Bennett, Z. Li, H. Choubisa, A. Hussey, and J. Godwin, Orb:A fast, scalable neural network potential, arXiv preprint arXiv:2410.22570 (2024). [45] M. Eissler, T. Korjakow, S. Ganscha, O. T. Unke, K.-R. M ̈uller, and S. Gugler, How simple can you go? an off- the-shelf transformer approach to molecular dynamics, J. Chem. Phys. 164 (2026). [46] X. Fu, Z. Wu, W. Wang, T. Xie, S. Keten, R. Gomez- Bombarelli, and T. Jaakkola, Forces are not enough: Benchmark and critical evaluation for machine learn- ing force fields with molecular simulations, Trans. Mach. Learn. Res. (2023). [47] F. Bigi, M. F. Langer, and M. Ceriotti, The dark side of the forces: assessing non-conservative force models for atomistic machine learning, in ICML, Proceedings of Ma- chine Learning Research, Vol. 267, edited by A. Singh, M. Fazel, D. Hsu, S. Lacoste-Julien, F. Berkenkamp, T. Maharaj, K. Wagstaff, and J. Zhu (PMLR, 2025) p. 4384–4414. [48] J. A. Keith, V. Vassilev-Galindo, B. Cheng, S. Chmiela, M. Gastegger, K.-R. M ̈uller, and A. Tkatchenko, Com- bining machine learning and computational chemistry for predictive insights into chemical systems, Chem. Rev. 121, 9816 (2021). [49] K. T. Sch ̈utt, H. E. Sauceda, P.-J. Kindermans, A. Tkatchenko, and K.-R. M ̈uller, SchNet–a deep learn- ing architecture for molecules and materials, J. Chem. Phys. 148 (2018). [50] N. Thomas, T. Smidt, S. Kearnes, L. Yang, L. Li, K. Kohlhoff, and P. Riley, Tensor field networks: Rotation- and translation-equivariant neural networks for 3D point clouds, arXiv preprint arXiv:1802.08219 (2018). [51] O. T. Unke, S. Chmiela, M. Gastegger, K. T. Sch ̈utt, H. E. Sauceda, and K.-R. M ̈uller, SpookyNet: Learning force fields with electronic degrees of freedom and nonlo- cal effects, Nat. Commun. 12, 7273 (2021). [52] X. Fu, B. M. Wood, L. Barroso-Luque, D. S. Levine, M. Gao, M. Dzamba, and C. L. Zitnick, Learning smooth and expressive interatomic potentials for physical prop- erty prediction, in ICML, Vol. 267 (2025) p. 17875– 17893. [53] B. M. Wood, M. Dzamba, X. Fu, M. Gao, M. Shuaibi, L. Barroso-Luque, K. Abdelmaqsoud, V. Gharakhanyan, J. R. Kitchin, D. S. Levine, K. Michel, A. Sriram, T. Co- hen, A. Das, A. Rizvi, S. J. Sahoo, Z. W. Ulissi, and C. L. Zitnick, UMA: A family of universal models for atoms, in NeurIPS (2026). [54] L. E. J. Brouwer, ̈ Uber abbildung von mannigfaltigkeiten, Mathematische Annalen 71, 97 (1911). [55] S. G. Krantz and H. R. Parks, The implicit function the- orem: history, theory, and applications (Springer Science & Business Media, 2002). [56] K. Sch ̈utt, P.-J. Kindermans, H. E. Sauceda Felix, S. Chmiela, A. Tkatchenko, and K.-R. M ̈uller, Schnet: A continuous-filter convolutional neural network for mod- eling quantum interactions, in NeurIPS , Vol. 30 (Curran Associates, Inc., 2017) p. 991–1001. [57] A. Grisafi, A. Fabrizio, B. Meyer, D. M. Wilkins, C. Corminboeuf, and M. Ceriotti, Equivariant graph neural networks for fast electron density estimation of molecules, liquids, and solids, npj Comput. Mater. 8, 183 (2022). [58] M. Esders, T. Schnake, J. Lederer, A. Kabylda, G. Mon- tavon, A. Tkatchenko, and K.-R. M ̈uller, Analyzing atomic interactions in molecules as learned by neural net- works, J. Chem. Theory Comput. 21, 714 (2025). [59] S. Chmiela,V. Vassilev-Galindo,O. T. Unke, A. Kabylda, H. E. Sauceda, A. Tkatchenko, and K.-R. M ̈uller, Accurate global machine learning force fields for molecules with hundreds of atoms, Sci. Adv. 9, eadf0873 (2023). [60] K. Kawaguchi, On the theory of implicit deep learning: Global convergence with implicit layers, in ICML (2021). 13 [61] S. Bai, V. Koltun, and Z. Kolter, Stabilizing equilib- rium models by jacobian regularization, in ICML, Vol. 139 (PMLR, 2021) p. 554–565. [62] S. Bai, Z. Geng, Y. Savani, and J. Z. Kolter, Deep equi- librium optical flow estimation, in Proceedings of the IEEE/CVF conference on computer vision and pattern recognition (2022) p. 610–620. [63] J. S. Spencer, D. Pfau, A. Botev, and W. M. C. Foulkes, Better, faster fermionic neural networks, arXiv preprint arXiv:2011.07125 (2020). [64] J. Hermann, Z. Sch ̈atzle, and F. No ́e, Deep-neural- network solution of the electronic Schr ̈odinger equation, Nat. Chem. 12, 891 (2020). [65] K. T. Sch ̈utt, M. Gastegger, A. Tkatchenko, K.-R. M ̈uller, and R. J. Maurer, Unifying machine learning and quantum chemistry with a deep neural network for molec- ular wavefunctions, Nat. Commun. 10, 5024 (2019). [66] F. Song and J. Feng, Neural network self-consistent fields for density functional theory, npj Comput. Mater. 10.1038/s41524-026-02110-0 (2026). [67] H. Zhang, C. Liu, Z. Wang, X. Wei, S. Liu, N. Zheng, B. Shao, and T.-Y. Liu, Self-consistency training for density-functional-theory Hamiltonian prediction, in ICML, Proceedings of Machine Learning Research, Vol. 235, edited by R. Salakhutdinov, Z. Kolter, K. Heller, A. Weller, N. Oliver, J. Scarlett, and F. Berkenkamp (PMLR, 2024) p. 59329–59357. [68] Z. Wang, C. Liu, N. Zou, H. Zhang, X. Wei, L. Huang, L. Wu, and B. Shao, Infusing self-consistency into density functional theory hamiltonian prediction via deep equi- librium models, in NeurIPS, NIPS ’24 (Curran Associates Inc., Red Hook, NY, USA, 2024). [69] M. Cranmer, S. Greydanus, S. Hoyer, P. Battaglia, D. Spergel, and S. Ho, Lagrangian neural networks, arXiv preprint arXiv:2003.04630 (2020). [70] S. Greydanus, M. Dzamba, and J. Yosinski, Hamiltonian neural networks, NeurIPS 32 (2019). [71] J. Sohl-Dickstein, E. Weiss, N. Maheswaranathan, and S. Ganguli, Deep unsupervised learning using nonequi- librium thermodynamics, in ICML (PMLR, 2015) p. 2256–2265. [72] Y. Song, J. Sohl-Dickstein, D. P. Kingma, A. Kumar, S. Ermon, and B. Poole, Score-based generative modeling through stochastic differential equations, in ICLR (2021). [73] A. Q. Nichol and P. Dhariwal, Improved denoising dif- fusion probabilistic models, in ICML (PMLR, 2021) p. 8162–8171. [74] W. Grathwohl, R. T. Chen, J. Bettencourt, I. Sutskever, and D. Duvenaud, FFJORD: Free-form continuous dy- namics for scalable reversible generative models, in ICLR (2019). [75] G. Papamakarios, E. Nalisnick, D. J. Rezende, S. Mo- hamed, and B. Lakshminarayanan, Normalizing flows for probabilistic modeling and inference, J. Mach. Learn. Res. 22, 1 (2021). [76] B. Zhang and R. Sennrich, Root mean square layer nor- malization, in NeurIPS, Vol. 32 (2019). [77] Y.-L. Liao, A. J. Hoffman, S. C. Shen, A. Duval, S. W. Norwood, and T. Smidt, EquiformerV3: Scaling effi- cient, expressive, and general SE(3)-equivariant graph attention transformers, arXiv preprint arXiv:2604.09130 (2026). [78] A. Girard, A fast ‘Monte-Carlo cross-validation’ proce- dure for large least squares problems with noisy data, Numer. Math. 56, 1–23 (1989). [79] I. Loshchilov and F. Hutter, Decoupled weight decay reg- ularization, in ICLR (2019). [80] K. Sch ̈utt, P. Kessel, M. Gastegger, K. A. Nicoli, A. Tkatchenko, and K.-R. M ̈uller, SchNetPack: A deep learning toolbox for atomistic systems, J. Chem. Theory Comput. 15, 448 (2018). [81] K. T. Sch ̈utt, S. S. Hessmann, N. W. Gebauer, J. Lederer, and M. Gastegger, SchNetPack 2.0: A neural network toolbox for atomistic machine learning, J. Chem. Phys. 158 (2023). [82] D. Hendrycks and K. Gimpel, Gaussian error linear units (GELUs), arXiv preprint arXiv:1606.08415 (2016). [83] Q. Li, Z. Han, and X.-M. Wu, Deeper insights into graph convolutional networks for semi-supervised learning, in Proceedings of the AAAI Conference on Artificial Intel- ligence, Vol. 32 (2018). [84] K. Oono and T. Suzuki, Graph neural networks expo- nentially lose expressive power for node classification, in ICLR (2020). [85] A. Griewank, On automatic differentiation and algorith- mic linearization, Pesq. Oper. 34, 621 (2014). [86] E. Hairer, G. Wanner, and S. P. Nørsett, Solving ordi- nary differential equations I: Nonstiff problems (Springer, 1993). [87] M. I. Budyko, The effect of solar radiation variations on the climate of the Earth, Tellus 21, 611 (1969). [88] K. McGuffie and A. Henderson-Sellers, A Climate Mod- elling Primer, 3rd ed. (John Wiley & Sons, 2005). [89] A. Szabo and N. S. Ostlund, Modern Quantum Chem- istry: Introduction to Advanced Electronic Structure The- ory (McGraw-Hill, 1989). [90] R. G. Parr and W. Yang, Density-Functional Theory of Atoms and Molecules, International Series of Mono- graphs on Chemistry, Vol. 16 (Oxford University Press, New York, 1989). [91] S. Greydanus, M. Dzamba, and J. Yosinski, Hamilto- nian neural networks, in NeurIPS , Vol. 32, edited by H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alch ́e- Buc, E. Fox, and R. Garnett (Curran Associates, Inc., 2019) p. 15379–15389. S1 Appendix S1: Long-range effects and generalization in cumulene molecules Energy increase through torsion (Relaxed at 0°) (Relaxed at 90°) (d) Force MAE (in eV/A, trained on - ) (c) Energy MAE (in eV, trained on - ) (b) Extrapolation to longer chain lengths (trained on - ) (a) Training & Testing on all chain lengths Implicit0.00580.00430.00360.00370.00590.03460.04270.0514 5Explicit0.00650.00380.00320.00380.00280.05030.06390.0502 1.75Explicit0.00670.00620.00460.00410.07050.11620.10280.1161 Implicit0.01090.00740.00720.00600.00560.02000.03110.0473 12Explicit0.00510.00350.00340.00340.00300.07380.10360.1138 Explicit-60.00770.00570.00430.00500.00510.03300.05190.0553 Implicit0.00590.00370.00440.00330.00290.09040.08800.1010 CHCHCHCHCHCHCHCHCutoffType Implicit0.00350.00350.00210.00320.00300.09480.23270.3788 5Explicit0.00550.00280.00210.00450.00500.18710.30770.3171 1.75Explicit0.00330.00330.00150.00250.27810.24060.23760.3635 Implicit0.01190.02150.01800.01400.01140.01380.05610.1567 12Explicit0.00090.00160.00460.00440.00670.20890.26220.4878 Explicit-60.00660.00530.00400.01030.01110.08820.18320.2070 Implicit0.00090.00070.00160.00270.00110.21360.29570.2383 CHCHCHCHCHCHCHCHCutoffType FIG. S1. Long-range generalization on the cumulene molecules C n H 4 (n = 2,..., 9). Each panel shows the molecular energy relative to the relaxed structure as a function of the torsion angle between the two terminal CH 2 groups; the black line is the GFN2-xTB label and colored lines are model predictions. Line colors denote the message-passing cutoff (1.75 ̊ A, 5 ̊ A and 12 ̊ A) and styles distinguish implicit (dashed) from three-layer explicit (dotted) models, with Explicit-6 (dash-dotted) a six-layer explicit baseline. (a) models trained on all eight molecules; almost all implicit and explicit variants reproduce the test curves, the explicit 1.75 ̊ A model being the notable exception (constant prediction on the longer chains). (b) extrapolation, with models trained only on C 2 H 4 to C 6 H 4 and evaluated on the larger, unseen chains C 7 H 4 to C 9 H 4 (bottom row). Implicit models generalize comparatively well to the unseen molecules, whereas explicit models fit them poorly and predict a qualitatively wrong energy decrease on the longest chains. (c & d) per-molecule and per-model mean absolute errors of energy (in eV) & force predictions (in eV/ ̊ A) in the extrapolation setting, the columns C 7 H 4 to C 9 H 4 are unseen during training. Green cells achieved chemical accuracy of 1 kcal/mol = 0.043 eV on the energy, or an analogous target force error of 1 kcal/mol/ ̊ A = 0.043 eV/ ̊ A [59]. Yellow cells are above 1× and red cells above 3× these target accuracies. Implicit modeling improves both energy and force profiles. An implicit model is equivalent to an explicit MLFF with an infinite stack of weight-sharing layers, and therefore has no finite horizon on how far apart two structures can interact in its prediction. In contrast, an explicit network can only propagate information over its ‘effective cutoff’, at most the number of layers times the neighborhood cutoff (3× 5 ̊ A). We test whether this ‘infinite’ effective cutoff, by communicating information across the whole system, improves predictions that depend on long-range effects, using the cumulene dataset [2]. The eight molecules C n H 4 (n = 2,..., 9) differ structurally only in the length n of the double-bonded carbon backbone (Fig. S1). Within each molecule’s test set only the torsion angle between the two terminal methylidene (CH 2 ) groups changes gradually, driving an increase in potential energy through a long-range electronic interaction. The longest interatomic distance – and thus the length scale over which this torsion must be communicated to predict the energy – grows with the chain, from 3.1 ̊ A on C 2 H 4 to 11.5 ̊ A on C 9 H 4 . Training on 4500 samples of each molecule S2 allows almost all implicit and explicit models to predict the test curve well (Fig. S1a). The exception is a three-layer explicit network with a 1.75 ̊ A cutoff. This cutoff is chosen to lie between the largest bonded (1.6 ̊ A) and smallest non- bonded (1.8 ̊ A) atom distance so that messages travel only along bonds. This short-range model fails once aggregating the methylene positions takes more than three steps: beyond five carbons its prediction is constant, near the per- molecule training mean. In contrast, an implicit model with the same 1.75 ̊ A cutoff predicts the energy near-perfectly at any backbone size, without oversmoothing [83, 84] even when up to ten propagation steps are required. We use the same dataset to demonstrate extrapolation capabilities of implicit models. By training only on cumulene molecules with up to six carbon atoms, we demonstrate that simply ‘turning MLFFs implicit’ without changing the underlying architecture improves their generalization to the larger chains of up to nine carbons. Dashed lines in the bottom row of Fig. S1b highlight that all implicit models fit the test set energy curves of unseen molecules comparatively well. Especially the default configuration with a message-passing cutoff of 5 ̊ A performs well with energy errors around 0.04 eV to 0.05 eV across the three unseen chains. Every implicit model performs better than its explicit counterpart (Fig. S1c). Implicit energy errors occur most at the difficult high torsion region where the ground-truth energy has a discontinuous cusp and a physical bond break might occur. In part energy errors can be corrected by a constant energy offset which might be admissible as no systems of such large size were seen during training. This shows up in the forces – the energy’s derivative, and the quantity that actually drives molecular dynamics – which the implicit models recover more faithfully than their explicit counterparts (Fig. S1d), even where the absolute energy is shifted. Explicit models, in contrast, fit C 7 H 4 poorly and give a qualitatively wrong prediction on C 8 H 4 and C 9 H 4 , modeling a decrease in potential energy while the ground-truth increases. These findings hold across the ablated cutoffs (5 ̊ A, 12 ̊ A and 1.75 ̊ A); notably the larger 12 ̊ A cutoff, which extends the explicit effective cutoff, hurts rather than helps generalization. The only visible direction of improvement is seen when increasing the layer count of the explicit model from three to six layers (Explicit-6, blue dash-dotted line), which moves closer to the ‘infinite-depth’ implicit modeling regime. The deep explicit model generalizes decently to chains with seven or eight carbons, but still predicts an incorrect energy decrease on the largest molecule. For an explicit model of any fixed depth one can construct a system whose electronic effects act over too long a range to be modeled. Implicit networks instead have an effective range that adapts to the geometry, and we observe no such breakdown across the chain lengths tested – suggesting that the inductive bias of implicit modeling improves generalization beyond what a larger explicit cutoff or depth provides. FIG. S2. Vibrational spectra for ethanol. The I-SO3net model predicts the 50 ps trajectory with, on average, 1.77 iterations per timestep. Dashed lines indicate peaks of ethanol’s vibrational spec- trum in the MD17 reference trajectory. Model Matched peaks Mean deviation (cm −1 ) I-SO3net109±8 1 layer1011±12 2 layer914±15 3 layer813±12 TABLE S1. Accuracy of spectral peaks. This table qualitatively compares the performance of ex- plicit SO3nets to an implicit SO3net on recreating physical observables. Appendix S2: Vibrational spectra Efficient implicit force fields, with less than 2 layer evaluations on average, maintain the ability of explicit models to reproduce physical observables like vibrational spectra, as visualized in Fig. S2 and surpass their accuracy in certain spectral peaks. Table S1 helps to analyze the spectra more quantitatively. Peaks are detected as all frequencies with a higher intensity than their surrounding 40 cm −1 sliding window. For a peak to ‘match’ its DFT reference, it must be within this window. The I-SO3net model captures the reference peaks most precisely, resulting in the highest number of matched peaks. Furthermore, it shows the lowest mean absolute deviation and variance from the DFT reference S3 peaks. While the explicit SO3net with one layer matches the same number of peaks as I-SO3net, it also exhibits many spurious oscillations, especially around the double peaks at 1400 cm −1 and 3000 cm −1 and the O-H peak at 3700 cm −1 . 1st step, 0 ̊ Torsion: Many iters. to converge 2nd step, 45 ̊ Torsion: Acceleration with fixed point reuse 3rd+ step, 90 ̊ Torsion: Quick convergence with linear extrapolation. h t * ( ) Reuse Linear extrapolation t 3 t 2 t 1 h(t 3 ) * h (t 2 ) (1) h (t 2 ) (0) h (0) (t 3 ) h (t 1 ) (0) h (t 1 ) (1) h(t 1 ) * h(t 2 ) * f f f f f (b) Continuity of fixed points along C 2 H 4 torsion trajectory (a) Spherical harmonic feature of C 2 H 4 evolving through solver iterations and warmstarts on three conformations FIG. S3. Supplement to Fig. 1(a) on the C 2 H 4 molecule.(a) Visualization of one spherical harmonics feature as a 3D density contour. The bottom row shows its geometry-invariant initialization from the atomic numbers and how successive iterations of f refine it to the fixed point. The previous fixed point is applied to the new geometry in the second row (following the black dashed arrow), warmstarting and accelerating the solver. Linear extrapolation from two preceding fixed points (the converging black dashed arrows, top left) warmstarts the third step at 90 ◦ torsion precisely. (b) Continuity of fixed-point trajectory: The fixed-point trajectory of C 2 H 4 projected onto the plane spanned by its first two principal components, overlaid on a contour plot of the potential energy surface sampled on a fine grid. 3rd+ MD step: Linear extrapolation 1 iteration 1st MD step: Many iterations to converge 2nd MD step: Reuse of fixed-point 2-5 iterations FIG. S4.Vector fields underlying the 3D streamplot Fig. 1(b).Three consecutive MD timesteps of I-PaiNN on azobenzene, showing the vector field (evaluated on a fine grid on the panel plane), warmstart initializations (red), and fixed- point iteration trajectories f (blue). The plane visualized in each panel is spanned by the three fixed points of the three steps. S4 Appendix S3: Supplement to Fig. 1a,b – Fixed-point trajectory details In the main text Fig. 1(a,b) sketches the continuity and our extrapolation on the trajectory of fixed points h ∗ (t) through time. To visualize continuity we choose a torsion trajectory of C 2 H 4 because the simplicity of the dynamics and the small size of the molecule allow for easy visualization. In the main text we stylize the fixed-point trajectory for reasons of space and style (preserving the core characteristics). Fig. S3 supplements the visualization of the exact C 2 H 4 fixed-point trajectory projected on the plane spanned by the two principal components. We also sample the energy readout of the plane in a fine grid revealing a contour plot of the potential energy surface (PES). The visualization is slightly flawed: it looks like torsion causes an initial dip in the energy prediction – it does not. The high-dimensional trajectory has monotonically increasing energy, only the two-dimensional slice on which the contour plot is visualized has an energy decrease in the place where the trajectory’s projection lands. Still, we can discern inner workings of the model: At around 45 ◦ the trajectory aligns with the first principal component, and with the energy gradient, which causes the steep increase in predicted energy starting from there on (see SI Fig. S1). The torsion changes on the C 2 H 4 are very small, making warmstarts and consecutive convergence too easy for an educational visualization. All vector fields, warmstarts, and iteration f shown in Fig. 1b therefore visualize real data from an I-PaiNN model on the azobenzene system. Fig. S4 provides a more detailed view of the three panels (which are stacked vertically in the main text Fig. 1(b)), including the vector field (blue, calculated on a fine grid on the plane) and the iterations f (red; in the second and third timestep they start in the middle of the panel because of the warmstarts). The plane visualized in all three panels (here and in the main text) is the plane spanned by the three fixed points of the three steps. Appendix S4: Supplement to Fig. 4 – Increased solver iterations at high potential energies T=100K 200K 300K 400K 500K 600K 1 3 5 6 4 2 Mean solver iterations No warmstart Fixed-point reuse Linear extrapolation Density of MD: I-PaiNN at 400K Density of dataset: DFT at 500K 80% 60% 40% 20% 0% Fraction of conformations (y-axis) of the quartile that falls into the iteration bucket (x-axis) (c) Iterations increase with temperature (Converging up to usual tolerance of 10 -2 ) (b) Fraction of conformations of an energy quartile (y-axis) that fall into range of iterations (x-axis) (a) Iterations as function of potential energy FIG. S5. Dependence of average solver iterations on potential energy and temperature on Ac-Ala3-NHMe (cf. Fig. 4). (a) mean iteration count (red, right axis) as a function of potential energy, overlaid on the energy density of the test MD at 400 K (blue) and the DFT training dataset at 500 K (green). (b) fraction of conformations in each energy quartile that fall into each iteration bucket. (c) Converging to a threshold of 10 −2 takes ∼5 iterations without warmstarts and ∼1 iterations with either type of warmstart, with only slight increases with the temperature. To differentiate the warmstarts better, Fig. 4d in the main text plots iterations to converge to a very fine threshold of 10 −3 . In the main text, Fig. 4 shows how the iteration count of implicit models (especially without warmstarts) increases visibly in regions of the Ramachandran plot that are known to be outside of the main energy minima of the potential energy surface. Fig. S5 supplements this analysis by focusing directly on the relationship of the potential energy of the system (or its instantaneous temperature) and the iterations that the implicit model invests to solve it. In Fig. S5a the distribution of potential energies sampled in the MD of the implicit models is shown in blue. We see that the MD22 trajectory from which the training set of the model is drawn samples higher energies (green): As MD22 is simulated at 500 K and our MD at 400 K, MD22 had more kinetic energy in the system, which (by equipartition) caused statistically higher potential energies, too. Notably, iterations, averaged over the 30 equiwidth bins of the plot, increase monotonically with potential energy except for noise caused by low sample counts at the edges of the blue distribution. In Fig. S5b we take a more detailed view. We categorize conformations into four quartiles based on their potential energy and color them accordingly. Next, we sort conformations by how many iterations the solvers took to predict their interatomic forces. Over 80 % of the conformations in the lowest energy quartile are solved within fewer than 3 iterations. For comparison, the overall average iteration count is 8.9, three times as high. In contrast, high potential energy conformations are strongly skewed to take higher iteration counts, suggesting that the model invests more compute in conformations of higher chemical diversity. S5 Fig. S5c plots the iteration counts needed to reach our usual threshold of 10 −2 in MD simulations run at different temperatures. High temperatures indicate high kinetic energy of the molecule, which is statistically linked (through equipartition) to increased potential energies. The correlation of iterations to potential energy is thus as expected: increases in temperature correlate with a minor increase in iteration counts. However, as Ac-Ala3-NHMe is very easy to solve using implicit models, the iteration counts of the warmstarted solvers stay very low, close to the minimum of one iteration per MD step at most temperatures. Fig. 4d in the main text plots convergence to a finer tolerance of 10 −3 where linear extrapolation distinguishes itself, maintaining low iteration counts even at high temperatures. Appendix S5: Accuracy tables Table S2 displays the test set accuracies of implicit and explicit MLFFs for each molecule and model type, separated into three subtables for the PaiNN, SO3net, and SchNet architectures. Only the force MAE for the implicit model using the equivariant merged layer norm (LN) [77] is shown in kcal/mol/ ̊ A. All other models on this architecture and molecule reference its performance and plot the percentual in- or decrease of their force MAE in shades from red to green. Implicit models using the simpler equivariant UL normalization slightly underperform LN on close to all tasks, likely because of the lower expressiveness of the UL norm. We train explicit models up to three layers on SO3net and SchNet, and up to six and eight layers with PaiNN on molecules of the MD22 and MD17 datasets, respectively. Explicit models that share parameters across all layers are plotted for PaiNN on MD17 in a gray, italic font and seem to mirror their explicit counterparts’ performance closely. Notably, the implicit model exceeds or matches the accuracy of its explicit counterparts up to five or six layers, while requiring only a fraction of their memory and compute using warmstarts. The only exception forms benzene, a very stable molecule with little variation in nuclear positions and a high number of symmetries in its circular structure: It is solved best by a two-layer explicit PaiNN but degrades again on high depths, possibly because of the oversmoothing effect [83, 84]. Explicit models seem to scale worse to high depths on larger molecules. We speculate that this is at least in part caused by training dynamics getting increasingly unstable with depth, an issue that we could not fix reliably even when restarting training with different random seeds. Implicit models with their conceptually infinite depth avoid this issue, possibly by training and regularizing the model for a more restrictive fixed-point representation, or by adding input-injection and normalization to the architectures. Table S3 summarizes the detailed tables by taking the mean over force MAEs of all molecules in a dataset and computing the relative in- or decrease over these aggregated scores. Across the three architectures, explicit models on the larger molecules of MD22 often catch up to the implicit accuracy with about three layers, bringing the necessary compute down to roughly the level of an implicit warmstarted force field. The memory consumption of implicit models in this setting is precisely one third of a three-layer explicit network of the same architecture, driving significant savings especially on large molecular systems, cf. Fig. S6. S6 TABLE S2. Force accuracy per molecule, with three sub-tables focusing on the PaiNN, SO3net, and SchNet architectures. In all tables the force MAE of the implicit model using the equivariant layer norm (LN ) is shown in kcal/mol/ ̊ A on a blue background. All other models show their increase (red) or decrease (green shades) in force MAE relative to the LN baseline. TABLE S3. Force accuracy per architecture and dataset, averaging the per-molecule force MAE shown above and using the same color scheme. The implicit LN is shown in kcal/mol/ ̊ A; other models note their in- or decrease relative to LN. S7 FIG. S6. Peak memory required for calculating forces as an energy derivative in implicit and explicit models of one, three, and five layers. Implicit models require roughly 1/K of the memory of an explicit model with K layers, as only the last application of f of the message-passing layer near the fixed point h ∗ = f (h ∗ , x) is required for implicit differentiation. The x-axis denotes the tested systems and their count of atoms (a) and message-passing edges (e) in the 5 ̊ A cutoff, on a conformation that’s randomly drawn from the datasets. Appendix S6: Memory consumption Fig. S6 documents the memory consumption of an implicit model and three explicit models with 1, 3, and 5 message-passing layers on the SchNet, PaiNN, and SO3net architectures. The x-axis denotes a range of molecules for which the measurement was performed, as well as the number of atoms and edges within a 5 ̊ A cutoff that drive the memory consumption. We use implicit differentiation to calculate forces as the derivative of the energy measured at the fixed-point representation h ∗ , which requires only local derivative information around the fixed point and not the entire fixed- point solver’s path, cf. Sec. IV B. In practice, this means that only the intermediate activations within the last iteration of the solver need to be kept in memory, no matter how many iterations are necessary to find the fixed point. This contrasts with explicit force fields which need to store and backpropagate through all of their message-passing layers (irrespective of whether the layers share learned parameters). In the largest measured system, the double-walled nanotube, this brings the memory consumption of a three- or five-layer explicit SO3net during force inference to 6.8 GB or 10.8 GB, respectively, while the implicit model requires only 2.8 GB. These measurements include the fixed-point history required for the linear warmstarts of implicit models. Their contribution and other factors such as the reduced parameter storage are minimal compared to the storage of message-passing activations. S8 Algorithm 1: ImplicitMLFF Implicit energy and force evaluation Require: Atomic positions R, atomic numbers Z, initial state h and adjoint ̄u, tolerances ε h ,ε ̄u Ensure: Energy E, forces F, and converged h, ̄u ▷ Forward (Calculate fixpoint h and E) 1: repeat 2:h prev ← h 3:h← f (h prev , R, Z) 4: until ∥h− h prev ∥≤ ε h 5: E ← f E (h) ▷ Backward (derive F implicitly at fixed point h) 6: w←∇ h E 7: repeat 8: ̄u prev ← ̄u 9:u← ̄u + w 10:( ̄u,−F)← u ⊤ ∂f/∂(h prev , R) 11: until ∥ ̄u− ̄u prev ∥≤ ε ̄u 12: return (E, F, h, ̄u) Algorithm 2: Warmstarts (Using a two-step history) State: Previous fixed points (h ∗ −∆t , h ∗ −2∆t ), ( ̄u ∗ −∆t , ̄u ∗ −2∆t ) Require: Atomic positions R, atomic numbers Z, Ensure: Forces F ▷ Warm-start via linear extrapolation 1: h← 2 h ∗ −∆t − h ∗ −2∆t 2: ̄u← 2 ̄u ∗ −∆t − ̄u ∗ −2∆t ▷ Implicit fixed point + differentiation 3: (E, F, h, ̄u)← ImplicitMLFF(R, Z, h, ̄u) ▷ Commit for next call 4: (h ∗ −2∆t , h ∗ −∆t )← (h ∗ −∆t , h) 5: ( ̄u ∗ −2∆t , ̄u ∗ −∆t )← ( ̄u ∗ −∆t , ̄u) 6: return F Algorithm 1 stores only the computation graph of the last iteration (h prev , R) f 7−→ h. The backward solver evaluates the VJP of this map wrt. both its inputs to solve equations 4 and 3 for F and ̄u ∗ in parallel. Algorithm 2 stores a rolling history of the previous two timesteps fixed points for h ∗ and ̄u ∗ and warmstarts Algorithm 1 on the next timestep with their linear extrapolation. Appendix S7: Algorithm of Implicit Conservative Force Fields Algorithm 1 summarizes the entire implicit energy-conservative MLFF. After receiving warmstarts for h and ̄u from Algorithm 2, it iteratively solves for both h ∗ to obtain the system’s potential energy E, and then for u ∗ to derive energy-conserving forces F =∇ R E. Notably, it decomposes the self-consistency equation (4) of the fixed-point adjoint u ∗ = ̄u ∗ + w into two parts, an intermediate variable ̄u ∗ := ( ∂f ∂h ∗ ) ⊤ u ∗ and the energy-head adjoint w := ∂f E ∂h ∗ . We capture the intermediate ̄u ∗ on preceding MD steps and extrapolate from them analogously to equation (6) for a guess of the current MD step ̄u≈ ̄u ∗ . This allows Line 10 to estimate the fixed-point adjoint u← ̄u + w in a single vector addition and to compute the improved estimates for ̄u ← u ∗ ∂f ∂h ∗ and F ← u ∗ ∂f ∂R in parallel, with a Vector Jacobian Product (VJP) of the stored last execution of f . Within a few iterations we obtain accurate forces F, as well as a nearly converged ̄u ≊ ̄u ∗ from which the next MD step is warmstarted. Compared to serially deriving u ∗ in equation (4) and then F in equation (3), parallel derivation using ̄u saves one VJP per MD step. Since implicit force fields often require only one forward and backward iteration each to converge sufficiently after warmstarts, and f and its VJP have a similar time complexity [85], saving one VJP can reduce the per-step cost of the entire implicit MLFF by up to 33 %. Finally, if f involves any computation regarding solely Z or R and not the iterate h, they should be computed once into intermediate embeddings, say h Z and h R , and fed into the remaining f at every iteration without recomputing them. Similarly during force derivation, each backward iteration should yield the cotangent dE/dh R , but only once convergence is detected should we backpropagate through the embedding from h R to R to derive atomic forces. This is analogous to cross-layer optimization in explicit networks. Appendix S8: Extrapolation from three or more fixed points The linear predictor in equation (6) can be extended to higher polynomial degree at negligible additional cost. A Taylor expansion of the equilibrium trajectory h ∗ (t) around time t gives the degree-k approximation T k (t + ∆t) = k X j=0 1 j! d j h ∗ dt j t ∆t j .(S81) The key idea behind Adams–Bashforth (AB) methods [86] is that k samples of a function’s first derivative at con- secutive time steps – computed from k+1 stored values, i.e. k predecessors – implicitly encode all derivatives up to order k−1, via the finite differences between samples. Rather than computing the higher derivatives in equation (S81) S9 explicitly, AB fits a degree-(k−1) polynomial interpolant through k first-order derivative samples and integrates it one step forward, yielding the prediction: h ∗ (t+∆t) = h ∗ (t) + ∆t k−1 X j=0 b j φ (t−j∆t) .(S82) Here,φ (t−j∆t) denotes the first time derivative of h ∗ at time t−j∆t. The coefficients b j are fixed rational numbers determined solely by the order k (Table S4); no matrix inversions or least-squares fits are required. TABLE S4. Adams–Bashforth coefficients b j in equation (S82). For j≥k, b j = 0. k b j for j<k 1 1 2 3 2 , − 1 2 3 23 12 , − 16 12 , 5 12 4 55 24 , − 59 24 , 37 24 , − 9 24 In the classical ODE setting, the derivativesφ (t−j∆t) are available directly from the evaluation of the ODE’s right- hand side. In our setting, the solver returns fixed-point values h ∗ (t) , h ∗ (t−∆t) ,... rather than time derivatives. We estimate the required first-order time derivatives from finite differences of the stored values. For the most recent point, only a backward difference can be computed, as h ∗ (t+∆t) is not yet available: φ (t) ≈ h ∗ (t) − h ∗ (t−∆t) ∆t .(S83) For all earlier points t− i∆t with i≥ 1, a central difference can be formed: φ (t−i∆t) ≈ h ∗ (t−(i−1)∆t) − h ∗ (t−(i+1)∆t) 2 ∆t .(S84) Substituting these derivative estimates into equation (S82) gives the degree-k warm-start for h ∗ , and analogously for ̄u. At k=1, this reduces to the linear predictor of equation (6). Higher degrees require storing k+1 previous fixed points but no additional model evaluations, making the computational overhead negligible. Appendix S9: Setup for iteration count measurements on aspirin To measure average iteration counts during MD on aspirin in Fig. 2, 10×5 ps NV E dynamics are simulated on the aspirin system at 500 K for each tested model and warmstart configuration. We set up all system geometries and momenta to be representative and have a broad coverage of the system ensemble. For this, we first sample 10 different random conformations of the MD17’s long aspirin trajectory and draw their atomic momenta from a normal distribution with a mean corresponding to the target temperature of 500 K. Finally, we apply three short equilibration phases, running Langevin dynamics thermostats with coupling strengths t LTC ∈ 50 fs, 200 fs and 1000 fs for t LTC each. This ensures that the total energy of the system corresponds to the target temperature, and momenta of every atom are pointing in a physically likely direction at the start of the measurement. We find these starting positions and momenta to be representative of the model ensemble, as the gathered statistics are consistent across the 10 replicas and across the length of the MD simulations. The fixed points h ∗ and ̄u ∗ are converged up to the required tolerance of 10 −2 during the simulations, which leads to the visualization of residuals in Fig. 2 (C.I.) to stop soon after crossing the red dotted threshold. Appendix S10: Setup for stability measurements on Ac-Ala3-NHMe For the stability experiments in Fig. 3B–C that focus on non-conservative energy drift, we measure the average temperature increase (per fs) across 15 simulations over 50 ps on the Ac-Ala3-NHMe molecule. We prepare each trajectory by relaxing a randomly drawn conformation from MD22, re-adding kinetic energy corresponding to a temperature of 600 K, and bringing it into thermal equilibrium at about 300 K using a short 5 ps NVE simulation at a fine solver tolerance of 10 −4 . S10 FIG. S7. Full illustration of an implicit model with a scalar-valued fixed point: The magnetic Ising model (mean-field Curie-Weiss) with h ∗ as the average magnetization, coupling strength 1 2 , and a varying external magnetic field x. The explicit function h (k+1) = tanh( 1 2 h (k) +x) – analogous to f in the main text – is plotted in blue. Where its input is equal to its output h (k+1) =h (k) (on the gray diagonal), it induces the implicit map x7→ h ∗ (x), which is visualized as the red dashed intersection. No analytical expression for h ∗ (x) exists, yet the (implicit) function varies smoothly with its input parameter x. The fixed point h ∗ (x) carries physical meaning as it is the average magnetization of the spins and has a well-defined derivative, the magnetic susceptibility χ = d dx h ∗ (x). This is analogous to the implicit energy prediction E and implicit derivation of forces F in equations (1) and (2). Because of its smooth shape, warmstarting with polynomial extrapolation like in equation (6) from neighbouring points of the red fixed-point trajectory h ∗ (x) can accelerate convergence. Three non-warmstarted solver trajectories (orange, olive, and wine) visualize how the convergence speed depends on the flatness of the underlying explicit function. Jacobian regularization as introduced in Sec. IV C aims to nudge learned parameters to regularize the spectrum of the blue plane, producing exactly such flat regions and fast convergence. Appendix S11: Implicit model illustration in 1D – The magnetic Ising model Fig. S7 visualizes the full mechanics of a one-dimensional implicit model in one plot. This includes most of the concepts discussed in the main text: Solver iteration, the implicit function, its smooth dependence on inputs, and a well-behaved spectrum as desiderata. In general terms, an implicit function is any function whose input depends on its output. As such, most phys- ical laws can be described in terms of implicit functions – often observations made in nature are produced as the balance, or equilibration, of opposing and interacting effects. Examples include the equilibrium temperature of a planet, where absorbed solar flux and emitted thermal radiation form a self-consistent equation for the ‘fixed-point’ temperature [87, 88]. In quantum chemistry, the Hartree-Fock and Kohn-Sham density functional theory equations are fixed-point problems where the molecular orbitals must diagonalize an operator that is itself constructed from those same orbitals [89, 90]. We analyze the Ising model, a simple, well-studied implicit model of the interaction of magnetic spins. Specifically, in the mean-field Curie-Weiss formulation the whole system reduces to a one-dimensional (scalar) fixed-point equation. Adopting our notation, it models the average magnetic dipole moment h ∗ for all the many spins, under the condition that they are only influenced by the mean-field of all the other spins h ∗ in the system and an external magnetic field x. This self-consistency and input dependence mirror the structure of I-MLFFs precisely, and just as the fixed point of an I-MLFF carries physical meaning in the form of the total energy, the fixed point of the Ising model corresponds to the average magnetization per spin. Similarly, the implicit derivative of the fixed point equates to the magnetic susceptibility of the system. These and further analogies, such as warmstarting and the effect of the shape of the underlying explicit function, are described in the caption of Fig. S7. S11 Appendix S12: Implicit Hamiltonian Neural Networks (I-HNN) h (0) u (0) h h (k) u (k) u * h * u h * Energy Energy Iterations (smoothed over time) Iterations (smoothed over time) (a) Algorithmic flowchart of I-HNN (b) Ground-truth trajectory (c) I-HNN trajectory (d) Ground-truth energy (e) I-HNN energy & iterations Time t 3 t 2 t 1 FIG. S8. Architecture and performance of an I-HNN on the gravitational two-body problem. (a) a schematic visualization of the I-HNN similar to 1. (b) the ground truth trajectory and the kinetic energy at each point. (c) the trajectory predicted by the I-HNN almost perfectly matches the ground truth and the iteration count correlates pointwise with the kinetic energy. (d) and (e) show the energy dynamics of the ground truth and the dynamics predicted by the I-HNN respectively and the iterations needed for the I-HNN to converge. For illustration purposes, the iteration count was smoothed symmetrically over 30 time steps with a Gaussian convolution. Implicit MLFFs have unique benefits in compute and memory costs as the sequential nature of predictions permits warmstarting fixed-point solvers, and the implicit derivation requires backpropagating only a single layer instead of an entire stack of neural network layers. The same techniques generalize naturally to diverse applications in physical simulation and image processing, cf. Sec. I, which we demonstrate here on a small Implicit Hamiltonian Neural Network (I-HNN) trained on the gravitational two-body task [91]. In contrast to HNNs, MLFFs as discussed in this study predict the potential energy E only from the particle positions R (and some constant inputs), and derive the forces F =−∇ R E. A separate integrator like Velocity-Verlet integrates the acting forces into the system momenta M and drives the system positions. This split is not possible for more general Hamiltonian systems, like for a charged particle in a magnetic field where both the positions and momenta influence the potential energy E(R, M) (as well as the kinetic energy K(R, M)). Here, simulation does not require predicting the potential energy but the more general Hamiltonian H = E + K and integrating the change in both of the system’s inputs simultaneously, [ ̇ R, ̇ M] = [∇ M H,∇ R H]. Our modeling techniques extend naturally to this domain. On the two-body gravity simulation, we adopt the small Multi-Layer Perceptron used as an HNN in [91] into an implicit structure analogous to I-MLFFs visualized in Fig. S8 (a): First, a linear layer embeds both inputs into h R,M . Then an affine layer (f ) with a pointwise tanh non-linearity is iterated h (k+1) = f (h (k) + h R,M ) until the relative residual is smaller than 10 −2 . The Hamiltonian H is predicted with a linear layer from h ∗ . Implicit differentiation analogous to SI Sec. S7 gives the instantaneous changes [ ̇ R, ̇ M] = [∇ M H,∇ R H], which are used to evolve the system in time. Using linear extrapolation from two previous time steps to warmstart both the forward and backward solver, the I-HNN needs on average 2.1 layer calls to converge, while maintaining the total energy of the system. Fig. S8 (c) and (e) illustrate how these iterations are distributed across the trajectory: More compute is invested in the ‘difficult’ regions of high kinetic energy (in which both implicit and explicit models tend to produce the most significant energy-conservation errors).