Faster Molecular Dynamics with Neural Network Potentials via Distilled Multiple Time-Stepping and Non-Conservative Forces
Authors: Nicolaï Gouraud, Côme Cattin, Thomas Plé, Olivier Adjoua, Louis Lagardère, Jean-Philip Piquemal
Organizations: Qubit Pharmaceuticals, Advanced Research Department, 75014 Paris, France · Sorbonne Université, Laboratoire de Chimie Théorique, UMR 7616 CNRS, 75005 Paris, France
Following our previous work (J. Phys. Chem. Lett., 2026, 17, 5, 1288-1295), we propose the DMTS-NC approach, a distilled multi-time-step (DMTS) strategy using non-conservative (NC) forces to further accelerate atomistic molecular dynamics simulations using foundation neural network models such as FeNNix-Bio1. There, a dual-level reversible reference system propagator algorithm (RESPA) formalism couples a target accurate conservative potential to a simplified distilled representation optimized for the production of non-conservative forces. Despite being non-conservative, the distilled architecture is designed to enforce key physical priors, such as equivariance under rotation and cancellation of atomic force components. These choices facilitate the distillation process and therefore improve drastically the robustness of simulation, significantly limiting abnormal discrepancies between the two models, thus achieving excellent agreement with the forces data. Overall, the DMTS-NC scheme is found to be more stable and efficient than its conservative counterpart with additional speedups reaching 15-30% over DMTS. Requiring no fine-tuning steps, it is easier to implement and can be pushed to the limit of the systems physical resonances to maintain accuracy while providing maximum efficiency. We obtain additional speedup by combining hydrogen mass repartitioning (HMR), High Hydrogen Friction (HHF) to further extended the largest timestep up to 10fs of our schemes while conserving stability and accuracy. As for DMTS, DMTS-NC is applicable to any neural network potential and can be applied to approaches that are computationally heavier than FeNNix-Bio1. We show a proof of principle applying the approach to the distillation of MACE-OFF23 with consequent speedups ranging from 3.66 to 5.64 compared to single timestep.
Figures & tables
Number of atoms
Numerical
Integrator
648
4800
12000
instabilities
STS
36.89
9.08
3.63
No
DMTS 1-5fs
69.54
30.17
11.95
No
DMTS 1-6fs
80.62
34.34
NaN
Data holes
DMTS 1-7fs
NaN
NaN
NaN
—
DMTS-NC 1-5fs
91.13
32.69
12.39
No
Table 1: Simulation performances for bulk water systems, expressed in nanoseconds per day. In DMTS-NC, the stability limit is reached around an external time step of Δ=6.5 fs, after what resonances start to appear. In contrast, the main limitation the with the (generic) conservative model is the presence of “data holes” already at Δ=6 fs.
Figure 1: Sampling properties of bulk water simulations. All data was collected from 2 nanosecond simulations, with one frame saved per picosecond. Figure (a) shows that radial distributions are not affected by the use of DMTS-NC compared to DMTS or STS, Figure (b) shows that temperatures fit the equilibrium theoretical distribution, and samples of potential energies, represented on a boxplot (Figure (c)) and a quantile-quantile plot (Figure (d)), highlight the low-bias identical-variance deviation of the distributions from STS to DMTS-NC, as expected for a multi-time-step method in a stable regime.
Speed (ns/day)
Numerical
Integrator
Phenol-lyso.
DHFR
instabilities
STS
2.55
1.85
No
DMTS 1.75-3.5
6.59
4.46
No
DMTS 1-4
NaN
NaN
—
Fine-tuned DMTS 2-4
7.45
4.97
No
DMTS-NC 1-4
7.20
4.89
No
Table 2: Speed performances for the two ligand-protein systems. Similarly to bulk water simulations, DMTS-NC allows to use a larger external time step than conservative DMTS, without relying on active learning.
System
Integrator
Temperature
Potential energy
(K)
(kcal/mol/atom)
Water boxes
216 mol.
STS
299.84 9.57
−8.049 0.042
DMTS-NC
299.87 9.56
−8.038 0.043
1600 mol.
STS
300.09 3.60
−8.047 0.017
DMTS-NC
299.93 3.54
−8.037 0.016
4000 mol.
STS
300.61 2.93
−8.043 0.014
Table 3: Average temperatures and potential energies ( ± empirical standard deviations). DMTS-NC simulations were performed at the highest stable time steps that we found, i.e. δ=1.25 fs– Δ=5 fs for the protein systems, and δ=1.2 fs– Δ=6 fs for bulk water.
Figure 2: Time evolution of the protein backbone RMSD (Figure (a)) and the ligand’s Distance to Bound Configuraiton (Figure (b)) during a 20ns simulation of the solvated lysosyme-phenol. The DMTS-NC integrator uses δ=1.25 fs– Δ=5 fs time steps.
Figure 3: Hydration free energy of 44 small molecules in NPT. Error bars are computed from three independent lambda-ABF runs for each integrator and each molecule.
Water boxes
Solvated prot.
System size / Name
1600 mol.
4000 mol.
Phenol-Lys.
DHFR
Speed (ns/day)
STS 1 fs
9.08
3.63
2.55
1.85
DMTS-NC
36.2
13.7
8.42
5.57
HHF
47.9
17.2
9.09
6.38
HHF Accel. w.r.t.
DMTS-NC
32.61 %
25.92 %
7.96 %
14.57 %
STS
5.28
4.74
3.56
3.45
Table 4: HHF speed performances and acceleration compared to DMTS-NC and STS.
System
Time steps
Temperature
Pot. Energy (kcal/mol/atom)
(fs)
(K)
DMTS-NC
STS
Water boxes
1600 mol.
1–10
300 4
−8.02 0.02
−8.05 0.02
4000 mol.
1–10
300 2
−8.02 0.01
−8.04 0.01
Solvated prot.
Phenol-Lys.
1–6
297 2
−7.87 0.01
−7.87 0.01
DHFR
1–7
298 2
−7.81 0.01
−7.82 0.01
Table 5: HHF average temperatures and potential energies.
Water boxes
Solvated prot.
System size / Name
216 mol.
1600 mol.
4000 mol.
Phenol-Lys.
DHFR
Speed (ns/day)
STS
7.40
1.04
0.38
0.31
0.22
DMTS-NC
37.59
5.88
2.07
1.14
1.00
Acceleration
5.08
5.64
5.45
3.66
4.62
Table 6: DMTS-NC performances with MACE-OFF23(S).
Figure 4: Hydrogen velocity auto-correlation spectra for various external time steps, in a system of 1600 water molecules. From left to right, the first peak corresponds to molecule libration and slow fluctuations in the hydrogen bond network, the second peak to angle bending oscillations and the third and fourth peaks to bond stretching motions. The ones on the inset graph are non-physical artifacts. As the external time step increases, these peaks shift to the left, eventually coupling with the bond-oscillation peak at the stability limit.
Figure 5: Mean-square displacement (Figure (a)) and diffusion coefficients (Figure (b)) for various time steps and thermostats.
Integrator
Diffusion coeff. (Å 2 /ps)
Loss w.r.t. STS (%)
STS
0.1672 0.0017
—
STS-HMR
0.1472 0.0012
12.0
DMTS 1-5
0.1525 0.0007
8.8
DMTS-NC 1-6
0.1461 0.0005
12.6
HHF 1-10 Langevin
0.1016 0.0003
39.2
HHF 1-10 FFL
0.1126 0.0007
32.6
Table 7: Diffusion coefficients for various MTS strategies. The loss represents the relative decrease in diffusion compared to the STS reference.
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.
Johannes Maeß, Leon Werner, J. Thorben Frank +5
BIFOLD, Berlin, Germany. · Machine Learning Group, TU Berlin, Berlin, Germany. · Department of Artificial Intelligence, Korea University, Seoul, South Korea. +1
Accurate prediction of energy and forces for 3D molecular systems is one of fundamental challenges at the core of AI for Science applications. Many powerful and data-efficient neural networks predict molecular energies and forces from single atomic configurations. However, one crucial aspect of the data generation process is rarely considered while learning these models i.e. Molecular Dynamics (MD) simulation. MD simulations generate time-ordered trajectories of atomic positions that fluctuate in energy and explore regions of the potential energy surface (e.g., under standard NVE/NVT ensembles), rather than being constructed to steadily lower the potential energy toward a minimum as in geometry relaxations. This work explores a novel way to leverage MD data, when available, to improve the performance of such predictors. We introduce a novel training strategy called FRAMES, that use an auxiliary loss function for exploiting the temporal relationships within MD trajectories. Counter-intuitively, on two atomistic benchmarks and a synthetic system we observe that minimal temporal information, captured by pairs of just two consecutive frames, is often sufficient to obtain the best performance, while adding longer trajectory sequences can introduce redundancy and degrade performance. On the widely used MD17 and ISO17 benchmarks, FRAMES significantly outperforms its Equiformer baseline, achieving highly competitive results in both energy and force accuracy. Our work not only presents a novel training strategy which improves the accuracy of the model, but also provides evidence that for distilling physical priors of atomic systems, more temporal data is not always better.
Ali Mollahosseini, Mohammed Haroon Dupty, Wee Sun Lee
School of Electrical and Computer Engineering, University of Tehran, Tehran, Iran · School of Computing, National University of Singapore, Singapore
GROMACS is a de-facto standard for classical Molecular Dynamics (MD). The rise of AI-driven interatomic potentials that pursue near-quantum accuracy at MD throughput now poses a significant challenge: embedding neural-network inference into multi-GPU simulations retaining high-performance. In this work, we integrate the MLIP framework DeePMD-kit into GROMACS, enabling domain-decomposed, GPU-accelerated inference across multi-node systems. We extend the GROMACS NNPot interface with a DeePMD backend, and we introduce a domain decomposition layer decoupled from the main simulation. The inference is executed concurrently on all processes, with two MPI collectives used each step to broadcast coordinates and to aggregate and redistribute forces. We train an in-house DPA-1 model (1.6 M parameters) on a dataset of solvated protein fragments. We validate the implementation on a small protein system, then we benchmark the GROMACS-DeePMD integration with a 15,668 atom protein on NVIDIA A100 and AMD MI250x GPUs up to 32 devices. Strong-scaling efficiency reaches 66% at 16 devices and 40% at 32; weak-scaling efficiency is 80% to 16 devices and reaches 48% (MI250x) and 40% (A100) at 32 devices. Profiling with the ROCm System profiler shows that >90% of the wall time is spent in DeePMD inference, while MPI collectives contribute <10%, primarily since they act as a global synchronization point. The principal bottlenecks are the irreducible ghost-atom cost set by the cutoff radius, confirmed by a simple throughput model, and load imbalance across ranks. These results demonstrate that production MD with near ab initio fidelity is feasible at scale in GROMACS.
Luca Pennati, Andong Hu, Ivy Peng +2
KTH Royal Institute of Technology, Stockholm, Sweden